text:
from typing import Tuple, Optional, Union
import torch
import numpy as np
import SimpleITK as sitk
from ....torchio import DATA, AFFINE
from ....data.subject import Subject
from .. import RandomTransform

class RandomSpike(RandomTransform):
    r"""Add random MRI spike artifacts.

    Args:
        num_spikes: Number of spikes :math:`n` presnet in k-space.
            If a tuple :math:`(a, b)` is provided, then
            :math:`n \sim \mathcal{U}(a, b) \cap \mathbb{N}`.
            Larger values generate more distorted images.
        intensity: Ratio :math:`r` between the spike intensity and the maximum
            of the spectrum.
            Larger values generate more distorted images.
        p: Probability that this transform will be applied.
        seed: See :py:class:`~torchio.transforms.augmentation.RandomTransform`.

    .. note:: The execution time of this transform does not depend on the
        number of spikes.
    """
    def __init__(
            self,
            num_spikes: Union[int, Tuple[int, int]] = 1,
            intensity: Union[float, Tuple[float, float]] = (1, 3),
            p: float = 1,
            seed: Optional[int] = None,
            ):
        super().__init__(p=p, seed=seed)
        self.intensity_range = self.parse_range(
            intensity, 'intensity_range')
        if isinstance(num_spikes, int):
            self.num_spikes_range = num_spikes, num_spikes
        else:
            self.num_spikes_range = num_spikes

    def apply_transform(self, sample: Subject) -> dict:
        random_parameters_images_dict = {}
        for image_name, image_dict in sample.get_images_dict().items():
            params = self.get_params(
                self.num_spikes_range,
                self.intensity_range,
            )
            spikes_positions_param, intensity_param = params
            random_parameters_dict = {
                'intensity': intensity_param,
                'spikes_positions': spikes_positions_param,
            }
            random_parameters_images_dict[image_name] = random_parameters_dict
            image_dict[DATA] = self.add_artifact(
                image_dict.as_sitk(),
                spikes_positions_param,
                intensity_param,
            )
            # Add channels dimension
            image_dict[DATA] = image_dict[DATA][np.newaxis, ...]
            image_dict[DATA] = torch.from_numpy(image_dict[DATA])
        sample.add_transform(self, random_parameters_images_dict)
        return sample

    @staticmethod
    def get_params(
            num_spikes_range: Tuple[int, int],
            intensity_range: Tuple[float, float],
            ) -> Tuple:
        ns_min, ns_max = num_spikes_range
        num_spikes_param = torch.randint(ns_min, ns_max + 1, (1,)).item()
        intensity_param = torch.FloatTensor(1).uniform_(*intensity_range)
        spikes_positions = torch.rand(num_spikes_param, 3).numpy()
        return spikes_positions, intensity_param.item()

    def add_artifact(
            self,
            image: sitk.Image,
            spikes_positions: np.ndarray,
            intensity_factor: float,
            ):
        array = sitk.GetArrayViewFromImage(image).transpose()
        spectrum = self.fourier_transform(array)
        shape = np.array(spectrum.shape)
        mid_shape = shape // 2
        indices = np.floor(spikes_positions * shape).astype(int)
        for index in indices:
            diff = index - mid_shape
            i, j, k = mid_shape + diff
            spectrum[i, j, k] = spectrum.max() * intensity_factor
            # If we wanted to add a pure cosine, we should add spikes to both
            # sides of k-space. However, having only one is a better
            # representation og the actual cause of the artifact in real
            # scans.
            #i, j, k = mid_shape - diff
            #spectrum[i, j, k] = spectrum.max() * intensity_factor
        result = np.real(self.inv_fourier_transform(spectrum))
        return result.astype(np.float32)

max_stars_repo_name: sergivalverde/torchio

max_issues_repo_name: sergivalverde/torchio

max_forks_repo_name: sergivalverde/torchio


text:
import argparse
import networkx as nx
from typing import Optional

from charge.repository import Repository
from charge.chargers import Charger, CDPCharger
from charge.collectors import HistogramCollector

_charger = None  # type: Optional[Charger]

def parse_arguments() -> str:
    """Parses command line arguments.

    Returns:
        The location of the repository to use, as specified using
                the -r command line option.
    """
    parser = argparse.ArgumentParser()
    parser.add_argument('-r', '--repository', type=str,
                        help='Path to Repository (zip file) to use for'
                        ' reference charges.')
    args = parser.parse_args()
    return args.repository

def init(repo_location: Optional[str] = None) -> None:
    """Create the charger from a repository.

    Gets the repository location from the command line if none is
    given.

    Args:
        repo_location: The path to the repository to use.

    """
    global _charger
    if not repo_location:
        repo_location = parse_arguments()
    repo = Repository.read(repo_location)
    _charger = CDPCharger(repo,
                          rounding_digits=3)

def charge(graph: nx.Graph, total_charge: int) -> None:
    """Adds charges to a molecule.

    Updates the graph object.

    Args:
        graph: The molecule to add charges to.
        total_charge: The total charge to solve for.
    """
    global _charger
    _charger.charge(graph, total_charge)

ext: py

reason: import networkx


text:
import gzip
import os
from urllib.request import urlretrieve

import numpy as np
from tqdm import tqdm

BASE_URL = "https://cpar.s3.amazonaws.com/"

def download_from_s3(file_name):

  download_url = BASE_URL + file_name
  save_path = os.path.join('data', file_name)
  urlretrieve(download_url, save_path)

def load_data(path=None):
    """Loads the cpar-char dataset.
    # Returns
        Tuple of Numpy arrays: `(x_train, y_train), (x_test, y_test)`.
    """
    files = ['digit_train-labels-idx1-ubyte.gz', 'digit_train-images-idx3-ubyte.gz',
             'digit_test-labels-idx1-ubyte.gz', 'digit_test-images-idx3-ubyte.gz']

    paths = []
    if path is None:
      if os.path.isdir('data') is not True:
        os.mkdir('data')
      for fname in tqdm(files):
          if os.path.exists(os.path.join('data', fname)) is False:
            download_from_s3(fname)
          paths.append(fname)

    with gzip.open(os.path.join('data', paths[0]), 'rb') as lbpath:
        y_train = np.frombuffer(lbpath.read(), np.uint8, offset=8)

    with gzip.open(os.path.join('data', paths[1]), 'rb') as imgpath:
        x_train = np.frombuffer(imgpath.read(), np.uint8,
                                offset=16).reshape(len(y_train), 28, 28)

    with gzip.open(os.path.join('data', paths[2]), 'rb') as lbpath:
        y_test = np.frombuffer(lbpath.read(), np.uint8, offset=8)

    with gzip.open(os.path.join('data', paths[3]), 'rb') as imgpath:
        x_test = np.frombuffer(imgpath.read(), np.uint8,
                               offset=16).reshape(len(y_test), 28, 28)

    return (x_train, y_train), (x_test, y_test)

reason: import numpy


text:
import numpy as np
from numpy.testing import assert_array_equal
from nose.tools import assert_equals, assert_almost_equals
from ...knn.item_strategies import AllPossibleItemsStrategy, ItemsNeighborhoodStrategy
from ....models.classes import  MatrixPreferenceDataModel, \
     MatrixBooleanPrefDataModel
from ..classes import MatrixFactorBasedRecommender

movies = {'Marcel Caraciolo': {'Lady in the Water': 2.5, 'Snakes on a Plane': 3.5,
 'Just My Luck': 3.0, 'Superman Returns': 3.5, 'You, Me and Dupree': 2.5,
 'The Night Listener': 3.0},
'Luciana Nunes': {'Lady in the Water': 3.0, 'Snakes on a Plane': 3.5,
 'Just My Luck': 1.5, 'Superman Returns': 5.0, 'The Night Listener': 3.0,
 'You, Me and Dupree': 3.5},
'Leopoldo Pires': {'Lady in the Water': 2.5, 'Snakes on a Plane': 3.0,
 'Superman Returns': 3.5, 'The Night Listener': 4.0},
'Lorena Abreu': {'Snakes on a Plane': 3.5, 'Just My Luck': 3.0,
 'The Night Listener': 4.5, 'Superman Returns': 4.0,
 'You, Me and Dupree': 2.5},
'Steve Gates': {'Lady in the Water': 3.0, 'Snakes on a Plane': 4.0,
 'Just My Luck': 2.0, 'Superman Returns': 3.0, 'The Night Listener': 3.0,
 'You, Me and Dupree': 2.0},
'Sheldom': {'Lady in the Water': 3.0, 'Snakes on a Plane': 4.0,
 'The Night Listener': 3.0, 'Superman Returns': 5.0, 'You, Me and Dupree': 3.5},
'Penny Frewman': {'Snakes on a Plane': 4.5, 'You, Me and Dupree': 1.0, 'Superman Returns': 4.0},
'Maria Gabriela': {}}

matrix_model = MatrixPreferenceDataModel(movies)
boolean_matrix_model = MatrixBooleanPrefDataModel(movies)

def test_create_MatrixFactorBasedRecommender():
    items_strategy = AllPossibleItemsStrategy()
    recsys = MatrixFactorBasedRecommender(
        model=matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)
    assert_equals(recsys.items_selection_strategy, items_strategy)
    assert_equals(recsys.model, matrix_model)
    assert_equals(recsys.capper, True)
    assert_equals(recsys.learning_rate, 0.01)
    assert_equals(recsys.regularization, 0.02)
    assert_equals(recsys.init_mean, 0.1)
    assert_equals(recsys.n_interations, 30)
    assert_equals(recsys.init_stdev, 0.1)
    assert_equals(recsys.with_preference, False)
    assert_equals(recsys.user_factors.shape, (8, 2))
    assert_equals(recsys.item_factors.shape, (6, 2))
    assert_equals(recsys._global_bias, 3.2285714285714286)

def test_all_other_items_MatrixFactorBasedRecommender():
    items_strategy = AllPossibleItemsStrategy()
    recsys = MatrixFactorBasedRecommender(
        model=matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)

    assert_array_equal(np.array(['Lady in the Water']), recsys.all_other_items('Lorena Abreu'))
    assert_array_equal(np.array([], dtype='|S'), recsys.all_other_items('Marcel Caraciolo'))
    assert_array_equal(np.array(['Just My Luck', 'Lady in the Water', 'Snakes on a Plane',
       'Superman Returns', 'The Night Listener', 'You, Me and Dupree']), recsys.all_other_items('Maria Gabriela'))

    recsys = MatrixFactorBasedRecommender(
        model=matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)

    assert_array_equal(np.array(['Lady in the Water']), recsys.all_other_items('Lorena Abreu'))
    assert_array_equal(np.array([], dtype='|S'), recsys.all_other_items('Marcel Caraciolo'))
    assert_array_equal(np.array(['Just My Luck', 'Lady in the Water', 'Snakes on a Plane',
       'Superman Returns', 'The Night Listener', 'You, Me and Dupree']), recsys.all_other_items('Maria Gabriela'))

    recsys = MatrixFactorBasedRecommender(
        model=matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)

    assert_array_equal(np.array(['Lady in the Water']), recsys.all_other_items('Lorena Abreu'))
    assert_array_equal(np.array([], dtype='|S'), recsys.all_other_items('Marcel Caraciolo'))
    assert_array_equal(np.array(['Just My Luck', 'Lady in the Water', 'Snakes on a Plane',
       'Superman Returns', 'The Night Listener', 'You, Me and Dupree']), recsys.all_other_items('Maria Gabriela'))

    recsys = MatrixFactorBasedRecommender(
        model=boolean_matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)

    assert_array_equal(np.array(['Lady in the Water']), recsys.all_other_items('Lorena Abreu'))
    assert_array_equal(np.array([], dtype='|S'), recsys.all_other_items('Marcel Caraciolo'))
    assert_array_equal(np.array(['Just My Luck', 'Lady in the Water', 'Snakes on a Plane',
       'Superman Returns', 'The Night Listener', 'You, Me and Dupree']), recsys.all_other_items('Maria Gabriela'))

    recsys = MatrixFactorBasedRecommender(
        model=boolean_matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)

    assert_array_equal(np.array(['Lady in the Water']), recsys.all_other_items('Lorena Abreu'))
    assert_array_equal(np.array([], dtype='|S'), recsys.all_other_items('Marcel Caraciolo'))
    assert_array_equal(np.array(['Just My Luck', 'Lady in the Water', 'Snakes on a Plane',
       'Superman Returns', 'The Night Listener', 'You, Me and Dupree']), recsys.all_other_items('Maria Gabriela'))

    recsys = MatrixFactorBasedRecommender(
        model=boolean_matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)

    assert_array_equal(np.array(['Lady in the Water']), recsys.all_other_items('Lorena Abreu'))
    assert_array_equal(np.array([], dtype='|S'), recsys.all_other_items('Marcel Caraciolo'))
    assert_array_equal(np.array(['Just My Luck', 'Lady in the Water', 'Snakes on a Plane',
       'Superman Returns', 'The Night Listener', 'You, Me and Dupree']), recsys.all_other_items('Maria Gabriela'))

def test_estimate_preference_MatrixFactorBasedRecommender():
    items_strategy = ItemsNeighborhoodStrategy()
    recsys = MatrixFactorBasedRecommender(
        model=matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)
    assert_almost_equals(3.5, recsys.estimate_preference('Marcel Caraciolo', 'Superman Returns'))
    assert_almost_equals(3.206, recsys.estimate_preference('Leopoldo Pires', 'You, Me and Dupree'), 1)

    recsys = MatrixFactorBasedRecommender(
        model=matrix_model,
        items_selection_strategy=items_strategy,
        n_features=3)
    assert_almost_equals(3.5, recsys.estimate_preference('Marcel Caraciolo', 'Superman Returns'))
    assert_almost_equals(3.21,
         recsys.estimate_preference(user_id='Leopoldo Pires', item_id='You, Me and Dupree'), 1)

    #With capper = False
    recsys = MatrixFactorBasedRecommender(
        model=matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2, capper=False)
    assert_almost_equals(3.23, recsys.estimate_preference('Leopoldo Pires', 'You, Me and Dupree'), 1)

    #Boolean Matrix Model
    recsys = MatrixFactorBasedRecommender(
        model=boolean_matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)
    assert_almost_equals(1.0, recsys.estimate_preference('Marcel Caraciolo', 'Superman Returns'))
    assert_almost_equals(0.72, recsys.estimate_preference('Leopoldo Pires', 'You, Me and Dupree'), 1)
    assert_almost_equals(0.73,
         recsys.estimate_preference(user_id='Leopoldo Pires', item_id='You, Me and Dupree'), 1)

    #With capper = False
    recsys = MatrixFactorBasedRecommender(
        model=boolean_matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2, capper=False)
    assert_almost_equals(0.73, recsys.estimate_preference('Leopoldo Pires', 'You, Me and Dupree'), 1)
    #Non-Preferences
    assert_almost_equals(0.7258, recsys.estimate_preference('Maria Gabriela', 'You, Me and Dupree'), 1)

def test_recommend_MatrixFactorBasedRecommender():
    items_strategy = ItemsNeighborhoodStrategy()
    #Empty Recommendation
    recsys = MatrixFactorBasedRecommender(
        model=matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)
    assert_array_equal(np.array([]), recsys.recommend('Marcel Caraciolo'))

    #Semi Recommendation
    recsys = MatrixFactorBasedRecommender(
        model=matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)
    assert_array_equal(np.array(['You, Me and Dupree', 'Just My Luck']), \
        recsys.recommend('Leopoldo Pires'))

    #Semi Recommendation
    assert_array_equal(np.array(['You, Me and Dupree']), \
        recsys.recommend('Leopoldo Pires', 1))

    #Empty Recommendation
    recsys = MatrixFactorBasedRecommender(
        model=matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)
    assert_array_equal(np.array([]), recsys.recommend('Maria Gabriela'))

    #Test with params update
    recsys.recommend(user_id='Maria Gabriela', n_features=2)
    assert_array_equal(np.array([]), recsys.recommend('Maria Gabriela'))

    #with_preference
    #recsys = MatrixFactorBasedRecommender(
    #    model=matrix_model,
    #    items_selection_strategy=items_strategy,
    #    n_features=2, with_preference=True)
    #assert_array_equal(np.array([('Just My Luck', 3.20597319063), \
    #            ('You, Me and Dupree', 3.14717875510)]), \
    #            recsys.recommend('Leopoldo Pires'))

    #Empty Recommendation
    recsys = MatrixFactorBasedRecommender(
        model=boolean_matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)
    assert_array_equal(np.array([]), recsys.recommend('Marcel Caraciolo'))

    #Semi Recommendation
    recsys = MatrixFactorBasedRecommender(
        model=boolean_matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)
    assert_array_equal(np.array(['Just My Luck', 'You, Me and Dupree']), \
    recsys.recommend('Leopoldo Pires'))

    #Semi Recommendation
    recsys = MatrixFactorBasedRecommender(
        model=boolean_matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)
    assert_array_equal(np.array(['You, Me and Dupree']), \
        recsys.recommend('Leopoldo Pires', 1))

    #Empty Recommendation
    recsys = MatrixFactorBasedRecommender(
        model=boolean_matrix_model,
        items_selection_strategy=items_strategy,
        n_features=2)
    assert_array_equal(np.array([]), recsys.recommend('Maria Gabriela'))

    #Test with params update
    recsys.recommend(user_id='Maria Gabriela', n_features=2)
    assert_array_equal(np.array([]), recsys.recommend('Maria Gabriela'))

    #with_preference
    #recsys = MatrixFactorBasedRecommender(
    #    model=boolean_matrix_model,
    #    items_selection_strategy=items_strategy,
    #    n_features=2)
    #assert_array_equal(np.array([('Just My Luck', 3.20597), \
    #            ('You, Me and Dupree', 3.1471)]), \
    #            recsys.recommend('Leopoldo Pires'))

max_stars_repo_name: CrazyCaspian/crab

max_issues_repo_name: qqwjq/crab

max_forks_repo_name: qqwjq/crab

reason: import numpy,from numpy


text:
import socket

import matplotlib.pyplot as plt
import numpy as np

from tools import config

# Bind the TCP connection with only on device
s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
s.bind((config.HOST, config.PORT))
s.listen(1)

# Connect with client
conn, addr = s.accept()
print("Connected by ", addr)

# Saved temperature
temp = []

while True:
    # Receive data from client
    data = conn.recv(1024)

    # Check whether the clock can be used or not
    if "end" in data.decode():
        break
    else:
        print(data.decode())
        temp.append(float(data.decode()))

print("The clock is over")
conn.close()

plt.scatter(np.arange(len(temp)), temp)
plt.show()


text:
import pybel
import openbabel
import numpy as np
import logging

# Class for representing box
class Box:
    """
    This class represents the dimensions based on the molecule geometry
    """

    def __init__(self, xyzr, h=0.5, ext=4):
        self.logger = logging.getLogger(__name__)
        self.xmin = []
        self.xmax = []
        self.ymin = []
        self.ymax = []
        self.zmin = []
        self.zmax = []
        self.h = h
        self.ext = ext
        self.nx = []
        self.ny = []
        self.nz = []
        self.x = []
        self.y = []
        self.z = []
        self.xyzr = xyzr

        """
        This function sets the box dimensions based on the molecule geometry.
        INPUT:
            h: grid spacing in x,y, and z directions of the mesh on to
                         which the molecule geometry is mapped.
            xyzr: array of x,y,z, and radius values for each atom in the molecule
            ext:  distance of the box edge from the atom closest to the edge.

        OUTPUT:
            pybelmol: the new pybel molecule object
            total_chg: total charge on the new molecule
            mol3Dfile:  the .mol file of the adduct
        """
        # Determine the edge positions of the box along the x-axis
        [xleft, xright] = self.getBoxEdges(0)

        # Determine the edge positions of the box along the x-axis
        [yleft, yright] = self.getBoxEdges(1)

        # Determine the edge positions of the box along the x-axis
        [zleft, zright] = self.getBoxEdges(2)

        nx = np.int(np.round((xright - xleft) / self.h) + 1)
        ny = np.int(np.round((yright - yleft) / self.h) + 1)
        nz = np.int(np.round((zright - zleft) / self.h) + 1)

        self.logger.debug("nx = %s ny = %s nz = %s", nx, ny, nz)

        # Correct the values of edge distance based on nx, ny, and nz values.
        xright = xleft + (nx - 1) * self.h
        yright = yleft + (ny - 1) * self.h
        zright = zleft + (nz - 1) * self.h

        self.xleft = xleft
        self.xright = xright
        self.yleft = yleft
        self.yright = yright
        self.zleft = zleft
        self.zright = zright

        self.nx = nx
        self.ny = ny
        self.nz = nz

        for i in range(0, nx):
            xi = self.xleft + i * self.h
            self.x.append(xi)

        for i in range(0, ny):
            yi = self.yleft + i * self.h
            self.y.append(yi)

        for i in range(0, nz):
            zi = self.zleft + i * self.h
            self.z.append(zi)

    # Function to determine the box length dimension along an axis direction
    def getBoxEdges(self, index):
          # Determine the edge positions of the box along the axis
        natoms = len(self.xyzr)
        patom_min = self.xyzr[0][index] - self.xyzr[0][3] - self.ext
        patom_max = self.xyzr[0][index] + self.xyzr[0][3] + self.ext

        for i in range(1, natoms):  # i:1, 2, ...,natoms - 1
            tmp = self.xyzr[i][index] + self.xyzr[i][3] + self.ext  # tmp = p + rad + ext
            if patom_max < tmp:
                patom_max = tmp
            tmp = self.xyzr[i][index] - self.xyzr[i][3] - self.ext  # tmp = p - rad - ext
            if patom_min > tmp:
                patom_min = tmp

        left_ex = patom_min / self.h
        right_ex = patom_max / self.h

        pleft = np.floor(left_ex) * self.h - self.ext
        pright = np.ceil(right_ex) * self.h + self.ext
        return pleft, pright

    # create atom accessibility map
    def createAtomAccessibiltyMap(self, prob_rad):
        """
        This function creates the molecule surface accessibility map of the new
        atom to be added in the molecule.
        INPUT:
            prob_rad: radius of the atom to be added

        OUTPUT:
            pybelmol: the new pybel molecule object
            total_chg: total charge on the new molecule
            mol3Dfile:  the .mol file of the adduct
        """

        # initialize the values of the accessibilty map to zero
        acc_map = np.zeros((self.nx, self.ny, self.nz))

        natoms = len(self.xyzr)

        # Create the  accessibility map by
        # rolling a spherical particle of radius (= probe radius) on the surface of each atom
        # to determine the domain accessible/in-accessible to the particle.
        for i in range(0, natoms):
            atom_radius = self.xyzr[i][3]
            atom_x = self.xyzr[i][0]
            atom_y = self.xyzr[i][1]
            atom_z = self.xyzr[i][2]

            r_atom = atom_radius + prob_rad

            x1 = np.floor((atom_x - self.xleft - r_atom) / self.h)
            x2 = np.ceil((atom_x - self.xleft + r_atom) / self.h)
            ix1 = np.int(x1)
            ix2 = np.int(x2)

            y1 = np.floor((atom_y - self.yleft - r_atom) / self.h)
            y2 = np.ceil((atom_y - self.yleft + r_atom) / self.h)
            iy1 = np.int(y1)
            iy2 = np.int(y2)

            z1 = np.floor((atom_z - self.zleft - r_atom) / self.h)
            z2 = np.ceil((atom_z - self.zleft + r_atom) / self.h)
            iz1 = np.int(z1)
            iz2 = np.int(z2)

            for ix in range(ix1, ix2 + 1):
                for iy in range(iy1, iy2 + 1):
                    for iz in range(iz1, iz2 + 1):

                        xdist = self.x[ix] - atom_x
                        ydist = self.y[iy] - atom_y
                        zdist = self.z[iz] - atom_z
                        rdist = np.sqrt(xdist**2 + ydist**2 + zdist**2)

                        if rdist < r_atom:
                            acc_map[ix, iy, iz] = 1.0
        return acc_map

def addAtomToMol(mol, atom, idx, covalent=True):
    """
    This function adds an atom (atom) to a molecule (mol) at the site
    near an atom (idx). It uses the covalent radius of the atoms
    to determine the accessible site for placing the the new atom.

    INPUT:
        mol - pybel molecule object
        atom - element symbol (string) of the atom to be added
        idx - atom number of the atom near which the new atom is
              placed
        covalent - True for covalent; otherwise non-covalent (van der Waal)
    OUTPUT:
        mol - pybel molecule object containing the new atom
    """

    logger = logging.getLogger(__name__)
    etab = openbabel.OBElementTable()     # element object

    atoms_radius = []
    xyzr = []

    for x in mol.atoms:
        atomic_num = x.atomicnum  # get atomic number of each atom
        if covalent is True:
            # use covalent radius
            atom_rad = etab.GetCovalentRad(atomic_num)

        else:
            # use van der Waals radius (but for now use covalent)
            atom_rad = etab.GetVdwRad(atomic_num)
        atoms_radius.append(atom_rad)

        a = x.coords + (atom_rad,)  # create x,y,z,radius array
        xyzr.append(a)

    # get the atomic number and covalent radius of the new atom to be added
    newatom_atomic_num = etab.GetAtomicNum(atom)

    if covalent is True:
        newatom_rad = etab.GetVdwRad(newatom_atomic_num)
    else:
        newatom_rad = etab.GetVdwRad(1)
    newatom_xyzr = [0, 0, 0, newatom_rad]

    # set box dimensions based on molecule size
    box = Box(xyzr, 0.5, 10)

    # create accessibilty map
    probe_radius = newatom_rad
    acc_map = box.createAtomAccessibiltyMap(probe_radius)
    logger.info("Finished creating atom accessibility map.")
    xc = xyzr[idx][0]
    yc = xyzr[idx][1]
    zc = xyzr[idx][2]
    rc = xyzr[idx][3]

    rdist = rc + probe_radius
    logger.info("rdist = %s", rdist)

    # randomly select a position for the new atom near the specified atom of the
    # molecule
    found = 0
    delr = 0.01
    ntrials = 5000
    ndist = 50
    count = 0

    logger.info("Selecting a position for the new atom near the atom located at: (%s, %s, %s)", xc, yc, zc)

    while (found == 0 and count < ndist):
        logger.debug("count #%s", count)
        rdist = rdist + delr
        ntrials = 5000
        trial = 0

        x1 = xc - rdist
        x2 = xc + rdist
        if x1 < box.xleft:
            x1 = box.xleft

        if x2 > box.xright:
            x2 = box.xright

        y1 = yc - rdist
        y2 = yc + rdist
        if y1 < box.yleft:
            y1 = box.yleft

        if y2 > box.yright:
            y2 = box.yright

        z1 = zc - rdist
        z2 = zc + rdist
        if z1 < box.zleft:
            z1 = box.zleft

        if z2 > box.zright:
            z2 = box.zright

        logger.debug("x1, y1, z1: %s %s %s", x1, y1, z1)
        logger.debug("x2, y2, z2: %s %s %s", x2, y2, z2)
        while trial < ntrials and found == 0:
            logger.debug("trial #: %s", trial)
            rand_x = x1 + np.random.ranf() * (x2 - x1)
            rand_y = y1 + np.random.ranf() * (y2 - y1)
            rand_z = z1 + np.random.ranf() * (z2 - z1)

            ix = np.int((rand_x - box.xleft) / box.h)
            iy = np.int((rand_y - box.yleft) / box.h)
            iz = np.int((rand_z - box.zleft) / box.h)

            if acc_map[ix, iy, iz] == 0:
                newatom_xyzr[0] = rand_x
                newatom_xyzr[1] = rand_y
                newatom_xyzr[2] = rand_z
                found = 1
                logger.info("Found a site: (%s, %s, %s)", rand_x, rand_y, rand_z)
                rd = np.sqrt((rand_x - xc) * (rand_x - xc) + (rand_y - yc) * (rand_y - yc) +
                             (rand_z - zc) * (rand_z - zc))
                logger.info("inter-atom distance (A): %s", rd)

            else:
                trial = trial + 1
                found = 0

        count = count + 1

    if found == 0:
        logger.critical("Could not find a site to place adduct")
        raise Exception("Could not find a site to place adduct")
    else:
        newatom = mol.OBMol.NewAtom()
        newatom.SetAtomicNum(newatom_atomic_num)
        newatom.SetVector(rand_x, rand_y, rand_z)

        xyzr.append(newatom_xyzr)

        if covalent is True:
            new_atomnum = len(mol.atoms)

            mol.OBMol.AddBond(idx + 1, new_atomnum, 1)

            # good for single bonds. For double and triple bonds, we have to
            # the appropriate covalent radius. For our purposes, we only need
            # single bond (for protonation).

    return mol

def removeAtomFromMol(mol, idx):
    """
    This function removes an atom from the pybel moleule object.
    """

    # mol should contain the atom coordinates
    atom = mol.OBMol.GetAtomById(idx)

    # delete the atom
    mol.OBMol.DeleteAtom(atom)

    return mol

def nearestHydrogen(mol, idx):
    logger = logging.getLogger(__name__)

    iatom = mol.atoms[idx]
    logger.debug('Starting atom: %s, type %s', idx, iatom.atomicnum)

    # get the neighboring atoms of the selected atom
    nbatoms = openbabel.OBAtomAtomIter(iatom.OBAtom)
    for nb in nbatoms:
        logger.debug('Connected to: %s, type %s', nb.GetId(), nb.GetAtomicNum())
        # if the neighboring atom is hydrogen
        if nb.GetAtomicNum() == 1:
            return nb.GetId()

    raise Exception('Hydrogen not found.')

max_stars_repo_name: pnnl/isicle

max_issues_repo_name: pnnl/isicle

max_forks_repo_name: pnnl/isicle


text:
import numpy as np
import scipy as sp
import scipy.stats
import numpy

def confidenceinterval(array):
    confidence = 0.95
    a = 1.0 * np.array(array)
    m = np.mean(a)
    fc = scipy.stats.sem(a)
    h = fc * sp.stats.t._ppf((1 + confidence) / 2., n - 1) / ((n - 1) ** 0.5)
    return m - h, m + h

n = 10
a = [71,59,58]
array = numpy.asarray(a)
print(confidenceinterval(array))

reason: import numpy,import scipy


text:
import numpy as np
from numpy.testing import assert_array_equal
from nose import with_setup
try:
    from nose.tools import assert_is
except ImportError:
    from landlab.testing.tools import assert_is
from nose.tools import assert_equal

from landlab.grid.raster_steepest_descent import (
    _calc_steepest_descent_across_adjacent_cells)

def setup_unit_grid():
    """Set up a test grid with unit spacing."""
    from landlab import RasterModelGrid
    globals()['rmg'] = RasterModelGrid(4, 5)
    globals()['values_at_nodes'] = np.arange(20, dtype=float)

def setup_non_unit_grid():
    """Set up a test grid with non-unit spacing."""
    from landlab import RasterModelGrid
    globals()['rmg'] = RasterModelGrid(4, 5, 2)
    globals()['values_at_nodes'] = np.arange(20, dtype=float)

def setup_3x3_grid():
    """Set up a grid of 3 rows and 3 columns."""
    from landlab import RasterModelGrid
    globals().update({
        'rmg_3x3': RasterModelGrid(3, 3),
        'values_at_nodes': np.flipud(np.array([6, 7, 8,
                                               3, 4, 5,
                                               0, 1, 2], dtype=float))
    })

@with_setup(setup_3x3_grid)
def test_scalar_arg():
    """Test scalar arg for nodes."""
    grad = _calc_steepest_descent_across_adjacent_cells(
        rmg_3x3, values_at_nodes, 0)
    assert_equal(grad, -3.)

    grad = _calc_steepest_descent_across_adjacent_cells(
        rmg_3x3, values_at_nodes, 0, method='d8')
    assert_equal(grad, -3.)

    values_at_nodes[2] = -10
    grad = _calc_steepest_descent_across_adjacent_cells(
        rmg_3x3, values_at_nodes, 0, method='d8')
    assert_equal(grad, - (4 + 10) / np.sqrt(2.))

@with_setup(setup_unit_grid)
def test_iterable():
    """Test iterable arg for nodes."""
    grad = _calc_steepest_descent_across_adjacent_cells(
        rmg, values_at_nodes, [0, 4])
    assert_array_equal(grad, [-5., -5.])

@with_setup(setup_unit_grid)
def test_scalar_arg_with_links():
    values = np.array([0, 1,  3, 6, 10,
                       0, 1,  3, 6, 10,
                       0, 1,  3, 5, 10,
                       0, 1, -3, 6, 10, ], dtype=float)
    (grad, node) = _calc_steepest_descent_across_adjacent_cells(
        rmg, values, (0, 4), return_node=True)
    assert_array_equal(grad, [-1, -6])
    assert_array_equal(node, [5, 17])

    values_at_nodes[2] = -10
    (grad, node) = _calc_steepest_descent_across_adjacent_cells(
        rmg, values_at_nodes, 0, method='d8', return_node=True)
    assert_equal(grad, - (6 + 10) / np.sqrt(2.))
    assert_equal(node, 2)

@with_setup(setup_unit_grid)
def test_node_id_in_direction_of_max():
    values = np.array([-1, 1,  3, 6, 10,
                       0, 1,  3, 6, 10,
                       0, 1,  3, 5, 10,
                       0, 1, -3, 6, 10, ], dtype=float)
    (_, node_ids) = _calc_steepest_descent_across_adjacent_cells(
        rmg, values, (0, 4), return_node=True)
    assert_array_equal(node_ids, [5, 17])

    (grads,
     node_ids) = _calc_steepest_descent_across_adjacent_cells(
        rmg, values, (0, 4), method='d8', return_node=True)
    assert_array_equal(node_ids, [0, 17])

@with_setup(setup_3x3_grid)
def test_node_in_direction_of_max():
    for node_id in [0, 1, 2, 3, 5, 6, 7, 8]:
        values = np.zeros(9)
        values[node_id] = -1
        (_, node) = _calc_steepest_descent_across_adjacent_cells(
            rmg_3x3, values, 0, return_node=True, method='d8')
        assert_array_equal(node, node_id)

@with_setup(setup_3x3_grid)
def test_node_in_direction_of_max_with_ties():
    values = np.zeros(9)
    (_, node) = _calc_steepest_descent_across_adjacent_cells(
        rmg_3x3, values, 0, return_node=True, method='d8')
    assert_array_equal(node, 5)

    for (node_id, expected) in zip([5, 7, 3, 1, 8, 6, 0],
                                   [7, 3, 1, 8, 6, 0, 2]):
        values[node_id] = 1
        (_, node) = _calc_steepest_descent_across_adjacent_cells(
            rmg_3x3, values, 0, return_node=True, method='d8')
        assert_array_equal(node, expected)

max_stars_repo_name: SiccarPoint/landlab

max_issues_repo_name: csherwood-usgs/landlab

max_forks_repo_name: csherwood-usgs/landlab


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import urllib.request
import sys
import time
import urllib
import os.path

from struct import pack

from sklearn import preprocessing

import numpy as np
import pandas as pd
import fact.io
import progressbar

N = 200000
factLevel2URL = "https://factdata.app.tu-dortmund.de/dl2/FACT-Tools/v1.1.1/"
gammaName = "gamma_simulations_facttools_dl2.hdf5"
protonName = "proton_simulations_facttools_dl2.hdf5"
outFolder = "./"

def reporthook(count, block_size, total_size):
    global start_time
    if count == 0:
        start_time = time.time()
        return
    duration = time.time() - start_time
    progress_size = int(count * block_size)
    speed = int(progress_size / (1024 * duration))
    percent = int(count * block_size * 100 / total_size)
    sys.stdout.write("\r...%d%%, %d MB, %d KB/s, %d seconds passed" %
                    (percent, progress_size / (1024 * 1024), speed, duration))
    sys.stdout.flush()

if not os.path.isfile(outFolder + "/" + gammaName):
    print("Fact gamma data not found, downloading it from", factLevel2URL + "\n")
    urllib.request.urlretrieve(factLevel2URL + gammaName, outFolder + "/" + gammaName, reporthook)

if not os.path.isfile(outFolder + "/" + protonName):
    print("Fact proton data not found, downloading it from", factLevel2URL + "\n")
    urllib.request.urlretrieve(factLevel2URL + protonName, outFolder + "/" + protonName, reporthook)

print("Loading ", outFolder + "/" + gammaName, " with FACT tools")
#gammaDF = fact.io.read_data(outFolder + "/" + gammaName, key='events', first=_first+shift, last=_last+shift)
gammaDF = fact.io.read_data(outFolder + "/" + gammaName, key='events')
gammaDF["label"] = False

#print("{} item in gammaDF".format(len(gammaDF)))

print("Loading ", outFolder + "/" + protonName, " with FACT tools")
protonDF = fact.io.read_data(outFolder + "/" + protonName, key='events')
protonDF["label"] = True

#print("{} items in protonDF".format(len(protonDF)))

#asdf
header = [
    "concentration_cog",
    "concentration_core",
    "concentration_one_pixel",
    "concentration_two_pixel",
    "leakage1",
    "leakage2",
    "size",
    "width",
    "length",
    "skewness_long",
    "skewness_trans",
    "kurtosis_long",
    "kurtosis_trans",
    "num_islands",
    "num_pixel_in_shower",
    "photoncharge_shower_mean",
    "photoncharge_shower_variance",
    "area",
    "log_size",
    "log_length",
    "size_area",
    "area_size_cut_var"
]

#print(protonDF.columns.tolist())
#gammaDF = gammaDF[gammaDF.columns.tolist()]
qualityCuts = "num_pixel_in_shower >= 10 & num_islands < 8 & length < 70 & width < 35 & leakage1 < 0.6 & leakage2 < 0.85"

print("Applying sanity filtering")
#gammaDF = gammaDF.replace([np.inf, -np.inf], np.nan).dropna(axis=0, how='any')
#protonDF = protonDF.replace([np.inf, -np.inf], np.nan).dropna(axis=0, how='any')

print("Applying quality cuts")
gammaDF = gammaDF.query(qualityCuts)
protonDF = protonDF.query(qualityCuts)

print("Num gammas after filtering+cuts: ", len(gammaDF))
print("Num protons after filtering+cuts: ", len(protonDF))
print("Generate some additional features taken from Jens' YML script")

# There is a "chained assignment" warning for the following computations
# However, I dont see a reason why this should be a problem. Thus we disable the warning
pd.options.mode.chained_assignment = None

gammaDF["area"] = gammaDF["width"]*gammaDF["length"]*np.pi
gammaDF["log_size"] = np.log(gammaDF["size"])
gammaDF["log_length"] = np.log(gammaDF["length"])
gammaDF["size_area"] = gammaDF["size"]/(gammaDF["width"]*gammaDF["length"]*np.pi)
gammaDF["area_size_cut_var"] = (gammaDF["width"]*gammaDF["length"]*np.pi) /(np.log(gammaDF["size"])**2)

protonDF["area"] = protonDF["width"]*protonDF["length"]*np.pi
protonDF["log_size"] = np.log(protonDF["size"])
protonDF["log_length"] = np.log(protonDF["length"])
protonDF["size_area"] = protonDF["size"]/(protonDF["width"]*protonDF["length"]*np.pi)
protonDF["area_size_cut_var"] = (protonDF["width"]*protonDF["length"]*np.pi) /(np.log(protonDF["size"])**2)

print("Selecting rows")
gammaDF = gammaDF[header + ["label"]]
protonDF = protonDF[header + ["label"]]

print("Shuffle data")
gammaDF = gammaDF.sample(frac=1).reset_index(drop=True)
protonDF = protonDF.sample(frac=1).reset_index(drop=True)

print("Sampling data")
#dfTrain = pd.concat([gammaDF.sample(n=int(NTrain/2)), protonDF.sample(n=int(NTrain/2))])
df = pd.concat([gammaDF.loc[0:int(N/2)-1,:], protonDF.loc[0:int(N/2)-1,:]])
df = df.sample(frac=1).reset_index(drop=True)
df.to_csv("data.csv", sep=",", index=False, header=True)


text:
import tensorflow as tf
# import matplotlib.pyplot as plt
#import numpy as np
from numpy.random import RandomState

def get_weight(shape, lamb):
	var=tf.Variable(tf.random_normal(shape),dtype = tf.float32)
	tf.add_to_collection('losses',tf.contrib.layers.l2_regularizer(lamb)(var))
	return var

batch_size=8
w1 = tf.Variable(tf.random_normal([2,3],stddev=1,seed=1))
#w2 = tf.Variable(tf.random_normal([3,1],stddev=1,seed=1))
x  = tf.placeholder(tf.float32,shape=(None,2),name='x-input')
y_ = tf.placeholder(tf.float32,shape=(None,1),name='y-input')
layer_dimension=[2,10,10,10,1]
n_layers=len(layer_dimension)

cur_layer=x
in_dimension=layer_dimension[0]

for i in range(1,n_layers):
	out_dimension = layer_dimension[i]
	weight = get_weight([in_dimension,out_dimension],0.001)
	bias   = tf.Variable(tf.constant(0.1,shape=[out_dimension]))
	cur_layer=tf.nn.relu(tf.matmul(cur_layer,weight)+bias)
	in_dimension = layer_dimension[i]

# #a = tf.matmul(x,w1)
# y = tf.matmul(x,w1)

# corss_entropy = -tf.reduce_mean(y_ * tf.log(tf.clip_by_value(y,1e-10,1.0)))
mse_loss = tf.reduce_mean(tf.square(y_ - cur_layer))
tf.add_to_collection('losses',mse_loss)
loss = tf.add_n(tf.get_collection('losses'))

train_step = tf.train.AdamOptimizer(learning_rate=0.001).minimize(loss)

rdm = RandomState(1)
dataset_size=128
X=rdm.rand(dataset_size,2)
Y=[[x1+x2+rdm.rand()/10.0-0.05] for (x1,x2) in X ]

with tf.Session() as sess:
	init_op=tf.global_variables_initializer()
	sess.run(init_op)

	print(sess.run(w1))

	STEPS=15000
	for i in range(STEPS):
		start=(i*batch_size) % dataset_size
		end= min(start+batch_size,dataset_size)

		sess.run(train_step,feed_dict={x:X[start:end],y_:Y[start:end]})

		if i % 1000 == 0:
			total_cross_entropy = sess.run(loss,feed_dict={x:X,y_:Y})
			print("%d steps,corss_entropy on all data is %g" % (i,total_cross_entropy))

	print(sess.run(w1))

	# print(sess.run(y,feed_dict={x:[[0.7,0.9],[0.1,0.4],[0.5,0.8]]}))
# print(a.graph is tf.get_default_graph())

# g1 = tf.Graph()
# with g1.as_default():
# 	v = tf.get_variable("a",shape=[1],initializer=tf.zeros_initializer())

# g2 = tf.Graph()
# with g2.as_default():
# 	v = tf.get_variable("a",shape=[1],initializer=tf.ones_initializer())

# with tf.Session(graph=g1) as sess:
# 	tf.global_variables_initializer().run()
# 	with tf.variable_scope("",reuse=True):
# 		print(sess.run(tf.get_variable('a')))

# with tf.Session(graph=g2) as sess:
# 	tf.global_variables_initializer().run()
# 	with tf.variable_scope("",reuse=True):
# 		print(sess.run(tf.get_variable('a')))

# a1 = tf.get_variable(name='a1', shape=[2,3], initializer=tf.random_normal_initializer(mean=0, stddev=1))
# a2 = tf.get_variable(name='a2', shape=[1], initializer=tf.constant_initializer(1))
# a3 = tf.get_variable(name='a3', shape=[2,3], initializer=tf.ones_initializer())

# # zbh: not done
# with tf.Session() as sess:
# 	tf.global_variables_initializer().run()
# 	sess.run()
#     #sess.run(tf.initialize_all_variables())
#     print(sess.run(a1))
#     print(sess.run(a2))
#     print(sess.run(a3))

max_stars_repo_name: zbhoscar/tfelement

max_issues_repo_name: zbhoscar/tfelement

max_forks_repo_name: zbhoscar/tfelement


text:
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
import cv2
import numpy as np

class Preprocessor(object):
    """Defines an abstract observation preprocessor function.

    Attributes:
        shape (obj): Shape of the preprocessed output.
    """

    def __init__(self, obs_space, options):
        self._obs_space = obs_space
        self._options = options
        self._init()

    def _init(self):
        pass

    def transform(self, observation):
        """Returns the preprocessed observation."""
        raise NotImplementedError

class AtariPixelPreprocessor(Preprocessor):
    def _init(self):
        self._grayscale = self._options.get("grayscale", False)
        self._zero_mean = self._options.get("zero_mean", True)
        self._dim = self._options.get("dim", 80)
        self._channel_major = self._options.get("channel_major", False)
        if self._grayscale:
            self.shape = (self._dim, self._dim, 1)
        else:
            self.shape = (self._dim, self._dim, 3)

        # channel_major requires (# in-channels, row dim, col dim)
        if self._channel_major:
            self.shape = self.shape[-1:] + self.shape[:-1]

    # TODO(ekl) why does this need to return an extra size-1 dim (the [None])
    def transform(self, observation):
        """Downsamples images from (210, 160, 3) by the configured factor."""
        scaled = observation[25:-25, :, :]
        if self._dim < 80:
            scaled = cv2.resize(scaled, (80, 80))
        # OpenAI: Resize by half, then down to 42x42 (essentially mipmapping).
        # If we resize directly we lose pixels that, when mapped to 42x42,
        # aren't close enough to the pixel boundary.
        scaled = cv2.resize(scaled, (self._dim, self._dim))
        if self._grayscale:
            scaled = scaled.mean(2)
            scaled = scaled.astype(np.float32)
            # Rescale needed for maintaining 1 channel
            scaled = np.reshape(scaled, [self._dim, self._dim, 1])
        if self._zero_mean:
            scaled = (scaled - 128) / 128
        else:
            scaled *= 1.0 / 255.0
        if self._channel_major:
            scaled = np.reshape(scaled, self.shape)
        return scaled

# TODO(rliaw): Also should include the deepmind preprocessor
class AtariRamPreprocessor(Preprocessor):
    def _init(self):
        self.shape = (128,)

    def transform(self, observation):
        return (observation - 128) / 128

class OneHotPreprocessor(Preprocessor):
    def _init(self):
        assert self._obs_space.shape == ()
        self.shape = (self._obs_space.n,)

    def transform(self, observation):
        arr = np.zeros(self._obs_space.n)
        arr[observation] = 1
        return arr

class NoPreprocessor(Preprocessor):
    def _init(self):
        self.shape = self._obs_space.shape

    def transform(self, observation):
        return observation

max_stars_repo_name: rickyHong/Ray-repl


text:
'''

Classifier        :   Decision Tree
Dataset           :   Iris Dataset
Features          :   Sepal Width, Sepal Length, Petal Width, Petal Length
Labels            :   Versicolor, Setosa, Virginica
Training Dataset  :   147 Entries
Testing Dataset   :   3 Entries
Author            :   Prasad Dangare
Function Name     :   main

'''

from sklearn.datasets import load_iris
import numpy as np
from sklearn import tree

def main():

    dataset = load_iris()

    print("Features of dataset")
    print(dataset.feature_names)
    print("Target names of dataset")
    print(dataset.target_names)

    # Indices of removed elements
    index = [1, 51, 101]

    # Training data with removed elements
    train_target = np.delete(dataset.target, index)
    train_feature = np.delete(dataset.data, index, axis = 0)

    # Testing data for testing on trainning data
    test_target = dataset.target[index]
    test_feature = dataset.data[index]

    # form dcision tree classifier
    obj = tree.DecisionTreeClassifier()

    # apply training data to form tree
    obj.fit(train_feature, train_target)

    # apply testing on trained data
    result = obj.predict(test_feature)

    print("result prediction by ML ", result)
    print("result expected ", test_target)

if __name__ == "__main__":
    main()


text:
from chord_recognition import MIDIChord
import numpy as np
import pretty_midi
import copy

from vocab import RemiVocab, OctupleVocab

from constants import (
  # vocab keys
  TIME_SIGNATURE_KEY,
  BAR_KEY,
  POSITION_KEY,
  INSTRUMENT_KEY,
  PITCH_KEY,
  VELOCITY_KEY,
  DURATION_KEY,
  TEMPO_KEY,
  CHORD_KEY,
  NOTE_DENSITY_KEY,
  MEAN_PITCH_KEY,
  MEAN_VELOCITY_KEY,
  MEAN_DURATION_KEY,
  # discretization parameters
  DEFAULT_POS_PER_QUARTER,
  DEFAULT_VELOCITY_BINS,
  DEFAULT_DURATION_BINS,
  DEFAULT_TEMPO_BINS,
  DEFAULT_NOTE_DENSITY_BINS,
  DEFAULT_MEAN_VELOCITY_BINS,
  DEFAULT_MEAN_PITCH_BINS,
  DEFAULT_MEAN_DURATION_BINS,
  DEFAULT_RESOLUTION
)

# define "Item" for general storage
class Item(object):
  def __init__(self, name, start, end, velocity=None, pitch=None, instrument=None):
    self.name = name
    self.start = start
    self.end = end
    self.velocity = velocity
    self.pitch = pitch
    self.instrument = instrument

  def __repr__(self):
    return 'Item(name={}, start={}, end={}, velocity={}, pitch={}, instrument={})'.format(
      self.name, self.start, self.end, self.velocity, self.pitch, self.instrument)

# define "Event" for event storage
class Event(object):
  def __init__(self, name, time, value, text):
    self.name = name
    self.time = time
    self.value = value
    self.text = text

  def __repr__(self):
    return 'Event(name={}, time={}, value={}, text={})'.format(
      self.name, self.time, self.value, self.text)

class InputRepresentation():
  def version():
    return 'v4'

  def __init__(self, file, do_extract_chords=True, strict=False):
    if isinstance(file, pretty_midi.PrettyMIDI):
      self.pm = file
    else:
      self.pm = pretty_midi.PrettyMIDI(file)

    if strict and len(self.pm.time_signature_changes) == 0:
      raise ValueError("Invalid MIDI file: No time signature defined")

    self.resolution = self.pm.resolution

    self.note_items = None
    self.tempo_items = None
    self.chords = None
    self.groups = None

    self._read_items()
    self._quantize_items()
    if do_extract_chords:
      self.extract_chords()
    self._group_items()

    if strict and len(self.note_items) == 0:
      raise ValueError("Invalid MIDI file: No notes found, empty file.")

  # read notes and tempo changes from midi (assume there is only one track)
  def _read_items(self):
    # note
    self.note_items = []
    for instrument in self.pm.instruments:
      pedal_events = [event for event in instrument.control_changes if event.number == 64]
      pedal_pressed = False
      start = None
      pedals = []
      for e in pedal_events:
        if e.value >= 64 and not pedal_pressed:
          pedal_pressed = True
          start = e.time
        elif e.value < 64 and pedal_pressed:
          pedal_pressed = False
          pedals.append(Item(name='Pedal', start=start, end=e.time))
          start = e.time

      notes = instrument.notes
      notes.sort(key=lambda x: (x.start, x.pitch))

      if instrument.is_drum:
        instrument_name = 'drum'
      else:
        instrument_name = instrument.program

      pedal_idx = 0
      for note in notes:
        pedal_candidates = [(i + pedal_idx, pedal) for i, pedal in enumerate(pedals[pedal_idx:]) if note.end >= pedal.start and note.start < pedal.end]
        if len(pedal_candidates) > 0:
          pedal_idx = pedal_candidates[0][0]
          pedal = pedal_candidates[-1][1]
        else:
          pedal = Item(name='Pedal', start=0, end=0)

        self.note_items.append(Item(
          name='Note',
          start=self.pm.time_to_tick(note.start),
          end=self.pm.time_to_tick(max(note.end, pedal.end)),
          velocity=note.velocity,
          pitch=note.pitch,
          instrument=instrument_name))
    self.note_items.sort(key=lambda x: (x.start, x.pitch))
    # tempo
    self.tempo_items = []
    times, tempi = self.pm.get_tempo_changes()
    for time, tempo in zip(times, tempi):
      self.tempo_items.append(Item(
        name='Tempo',
        start=time,
        end=None,
        velocity=None,
        pitch=int(tempo)))
    self.tempo_items.sort(key=lambda x: x.start)
    # expand to all beat
    max_tick = self.pm.time_to_tick(self.pm.get_end_time())
    existing_ticks = {item.start: item.pitch for item in self.tempo_items}
    wanted_ticks = np.arange(0, max_tick+1, DEFAULT_RESOLUTION)
    output = []
    for tick in wanted_ticks:
      if tick in existing_ticks:
        output.append(Item(
          name='Tempo',
          start=self.pm.time_to_tick(tick),
          end=None,
          velocity=None,
          pitch=existing_ticks[tick]))
      else:
        output.append(Item(
          name='Tempo',
          start=self.pm.time_to_tick(tick),
          end=None,
          velocity=None,
          pitch=output[-1].pitch))
    self.tempo_items = output

  # quantize items
  def _quantize_items(self):
    ticks = self.resolution / DEFAULT_POS_PER_QUARTER
    # grid
    end_tick = self.pm.time_to_tick(self.pm.get_end_time())
    grids = np.arange(0, end_tick, ticks)
    # process
    for item in self.note_items:
      index = np.searchsorted(grids, item.start, side='right')
      if index > 0:
        index -= 1
      shift = round(grids[index]) - item.start
      item.start += shift
      item.end += shift

  def get_end_tick(self):
    return self.pm.time_to_tick(self.pm.get_end_time())

  # extract chord
  def extract_chords(self):
    method = MIDIChord(self.pm)
    chords = method.extract()
    output = []
    for chord in chords:
      output.append(Item(
        name='Chord',
        start=self.pm.time_to_tick(chord[0]),
        end=self.pm.time_to_tick(chord[1]),
        velocity=None,
        pitch=chord[2].split('/')[0]))
    if len(output) == 0 or output[0].start > 0:
      if len(output) == 0:
        end = self.pm.time_to_tick(self.pm.get_end_time())
      else:
        end = output[0].start
      output.append(Item(
        name='Chord',
        start=0,
        end=end,
        velocity=None,
        pitch='N:N'
      ))
    self.chords = output
    return self.chords

  # group items
  def _group_items(self):
    if self.chords:
      items = self.chords + self.tempo_items + self.note_items
    else:
      items = self.tempo_items + self.note_items

    def _get_key(item):
      type_priority = {
        'Chord': 0,
        'Tempo': 1,
        'Note': 2
      }
      return (
        item.start, # order by time
        type_priority[item.name], # chord events first, then tempo events, then note events
        -1 if item.instrument == 'drum' else item.instrument, # order by instrument
        item.pitch # order by note pitch
      )

    items.sort(key=_get_key)
    downbeats = self.pm.get_downbeats()
    downbeats = np.concatenate([downbeats, [self.pm.get_end_time()]])
    self.groups = []
    for db1, db2 in zip(downbeats[:-1], downbeats[1:]):
      db1, db2 = self.pm.time_to_tick(db1), self.pm.time_to_tick(db2)
      insiders = []
      for item in items:
        if (item.start >= db1) and (item.start < db2):
          insiders.append(item)
      overall = [db1] + insiders + [db2]
      self.groups.append(overall)

    # Trim empty groups from the beginning and end
    for idx in [0, -1]:
      while len(self.groups) > 0:
        group = self.groups[idx]
        notes = [item for item in group[1:-1] if item.name == 'Note']
        if len(notes) == 0:
          self.groups.pop(idx)
        else:
          break

    return self.groups

  def _get_time_signature(self, start):
    # This method assumes that time signature changes don't happen within a bar
    # which is a convention that commonly holds
    time_sig = None
    for curr_sig, next_sig in zip(self.pm.time_signature_changes[:-1], self.pm.time_signature_changes[1:]):
      if self.pm.time_to_tick(curr_sig.time) <= start and self.pm.time_to_tick(next_sig.time) > start:
        time_sig = curr_sig
        break
    if time_sig is None:
      time_sig = self.pm.time_signature_changes[-1]
    return time_sig

  def _get_ticks_per_bar(self, start):
    time_sig = self._get_time_signature(start)
    quarters_per_bar = 4 * time_sig.numerator / time_sig.denominator
    return self.pm.resolution * quarters_per_bar

  def _get_positions_per_bar(self, start=None, time_sig=None):
    if time_sig is None:
      time_sig = self._get_time_signature(start)
    quarters_per_bar = 4 * time_sig.numerator / time_sig.denominator
    positions_per_bar = int(DEFAULT_POS_PER_QUARTER * quarters_per_bar)
    return positions_per_bar

  def tick_to_position(self, tick):
    return round(tick / self.pm.resolution * DEFAULT_POS_PER_QUARTER)

  # item to event
  def get_remi_events(self):
    events = []
    n_downbeat = 0
    current_chord = None
    current_tempo = None
    for i in range(len(self.groups)):
      bar_st, bar_et = self.groups[i][0], self.groups[i][-1]
      n_downbeat += 1
      positions_per_bar = self._get_positions_per_bar(bar_st)
      if positions_per_bar <= 0:
        raise ValueError('Invalid REMI file: There must be at least 1 position per bar.')

      events.append(Event(
        name=BAR_KEY,
        time=None,
        value='None',
        text='{}'.format(n_downbeat)))

      time_sig = self._get_time_signature(bar_st)
      events.append(Event(
        name=TIME_SIGNATURE_KEY,
        time=None,
        value='{}/{}'.format(time_sig.numerator, time_sig.denominator),
        text='{}/{}'.format(time_sig.numerator, time_sig.denominator)
      ))

      if current_chord is not None:
        events.append(Event(
          name=POSITION_KEY,
          time=0,
          value='{}'.format(0),
          text='{}/{}'.format(1, positions_per_bar)))
        events.append(Event(
          name=CHORD_KEY,
          time=current_chord.start,
          value=current_chord.pitch,
          text='{}'.format(current_chord.pitch)))

      if current_tempo is not None:
        events.append(Event(
          name=POSITION_KEY,
          time=0,
          value='{}'.format(0),
          text='{}/{}'.format(1, positions_per_bar)))
        tempo = current_tempo.pitch
        index = np.argmin(abs(DEFAULT_TEMPO_BINS-tempo))
        events.append(Event(
          name=TEMPO_KEY,
          time=current_tempo.start,
          value=index,
          text='{}/{}'.format(tempo, DEFAULT_TEMPO_BINS[index])))

      quarters_per_bar = 4 * time_sig.numerator / time_sig.denominator
      ticks_per_bar = self.pm.resolution * quarters_per_bar
      flags = np.linspace(bar_st, bar_st + ticks_per_bar, positions_per_bar, endpoint=False)
      for item in self.groups[i][1:-1]:
        # position
        index = np.argmin(abs(flags-item.start))
        pos_event = Event(
          name=POSITION_KEY,
          time=item.start,
          value='{}'.format(index),
          text='{}/{}'.format(index+1, positions_per_bar))

        if item.name == 'Note':
          events.append(pos_event)
          # instrument
          if item.instrument == 'drum':
            name = 'drum'
          else:
            name = pretty_midi.program_to_instrument_name(item.instrument)
          events.append(Event(
            name=INSTRUMENT_KEY,
            time=item.start,
            value=name,
            text='{}'.format(name)))
          # pitch
          events.append(Event(
            name=PITCH_KEY,
            time=item.start,
            value='drum_{}'.format(item.pitch) if name == 'drum' else item.pitch,
            text='{}'.format(pretty_midi.note_number_to_name(item.pitch))))
          # velocity
          velocity_index = np.argmin(abs(DEFAULT_VELOCITY_BINS - item.velocity))
          events.append(Event(
            name=VELOCITY_KEY,
            time=item.start,
            value=velocity_index,
            text='{}/{}'.format(item.velocity, DEFAULT_VELOCITY_BINS[velocity_index])))
          # duration
          duration = self.tick_to_position(item.end - item.start)
          index = np.argmin(abs(DEFAULT_DURATION_BINS-duration))
          events.append(Event(
            name=DURATION_KEY,
            time=item.start,
            value=index,
            text='{}/{}'.format(duration, DEFAULT_DURATION_BINS[index])))
        elif item.name == 'Chord':
          if current_chord is None or item.pitch != current_chord.pitch:
            events.append(pos_event)
            events.append(Event(
              name=CHORD_KEY,
              time=item.start,
              value=item.pitch,
              text='{}'.format(item.pitch)))
            current_chord = item
        elif item.name == 'Tempo':
          if current_tempo is None or item.pitch != current_tempo.pitch:
            events.append(pos_event)
            tempo = item.pitch
            index = np.argmin(abs(DEFAULT_TEMPO_BINS-tempo))
            events.append(Event(
              name=TEMPO_KEY,
              time=item.start,
              value=index,
              text='{}/{}'.format(tempo, DEFAULT_TEMPO_BINS[index])))
            current_tempo = item

    return [{ 'name': e.name, 'value': e.value } for e in events]

  def get_description(self,
                      omit_time_sig=False,
                      omit_instruments=False,
                      omit_chords=False,
                      omit_meta=False):
    events = []
    n_downbeat = 0
    current_chord = None

    for i in range(len(self.groups)):
      bar_st, bar_et = self.groups[i][0], self.groups[i][-1]
      n_downbeat += 1
      time_sig = self._get_time_signature(bar_st)
      positions_per_bar = self._get_positions_per_bar(time_sig=time_sig)
      if positions_per_bar <= 0:
        raise ValueError('Invalid REMI file: There must be at least 1 position in each bar.')

      events.append(Event(
        name=BAR_KEY,
        time=None,
        value='{}'.format(n_downbeat),
        text='{}'.format(n_downbeat)))

      if not omit_time_sig:
        events.append(Event(
          name=TIME_SIGNATURE_KEY,
          time=None,
          value='{}/{}'.format(time_sig.numerator, time_sig.denominator),
          text='{}/{}'.format(time_sig.numerator, time_sig.denominator),
        ))

      if not omit_meta:
        notes = [item for item in self.groups[i][1:-1] if item.name == 'Note']
        n_notes = len(notes)
        velocities = np.array([item.velocity for item in notes])
        pitches = np.array([item.pitch for item in notes])
        durations = np.array([item.end - item.start for item in notes])

        note_density = n_notes/positions_per_bar
        index = np.argmin(abs(DEFAULT_NOTE_DENSITY_BINS-note_density))
        events.append(Event(
          name=NOTE_DENSITY_KEY,
          time=None,
          value=index,
          text='{:.2f}/{:.2f}'.format(note_density, DEFAULT_NOTE_DENSITY_BINS[index])
        ))

        # will be 0 if there's no notes
        mean_velocity = velocities.mean() if len(velocities) > 0 else np.nan
        index = np.argmin(abs(DEFAULT_MEAN_VELOCITY_BINS-mean_velocity))
        events.append(Event(
          name=MEAN_VELOCITY_KEY,
          time=None,
          value=index if mean_velocity != np.nan else 'NaN',
          text='{:.2f}/{:.2f}'.format(mean_velocity, DEFAULT_MEAN_VELOCITY_BINS[index])
        ))

        # will be 0 if there's no notes
        mean_pitch = pitches.mean() if len(pitches) > 0 else np.nan
        index = np.argmin(abs(DEFAULT_MEAN_PITCH_BINS-mean_pitch))
        events.append(Event(
          name=MEAN_PITCH_KEY,
          time=None,
          value=index if mean_pitch != np.nan else 'NaN',
          text='{:.2f}/{:.2f}'.format(mean_pitch, DEFAULT_MEAN_PITCH_BINS[index])
        ))

        # will be 1 if there's no notes
        mean_duration = durations.mean() if len(durations) > 0 else np.nan
        index = np.argmin(abs(DEFAULT_MEAN_DURATION_BINS-mean_duration))
        events.append(Event(
          name=MEAN_DURATION_KEY,
          time=None,
          value=index if mean_duration != np.nan else 'NaN',
          text='{:.2f}/{:.2f}'.format(mean_duration, DEFAULT_MEAN_DURATION_BINS[index])
        ))

      if not omit_instruments:
        instruments = set([item.instrument for item in notes])
        for instrument in instruments:
          instrument = pretty_midi.program_to_instrument_name(instrument) if instrument != 'drum' else 'drum'
          events.append(Event(
            name=INSTRUMENT_KEY,
            time=None,
            value=instrument,
            text=instrument
          ))

      if not omit_chords:
        chords = [item for item in self.groups[i][1:-1] if item.name == 'Chord']
        if len(chords) == 0 and current_chord is not None:
          chords = [current_chord]
        elif len(chords) > 0:
          if chords[0].start > bar_st and current_chord is not None:
            chords.insert(0, current_chord)
          current_chord = chords[-1]

        for chord in chords:
          events.append(Event(
            name=CHORD_KEY,
            time=None,
            value=chord.pitch,
            text='{}'.format(chord.pitch)
          ))

    return [f'{e.name}_{e.value}' for e in events]

#############################################################################################
# WRITE MIDI
#############################################################################################

def remi2midi(events, bpm=120, time_signature=(4, 4), polyphony_limit=16):
  vocab = RemiVocab()

  def _get_time(bar, position, bpm=120, positions_per_bar=48):
    abs_position = bar*positions_per_bar + position
    beat = abs_position / DEFAULT_POS_PER_QUARTER
    return beat/bpm*60

  def _get_time(reference, bar, pos):
    time_sig = reference['time_sig']
    num, denom = time_sig.numerator, time_sig.denominator
    # Quarters per bar, assuming 4 quarters per whole note
    qpb = 4 * num / denom
    ref_pos = reference['pos']
    d_bars = bar - ref_pos[0]
    d_pos = (pos - ref_pos[1]) + d_bars*qpb*DEFAULT_POS_PER_QUARTER
    d_quarters = d_pos / DEFAULT_POS_PER_QUARTER
    # Convert quarters to seconds
    dt = d_quarters / reference['tempo'] * 60
    return reference['time'] + dt

  # time_sigs = [event.split('_')[-1].split('/') for event in events if f"{TIME_SIGNATURE_KEY}_" in event]
  # time_sigs = [(int(num), int(denom)) for num, denom in time_sigs]

  tempo_changes = [event for event in events if f"{TEMPO_KEY}_" in event]
  if len(tempo_changes) > 0:
    bpm = DEFAULT_TEMPO_BINS[int(tempo_changes[0].split('_')[-1])]

  pm = pretty_midi.PrettyMIDI(initial_tempo=bpm)
  num, denom = time_signature
  pm.time_signature_changes.append(pretty_midi.TimeSignature(num, denom, 0))
  current_time_sig = pm.time_signature_changes[0]

  instruments = {}

  # Use implicit timeline: keep track of last tempo/time signature change event
  # and calculate time difference relative to that
  last_tl_event = {
    'time': 0,
    'pos': (0, 0),
    'time_sig': current_time_sig,
    'tempo': bpm
  }

  bar = -1
  n_notes = 0
  polyphony_control = {}
  for i, event in enumerate(events):
    if not bar in polyphony_control:
      polyphony_control[bar] = {}

    if f"{BAR_KEY}_" in events[i]:
      # Next bar is starting
      bar += 1
      polyphony_control[bar] = {}

      if i+1 < len(events) and f"{TIME_SIGNATURE_KEY}_" in events[i+1]:
        num, denom = events[i+1].split('_')[-1].split('/')
        num, denom = int(num), int(denom)
        current_time_sig = last_tl_event['time_sig']
        if num != current_time_sig.numerator or denom != current_time_sig.denominator:
          time = _get_time(last_tl_event, bar, 0)
          time_sig = pretty_midi.TimeSignature(num, denom, time)
          pm.time_signature_changes.append(time_sig)
          last_tl_event['time'] = time
          last_tl_event['pos'] = (bar, 0)
          last_tl_event['time_sig'] = time_sig

    elif i+1 < len(events) and \
        f"{POSITION_KEY}_" in events[i] and \
        f"{TEMPO_KEY}_" in events[i+1]:
      position = int(events[i].split('_')[-1])
      tempo_idx = int(events[i+1].split('_')[-1])
      tempo = DEFAULT_TEMPO_BINS[tempo_idx]

      if tempo != last_tl_event['tempo']:
        time = _get_time(last_tl_event, bar, position)
        last_tl_event['time'] = time
        last_tl_event['pos'] = (bar, position)
        last_tl_event['tempo'] = tempo

    elif i+4 < len(events) and \
        f"{POSITION_KEY}_" in events[i] and \
        f"{INSTRUMENT_KEY}_" in events[i+1] and \
        f"{PITCH_KEY}_" in events[i+2] and \
        f"{VELOCITY_KEY}_" in events[i+3] and \
        f"{DURATION_KEY}_" in events[i+4]:
      # get position
      position = int(events[i].split('_')[-1])
      if not position in polyphony_control[bar]:
        polyphony_control[bar][position] = {}

      # get instrument
      instrument_name = events[i+1].split('_')[-1]
      if instrument_name not in polyphony_control[bar][position]:
        polyphony_control[bar][position][instrument_name] = 0
      elif polyphony_control[bar][position][instrument_name] >= polyphony_limit:
        # If number of notes exceeds polyphony limit, omit this note
        continue

      if instrument_name not in instruments:
        if instrument_name == 'drum':
          instrument = pretty_midi.Instrument(0, is_drum=True)
        else:
          program = pretty_midi.instrument_name_to_program(instrument_name)
          instrument = pretty_midi.Instrument(program)
        instruments[instrument_name] = instrument
      else:
        instrument = instruments[instrument_name]

      # get pitch
      pitch = int(events[i+2].split('_')[-1])
      # get velocity
      velocity_index = int(events[i+3].split('_')[-1])
      velocity = min(127, DEFAULT_VELOCITY_BINS[velocity_index])
      # get duration
      duration_index = int(events[i+4].split('_')[-1])
      duration = DEFAULT_DURATION_BINS[duration_index]
      # create not and add to instrument
      start = _get_time(last_tl_event, bar, position)
      end = _get_time(last_tl_event, bar, position + duration)
      note = pretty_midi.Note(velocity=velocity,
                              pitch=pitch,
                              start=start,
                              end=end)
      instrument.notes.append(note)
      n_notes += 1
      polyphony_control[bar][position][instrument_name] += 1

  for instrument in instruments.values():
    pm.instruments.append(instrument)
  return pm

max_stars_repo_name: dvruette/MuseMorphose

max_issues_repo_name: dvruette/MuseMorphose

max_forks_repo_name: dvruette/MuseMorphose


text:
import numpy as np
import cv2

# Identify pixels above the threshold
# Threshold of RGB > 160 does a nice job of identifying ground pixels only
def color_thresh(img, rgb_thresh=(160, 160, 160), rock_thresh_lower=(20, 100, 100), rock_thresh_upper=(30, 255, 255)):
    # Terrain
    color_select_terrain = np.zeros_like(img[:, :, 0])
    thresh_terrain = (img[:, :, 0] > rgb_thresh[0]) & (img[:, :, 1] > rgb_thresh[1]) & (img[:, :, 2] > rgb_thresh[2])
    color_select_terrain[thresh_terrain] = 1

    # Obstacles
    color_select_obstacles = np.zeros_like(img[:, :, 0])
    thresh_obstacles = (img[:, :, 0] < rgb_thresh[0]) & (img[:, :, 1] < rgb_thresh[1]) & (img[:, :, 2] < rgb_thresh[2])
    color_select_obstacles[thresh_obstacles] = 1

    # Rocks
    color_select_rock = np.zeros_like(img[:, :, 0])
    hsv = cv2.cvtColor(img, cv2.COLOR_RGB2HSV)
    thresh_rock = (hsv[:, :, 0] < rock_thresh_upper[0]) & (hsv[:, :, 0] > rock_thresh_lower[0]) & \
                  (hsv[:, :, 1] < rock_thresh_upper[1]) & (hsv[:, :, 1] > rock_thresh_lower[1]) & \
                  (hsv[:, :, 2] < rock_thresh_upper[2]) & (hsv[:, :, 2] > rock_thresh_lower[2])
    color_select_rock[thresh_rock] = 1
    return color_select_terrain, color_select_obstacles, color_select_rock

# Define a function to convert from image coords to rover coords
def rover_coords(binary_img):
    # Identify nonzero pixels
    ypos, xpos = binary_img.nonzero()
    # Calculate pixel positions with reference to the rover position being at the
    # center bottom of the image.
    x_pixel = -(ypos - binary_img.shape[0]).astype(np.float)
    y_pixel = -(xpos - binary_img.shape[1] / 2).astype(np.float)
    return x_pixel, y_pixel

# Define a function to convert to radial coords in rover space
def to_polar_coords(x_pixel, y_pixel):
    # Convert (x_pixel, y_pixel) to (distance, angle)
    # in polar coordinates in rover space
    # Calculate distance to each pixel
    dist = np.sqrt(x_pixel ** 2 + y_pixel ** 2)
    # Calculate angle away from vertical for each pixel
    angles = np.arctan2(y_pixel, x_pixel)
    return dist, angles

# Define a function to map rover space pixels to world space
def rotate_pix(xpix, ypix, yaw):
    # Convert yaw to radians
    yaw_rad = yaw * np.pi / 180
    xpix_rotated = (xpix * np.cos(yaw_rad)) - (ypix * np.sin(yaw_rad))
    ypix_rotated = (xpix * np.sin(yaw_rad)) + (ypix * np.cos(yaw_rad))
    # Return the result
    return xpix_rotated, ypix_rotated

def translate_pix(xpix_rot, ypix_rot, xpos, ypos, scale):
    # Apply a scaling and a translation
    xpix_translated = (xpix_rot / scale) + xpos
    ypix_translated = (ypix_rot / scale) + ypos
    # Return the result
    return xpix_translated, ypix_translated

# Define a function to apply rotation and translation (and clipping)
# Once you define the two functions above this function should work
def pix_to_world(xpix, ypix, xpos, ypos, yaw, world_size, scale):
    # Apply rotation
    xpix_rot, ypix_rot = rotate_pix(xpix, ypix, yaw)
    # Apply translation
    xpix_tran, ypix_tran = translate_pix(xpix_rot, ypix_rot, xpos, ypos, scale)
    # Perform rotation, translation and clipping all at once
    x_pix_world = np.clip(np.int_(xpix_tran), 0, world_size - 1)
    y_pix_world = np.clip(np.int_(ypix_tran), 0, world_size - 1)
    # Return the result
    return x_pix_world, y_pix_world

# Define a function to perform a perspective transform
def perspect_transform(img, src, dst):
    M = cv2.getPerspectiveTransform(src, dst)
    warped = cv2.warpPerspective(img, M, (img.shape[1], img.shape[0]))  # keep same size as input image
    return warped

# Apply the above functions in succession and update the Rover state accordingly
def perception_step(Rover):
    # Perform perception steps to update Rover()
    # NOTE: camera image is coming to you in Rover.img

    # 1) Define source and destination points for perspective transform
    dst_size = 5
    bottom_offset = 6
    source = np.float32([[14, 140], [301, 140], [200, 96], [118, 96]])
    destination = np.float32([[Rover.img.shape[1] / 2 - dst_size, Rover.img.shape[0] - bottom_offset],
                              [Rover.img.shape[1] / 2 + dst_size, Rover.img.shape[0] - bottom_offset],
                              [Rover.img.shape[1] / 2 + dst_size, Rover.img.shape[0] - 2 * dst_size - bottom_offset],
                              [Rover.img.shape[1] / 2 - dst_size, Rover.img.shape[0] - 2 * dst_size - bottom_offset]])

    # 2) Apply perspective transform
    warped = perspect_transform(Rover.img, source, destination)

    # 3) Apply color threshold to identify navigable terrain/obstacles/rock samples
    threshed_terrain, threshed_obstacles, threshed_rock = color_thresh(warped)

    # 4) Update Rover.vision_image (this will be displayed on left side of screen)
    Rover.vision_image[:, :, 0] = threshed_obstacles * 255
    Rover.vision_image[:, :, 1] = threshed_rock * 255
    Rover.vision_image[:, :, 2] = threshed_terrain * 255

    # 5) Convert map image pixel values to rover-centric coords
    xpix_terrain, ypix_terrain = rover_coords(threshed_terrain)
    xpix_obstacles, ypix_obstacles = rover_coords(threshed_obstacles)
    xpix_rock, ypix_rock = rover_coords(threshed_rock)

    # 6) Convert rover-centric pixel values to world coordinates
    scale = 10
    x_world_terrain, y_world_terrain = pix_to_world(xpix_terrain, ypix_terrain, Rover.pos[0], Rover.pos[1], Rover.yaw,
                                                    Rover.worldmap.shape[0], scale)
    x_world_obstacles, y_world_obstacles = pix_to_world(xpix_obstacles, ypix_obstacles, Rover.pos[0], Rover.pos[1],
                                                        Rover.yaw, Rover.worldmap.shape[0], scale)
    x_world_rock, y_world_rock = pix_to_world(xpix_rock, ypix_rock, Rover.pos[0], Rover.pos[1], Rover.yaw,
                                              Rover.worldmap.shape[0], scale)

    # 7) Update Rover worldmap (to be displayed on right side of screen)
    Rover.worldmap[y_world_obstacles, x_world_obstacles, 0] += 10  # Favor navigable terrain
    Rover.worldmap[y_world_rock, x_world_rock, 1] += 1
    Rover.worldmap[y_world_terrain, x_world_terrain, 2] += 1

    # 8) Convert rover-centric pixel positions to polar coordinates
    # Update Rover pixel distances and angles
    # Rover.nav_dists = rover_centric_pixel_distances
    # Rover.nav_angles = rover_centric_angles
    Rover.nav_dists, Rover.nav_angles = to_polar_coords(xpix_terrain, ypix_terrain)

    return Rover


text:
import os
import time
from typing import Generator

import docker
import numpy as np
import pytest

from tempo.aio.model import Model
from tempo.aio.pipeline import Pipeline
from tempo.aio.utils import model, pipeline
from tempo.seldon import SeldonDockerRuntime
from tempo.serve.metadata import DockerOptions, ModelFramework
from tempo.serve.model import Model as _Model
from tempo.serve.pipeline import PipelineModels

@pytest.fixture
def runtime() -> SeldonDockerRuntime:
    return SeldonDockerRuntime(DockerOptions())

@pytest.fixture
def sklearn_model(sklearn_model: _Model, runtime: SeldonDockerRuntime) -> Generator[Model, None, None]:
    model = Model(
        name=sklearn_model.details.name,
        platform=sklearn_model.details.platform,
        uri=sklearn_model.details.uri,
        local_folder=sklearn_model.details.local_folder,
        protocol=sklearn_model.model_spec.protocol,
        runtime_options=sklearn_model.model_spec.runtime_options,
    )

    runtime.deploy(model)
    runtime.wait_ready(model)
    time.sleep(2)

    yield model

    try:
        runtime.undeploy(model)
    except docker.errors.NotFound:
        # TODO: Should undeploy be idempotent as well?
        # Ignore if the model has already been undeployed
        pass

@pytest.fixture
def custom_model() -> Model:
    @model(name="custom-model", platform=ModelFramework.Custom)
    async def _custom_model(payload: np.ndarray) -> np.ndarray:
        return payload.sum(keepdims=True)

    return _custom_model

@pytest.fixture
def inference_pipeline(
    sklearn_model, runtime: SeldonDockerRuntime, pipeline_conda_yaml: str
) -> Generator[Pipeline, None, None]:
    @pipeline(
        name="inference-pipeline",
        models=PipelineModels(sklearn=sklearn_model),
        local_folder=os.path.dirname(pipeline_conda_yaml),
    )
    async def _pipeline(payload: np.ndarray) -> np.ndarray:
        res1 = await _pipeline.models.sklearn(payload)
        if res1[0][0] > 0.7:
            return res1

        return res1.sum(keepdims=True)

    _pipeline.save(save_env=True)
    runtime.deploy(_pipeline)
    runtime.wait_ready(_pipeline)
    time.sleep(8)

    yield _pipeline

    try:
        runtime.undeploy(_pipeline)
    except docker.errors.NotFound:
        # TODO: Should undeploy be idempotent as well?
        # Ignore if the model has already been undeployed
        pass

max_stars_repo_name: outerbounds/tempo

max_issues_repo_name: outerbounds/tempo

max_forks_repo_name: outerbounds/tempo


text:
# Script to resample a raster to a larger pixel size using byte sequences.

import os
import numpy as np
from osgeo import gdal

# Don't forget to change the folder.
os.chdir(r'D:\osgeopy-data\Landsat\Washington')

# Open the input raster.
in_ds = gdal.Open('nat_color.tif')

# Computer the number of output rows and columns (half the input numbers
# because we're making the pixels twice as big).
out_rows = int(in_ds.RasterYSize / 2)
out_columns = int(in_ds.RasterXSize / 2)
num_bands = in_ds.RasterCount

# Create the output raster using the computed dimensions.
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create('nat_color_resampled.tif',
    out_columns, out_rows, num_bands)

# Change the geotransform so it reflects the larger cell size before
# setting it onto the output.
out_ds.SetProjection(in_ds.GetProjection())
geotransform = list(in_ds.GetGeoTransform())
geotransform[1] *= 2
geotransform[5] *= 2
out_ds.SetGeoTransform(geotransform)

# Read in the data for all bands, but have gdal resample it so that it has
# the specified number of rows and columns instead of the numbers that the
# input has. This effectively resizes the pixels.
data = in_ds.ReadRaster(
    buf_xsize=out_columns, buf_ysize=out_rows)

# Write the data to the output raster.
out_ds.WriteRaster(0, 0, out_columns, out_rows, data)

# Compute statistics and build overviews.
out_ds.FlushCache()
for i in range(num_bands):
    out_ds.GetRasterBand(i + 1).ComputeStatistics(False)
out_ds.BuildOverviews('average', [2, 4, 8, 16])

del out_ds

max_stars_repo_name: hohsieh/osgeopy-code

max_issues_repo_name: sthagen/osgeopy-code

max_forks_repo_name: sthagen/osgeopy-code


text:
import networkx as nx
from noisy_graphs.noisy_graph import NoisyGraph

def create_data_path_file(data_path: str):
    header = "exp_name,sigma_mean,sigma_variance,uncertainty_mean,uncertainty_variance,dc_distance," \
             "dc_correlation,dc_mean_change,bc_distance,bc_correlation,bc_mean_change,cc_distance," \
             "cc_correlation,cc_mean_change,ec_distance,ec_correlation,ec_mean_change\n"
    f = open(data_path, "w")
    f.write(header)
    f.close()

def perform_experiment(original_graph: nx.Graph, ftrp: float, exp_name: str, data_path: str):
    print(exp_name)

    # removing graph isolates
    # NoisyGraph is not intended to deal with isolated nodes
    original_graph.remove_nodes_from(list(nx.isolates(original_graph)))

    # constructing noisy graph
    noisy_graph = NoisyGraph(ftrp=ftrp)
    noisy_graph.construct_graph(original_graph)

    # algorithm compliance
    sigma_mean, sigma_variance = noisy_graph.get_sigmas_profile()

    # uncertainty
    uncertainty_mean, uncertainty_variance = noisy_graph.get_uncertainty_profile()

    # centrality_metrics
    dc_distance, dc_correlation, dc_mean_change = noisy_graph.degree_centrality_profile(original_graph)
    bc_distance, bc_correlation, bc_mean_change = noisy_graph.betweenness_profile(original_graph)
    cc_distance, cc_correlation, cc_mean_change = noisy_graph.closeness_profile(original_graph)
    ec_distance, ec_correlation, ec_mean_change = noisy_graph.eigenvector_centrality_profile(original_graph)

    # create raw file with original edges and noisy edges
    f = open(f"raw_data/{exp_name[:2]}.txt", "a")
    f.write(f"Experimental conditions: {exp_name}\n")
    f.write(f"Real edges: {noisy_graph.edges_if(real=True)}\n")
    f.write(f"Fake edges: {noisy_graph.edges_if(real=False)}\n\n")
    f.close()

    # add result to csv
    result = f"{exp_name},{sigma_mean},{sigma_variance},{uncertainty_mean},{uncertainty_variance},"
    result += f"{dc_distance},{dc_correlation},{dc_mean_change},"
    result += f"{bc_distance},{bc_correlation},{bc_mean_change},"
    result += f"{cc_distance},{cc_correlation},{cc_mean_change},"
    result += f"{ec_distance},{ec_correlation},{ec_mean_change}"
    result += "\n"

    f = open(data_path, "a")
    f.write(result)
    f.close()

max_stars_repo_name: lechugaa/negative-graphs

max_issues_repo_name: lechugaa/negative-graphs

max_forks_repo_name: lechugaa/negative-graphs


text:
import sys
import unittest as ut

import numpy as np
import numpy.testing as nt
import xarray as xr

# Import from directory structure if coverage test, or from installed
# packages otherwise
if "--cov" in str(sys.argv):
    from src.geocat.f2py import moc_globe_atl
else:
    from geocat.f2py import moc_globe_atl

# Dimensions
nyaux = 3  # nyaux = lat_aux_grid->shape[0]
kdep = 4  # kdep  = a_wvel/a_bolus/a_submeso->shape[0]
nlat = 5  # nlat  = a_wvel/a_bolus/a_submeso->shape[1] AND tlat->shape[0] AND rmlak->shape[1]
mlon = 6  # mlon  = a_wvel/a_bolus/a_submeso->shape[2] AND tlat->shape[1] AND rmlak->shape[2]

kdepnyaux2 = 2 * kdep * nyaux
nlatmlon = nlat * mlon
kdepnlatmlon = kdep * nlatmlon

# Generate arbitrary data
tmp_a = np.linspace(1, kdepnlatmlon, num=kdepnlatmlon).reshape(
    (kdep, nlat, mlon))

# nan input
tmp_a_nan = tmp_a.copy()
tmp_a_nan[0, 1, 1] = np.nan
tmp_a_nan[0, 2, 3] = np.nan
tmp_a_nan[1, 1, 1] = np.nan
tmp_a_nan[1, 2, 3] = np.nan

# nan input
tmp_a_msg = tmp_a_nan.copy()
tmp_a_msg[np.isnan(tmp_a_msg)] = -99

# Generate test data
lat_aux_grid = np.asarray([5, 10, 15])
rmlak = np.ones((2, nlat, mlon))
t_lat = np.linspace(1, nlatmlon, num=nlatmlon).reshape((nlat, mlon))

# Generate expected NCL outputs
ncl_truth = [
    0.00, 35.00, 60.00, 0.00, 185.00, 210.00, 0.00, 335.00, 360.00, 0.00,
    485.00, 510.00, 0.00, 35.00, 60.00, 0.00, 185.00, 210.00, 0.00, 335.00,
    360.00, 0.00, 485.00, 510.00, 0.00, 35.00, 60.00, 0.00, 185.00, 210.00,
    0.00, 335.00, 360.00, 0.00, 485.00, 510.00, 0.00, 35.00, 60.00, 0.00,
    185.00, 210.00, 0.00, 335.00, 360.00, 0.00, 485.00, 510.00, 0.00, 35.00,
    60.00, 0.00, 185.00, 210.00, 0.00, 335.00, 360.00, 0.00, 485.00, 510.00,
    0.00, 35.00, 60.00, 0.00, 185.00, 210.00, 0.00, 335.00, 360.00, 0.00,
    485.00, 510.00
]
ncl_truth = np.reshape(ncl_truth, (3, 2, kdep, nyaux))

ncl_truth_msg = [
    0.00, 27.00, 60.00, 0.00, 147.00, 210.00, 0.00, 335.00, 360.00, 0.00,
    485.00, 510.00, 0.00, 27.00, 60.00, 0.00, 147.00, 210.00, 0.00, 335.00,
    360.00, 0.00, 485.00, 510.00, 0.00, 27.00, 60.00, 0.00, 147.00, 210.00,
    0.00, 335.00, 360.00, 0.00, 485.00, 510.00, 0.00, 27.00, 60.00, 0.00,
    147.00, 210.00, 0.00, 335.00, 360.00, 0.00, 485.00, 510.00, 0.00, 27.00,
    60.00, 0.00, 147.00, 210.00, 0.00, 335.00, 360.00, 0.00, 485.00, 510.00,
    0.00, 27.00, 60.00, 0.00, 147.00, 210.00, 0.00, 335.00, 360.00, 0.00,
    485.00, 510.00
]
ncl_truth_msg = np.reshape(ncl_truth_msg, (3, 2, kdep, nyaux))

class Test_Moc_Globe_Atl(ut.TestCase):

    def test_moc_globe_atl_float64(self):

        out_arr = moc_globe_atl(lat_aux_grid, tmp_a.astype(np.float64),
                                tmp_a.astype(np.float64),
                                tmp_a.astype(np.float64), t_lat, rmlak)

        nt.assert_array_almost_equal(ncl_truth, out_arr)
        nt.assert_equal((3, 2, kdep, nyaux), out_arr.shape)

    def test_moc_globe_atl_float64_xr(self):

        out_arr = moc_globe_atl(xr.DataArray(lat_aux_grid),
                                xr.DataArray(tmp_a.astype(np.float64)),
                                xr.DataArray(tmp_a.astype(np.float64)),
                                xr.DataArray(tmp_a.astype(np.float64)),
                                xr.DataArray(t_lat), xr.DataArray(rmlak))

        nt.assert_array_almost_equal(ncl_truth, out_arr)
        nt.assert_equal((3, 2, kdep, nyaux), out_arr.shape)

    def test_moc_globe_atl_float32(self):

        out_arr = moc_globe_atl(lat_aux_grid, tmp_a.astype(np.float32),
                                tmp_a.astype(np.float32),
                                tmp_a.astype(np.float32), t_lat, rmlak)

        nt.assert_array_almost_equal(ncl_truth, out_arr)
        nt.assert_equal((3, 2, kdep, nyaux), out_arr.shape)

    def test_moc_globe_atl_int64(self):

        out_arr = moc_globe_atl(lat_aux_grid, tmp_a.astype(np.int64),
                                tmp_a.astype(np.int64), tmp_a.astype(np.int64),
                                t_lat, rmlak)

        nt.assert_array_almost_equal(ncl_truth, out_arr)
        nt.assert_equal((3, 2, kdep, nyaux), out_arr.shape)

    def test_moc_globe_atl_msg_99(self):

        out_arr = moc_globe_atl(lat_aux_grid,
                                tmp_a_msg,
                                tmp_a_msg,
                                tmp_a_msg,
                                t_lat,
                                rmlak,
                                msg=-99)

        nt.assert_array_almost_equal(ncl_truth_msg, out_arr)
        nt.assert_equal((3, 2, kdep, nyaux), out_arr.shape)

    def test_moc_globe_atl_msg_nan(self):

        out_arr = moc_globe_atl(lat_aux_grid,
                                tmp_a_nan,
                                tmp_a_nan,
                                tmp_a_nan,
                                t_lat,
                                rmlak,
                                msg=np.nan)

        nt.assert_array_almost_equal(ncl_truth_msg, out_arr)
        nt.assert_equal((3, 2, kdep, nyaux), out_arr.shape)


text:
import multiprocessing as mp

import numpy as np
import pytest

from Starfish.grid_tools import (
    chunk_list,
    air_to_vacuum,
    vacuum_to_air,
    vacuum_to_air_SLOAN,
    idl_float,
)

class TestChunking:
    def test_chunk_list_shape(self, grid_points):
        chunked = chunk_list(grid_points)
        assert len(chunked) == mp.cpu_count()
        chunked = chunk_list(grid_points, 3)
        assert len(chunked) == 3

class TestWavelengthUtils:
    def test_air_to_vacuum(self):
        wavelengths = np.linspace(1e4, 5e4, 1000)
        outputs = air_to_vacuum(wavelengths)
        assert len(outputs) == len(wavelengths)

    def test_vacuum_to_air(self):
        wavelengths = np.linspace(1e4, 5e4, 1000)
        outputs = vacuum_to_air(wavelengths)
        assert len(outputs) == len(wavelengths)

    def test_vacuum_to_air_SLOAN(self):
        wavelengths = np.linspace(1e4, 5e4, 1000)
        outputs = vacuum_to_air_SLOAN(wavelengths)
        assert len(outputs) == len(wavelengths)

    @pytest.mark.parametrize(
        "wavelengths", [np.linspace(1e4, 5e4, 1000), np.linspace(1e5, 5e5, 1000)]
    )
    def test_atv_vta_regression(self, wavelengths):
        np.testing.assert_array_almost_equal(
            wavelengths, vacuum_to_air(air_to_vacuum(wavelengths)), 2
        )

    @pytest.mark.parametrize(
        "wavelengths", [np.linspace(1e4, 5e4, 1000), np.linspace(1e5, 5e5, 1000)]
    )
    def test_atv_vta_sloan_regression(self, wavelengths):
        np.testing.assert_array_almost_equal(
            wavelengths, vacuum_to_air_SLOAN(air_to_vacuum(wavelengths)), 0
        )

@pytest.mark.parametrize(
    "idl, num", [("1D4", 1e4), ("1.0", 1.0), ("1D-4", 1e-4), ("1d0", 1.0)]
)
def test_idl_float(idl, num):
    np.testing.assert_almost_equal(idl_float(idl), num)

max_stars_repo_name: spencerhurt/Starfish

max_issues_repo_name: spencerhurt/Starfish

max_forks_repo_name: spencerhurt/Starfish


text:
import pytest
from pytest import approx
from lazydiff.vars import Var
from lazydiff import ops
from lazydiff import regression

import numpy as np
from sklearn.datasets import make_regression
from sklearn.linear_model import LinearRegression, Lasso, Ridge, ElasticNet
from sklearn.preprocessing import PolynomialFeatures

np.random.seed(1)

dim = 1
X,y,true_coef = make_regression(n_samples = 100, n_features = dim, n_informative = dim, bias = 10, \
                                coef = True, noise = 0, random_state=1)
def test_MSE():
    m = Var(np.ones(X.shape[1]))
    b = Var(0)
    loss = regression.MSE(X, y, m, b)
    assert loss.val == approx(np.sum((X.sum(axis=1)-y)**2)/X.shape[0])

def test_MSE_regularized():
    m = Var(np.ones(X.shape[1]))
    b = Var(0)
    loss = regression.MSE_regularized(X, y, m, b)
    numeric_loss = np.sum((X.sum(axis=1)-y)**2)/(2*X.shape[0]) + np.linalg.norm(m.val,1)
    assert loss.val == approx(numeric_loss)

def test_lasso_loss():
    m = Var(np.ones(X.shape[1]))
    b = Var(0)
    loss = regression.lasso_loss(X, y, m, b)
    numeric_loss = np.sum((X.sum(axis=1)-y)**2)/(2*X.shape[0]) + np.linalg.norm(m.val,1)
    assert loss.val == approx(numeric_loss)

def test_ridge_loss():
    m = Var(np.ones(X.shape[1]))
    b = Var(0)
    loss = regression.ridge_loss(X,y,m,b)
    numeric_loss = np.sum((X.sum(axis=1)-y)**2) + np.linalg.norm(m.val)**2
    assert loss.val == approx(numeric_loss)

def test_elastic_loss():
    m = Var(np.ones(X.shape[1]))
    b = Var(0)
    loss = regression.elastic_loss(X, y, m, b)
    numeric_loss = np.sum((X.sum(axis=1)-y)**2)/(2*X.shape[0]) + 0.5*np.linalg.norm(m.val,1) + 0.5*0.5*np.linalg.norm(m.val)**2
    assert loss.val == approx(numeric_loss)

def test_gradient_descent():
    m = Var(np.ones(X.shape[1]))
    b = Var(0)
    new_m, new_b, new_loss = regression.gradient_descent(X, y, regression.MSE, m, b)
    numeric_mse = np.sum((X.sum(axis=1)-y)**2)/(X.shape[0])
    # manually doing gradient descent
    n = X.shape[0]
    m_grad = 2/n*np.sum(-X.reshape(-1)*(y-((m.val*X).sum(axis=1))))
    manual_m = m.val-0.1*m_grad
    b_grad = 2/n*np.sum(-1*(y-((m.val*X).sum(axis=1)+b.val)))
    manual_b = b.val-0.1*b_grad

    assert new_loss.val == approx(numeric_mse)
    assert new_m.val == approx(manual_m, abs = 1e-3)
    assert new_b.val == approx(manual_b, abs = 1e-3)

def test_lasso():
    m_lasso = Var(np.ones(X.shape[1]))
    b_lasso = Var(0)
    earlyStop = 0 #1e-8
    forward = True
    epochs = 300
    lr = 0.1
    m_lasso, b_lasso, loss = regression.iterative_regression(X, y, m_lasso, b_lasso, regression.lasso_loss, lr, epochs,
                                                        earlyStop, forward)
    clf_l1 = Lasso().fit(X,y)
    assert m_lasso.val == approx(clf_l1.coef_, abs=1e-3)
    assert b_lasso.val == approx(clf_l1.intercept_, abs=1e-3)

def test_ridge():
    m_ridge = Var(np.ones(X.shape[1]))
    b_ridge = Var(0)
    earlyStop = 0 #1e-8
    forward = False
    epochs = 300
    lr = 0.001
    m_ridge, b_ridge, loss = regression.iterative_regression(X, y, m_ridge, b_ridge, regression.ridge_loss, lr, epochs,
                                                        earlyStop, forward)
    clf_l2 = Ridge().fit(X,y)
    assert m_ridge.val == approx(clf_l2.coef_, abs=1e-3)
    assert b_ridge.val == approx(clf_l2.intercept_, abs=1e-3)

def test_elastic():
    m_el = Var(np.ones(X.shape[1]))
    b_el = Var(0)
    earlyStop = 0 #1e-8
    forward = False
    epochs = 300
    lr = 0.1
    m_el, b_el, loss = regression.iterative_regression(X, y, m_el, b_el, regression.elastic_loss, lr, epochs,
                                                        earlyStop, forward)
    clf_el = ElasticNet().fit(X,y)
    assert m_el.val == approx(clf_el.coef_, abs=1e-3)
    assert b_el.val == approx(clf_el.intercept_, abs=1e-3)

def test_polynomial():
    X = np.linspace(-10,10,20).reshape(-1,1)
    y = X.reshape(-1)**2
    X = PolynomialFeatures(2, include_bias = False).fit_transform(X)
    m = Var(np.random.rand(X.shape[1]), seed = 1.0)
    b = Var(0)
    earlyStop = 0
    forward = False
    m, b, loss = regression.iterative_regression(X, y, m, b, regression.MSE, 0.0001, 10000,
                                                        earlyStop, forward)
    predict = np.sum(m.val*X, axis=1)+b.val
    assert predict == approx(y, abs=1e-2)

def test_store():
    m = Var(np.ones(X.shape[1]))
    b = Var(0)
    earlyStop = 1e-8
    forward = True
    history = {}
    m, b, loss = regression.iterative_regression(X, y, m, b, regression.MSE, 0.1, 100,
                                                    earlyStop, forward, history)
    clf = LinearRegression().fit(X,y)
    assert m.val == approx(clf.coef_, abs=1e-3)
    assert b.val == approx(clf.intercept_, abs=1e-3)
    assert history['m'][-1] == m.val
    assert history['b'][-1] == b.val
    assert history['loss'][-1] == loss.val