text:
import os
import pickle
from functools import lru_cache
import importlib.util

import numpy as np
from util import cvimage as Image

import config

class ResourceArchiveIndex:
    def __init__(self, archive, archive_path):
        self.archive = archive
        self.archive_path = archive_path
    def open(self):
        return self.archive.open(self.archive_path, 'r')
    def __hash__(self):
        return hash((ResourceArchiveIndex, self.archive, self.archive_path))

class FileSystemIndex:
    def __init__(self, path):
        self.path = path
    def open(self):
        return open(self.path, 'rb')
    def __hash__(self):
        return hash((FileSystemIndex, self.path))

if config.use_archived_resources:
    import zipfile
    root = 'resources/imgreco'
    archive = zipfile.ZipFile(open(config.resource_archive, 'rb'), 'r')

    filelist = archive.namelist()

    def get_path(names):
        return '/'.join([root, *names])

    def _open_file(path):
        return archive.open(path)

    def _get_index(names):
        archive_path = get_path(names)
        return ResourceArchiveIndex(archive, archive_path)

    def get_entries(base):
        prefix = 'resources/imgreco/' + base + '/'
        dirs = []
        files = []
        for name in filelist:
            if name.startswith(prefix):
                name = name[len(prefix):]
                if len(name) == 0:
                    continue
                elif name[-1] == '/':
                    dirs.append(name[:-1])
                elif '/' in name:
                    continue
                else:
                    files.append(name)
        return dirs, files
else:
    root = os.path.join(config.resource_root, 'imgreco')

    def get_path(names):
        return os.path.join(root, *names)

    def _get_index(names):
        path = get_path(names)
        return FileSystemIndex(os.path.join(root, path))

    def _open_file(path):
        return open(path, 'rb')

    def get_entries(base):
        findroot = get_path(base.split('/'))
        for _, dirs, files in os.walk(findroot):
            return (dirs, files)
        return ([], [])

def resolve(respath):
    names = respath.split('/')
    return _get_index(names)

def open_file(respath):
    if hasattr(respath, 'open'):
        return respath.open()
    return resolve(respath).open()

def load_image(name, mode=None):
    im = Image.open(open_file(name))
    if mode is not None and im.mode != mode:
        im = im.convert(mode)
    return im

@lru_cache(maxsize=None)
def load_image_cached(name, mode=None):
    return load_image(name, mode)

def load_image_as_ndarray(name):
    return np.asarray(load_image(name))

def load_pickle(name):
    with open_file(name) as f:
        result = pickle.load(f)
    return result

def load_minireco_model(name, filter_chars=None):
    model = load_pickle(name)
    if filter_chars is not None:
        model['data'] = [x for x in model['data'] if x[0] in filter_chars]
        model['chars'] = [x[0] for x in model['data']]
    return model

max_stars_repo_name: HTLXMC/ArknightsAutoHelper

max_issues_repo_name: HTLXMC/ArknightsAutoHelper

max_forks_repo_name: HTLXMC/ArknightsAutoHelper


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright © 2020 Michael J. Hayford
"""Dashboards constructed with ipywidgets.

.. Created on Tue Apr 28 21:15:20 2020

.. codeauthor: Michael J. Hayford
"""

import numpy as np
import ipywidgets as widgets

def create_focus_dashboard(figs, ray_data_items, foc, dfoc_rng, shift_rng,
                           on_axis_pt, continuous_update=True):
    defocus = widgets.FloatSlider(min=-dfoc_rng, max=+dfoc_rng,
                                  step=.01*dfoc_rng,
                                  description='defocus', value=foc,
                                  readout_format='.4f',
                                  continuous_update=continuous_update)
    x_shift = widgets.FloatSlider(min=-shift_rng, max=+shift_rng,
                                  step=.01*shift_rng,
                                  description='x shift', value=on_axis_pt[0],
                                  readout_format='.4f',
                                  continuous_update=continuous_update)
    y_shift = widgets.FloatSlider(min=-shift_rng, max=+shift_rng,
                                  step=.01*shift_rng,
                                  description='y shift', value=on_axis_pt[1],
                                  readout_format='.4f',
                                  continuous_update=continuous_update)

    def slider_update(change):
        dfoc_val = defocus.value
        dx = x_shift.value
        dy = y_shift.value

        # apply changes to fans and grids
        for ray_data in ray_data_items:
            ray_data.foc = dfoc_val
            ray_data.image_pt_2d = np.array([dx, dy])
            ray_data.update_data(build='update')

        # update and plot results
        for fig in figs:
            fig.clf()
            fig.plot()

    defocus.observe(slider_update, names='value')
    x_shift.observe(slider_update, names='value')
    y_shift.observe(slider_update, names='value')

    return defocus, x_shift, y_shift

def create_mirror_tilt_dashboard(mirror, app_mgr, figs, ray_data_items,
                                 foc, tilt_rng, shift_rng, oa_ray=None,
                                 continuous_update=True):
    tilt0 = mirror.decenter.euler[0]
    alpha_tilt = widgets.FloatSlider(min=tilt0-tilt_rng, max=tilt0+tilt_rng,
                                     step=.05*tilt_rng,
                                     description='alpha tilt', value=tilt0,
                                     readout_format='.4f',
                                     continuous_update=continuous_update)
    x_shift = widgets.FloatSlider(min=-shift_rng, max=+shift_rng,
                                  step=.05*shift_rng, description='x shift',
                                  value=0., readout_format='.4f',
                                  continuous_update=continuous_update)
    y_shift = widgets.FloatSlider(min=-shift_rng, max=+shift_rng,
                                  step=.05*shift_rng, description='y shift',
                                  value=0., readout_format='.4f',
                                  continuous_update=continuous_update)

    def make_slider_callback(changer):
        def slider_update(change):
            changer.set(change['new'])

            app_mgr.refresh_gui(build='update')

            if oa_ray is not None:
                oa_ray.update_data(build='rebuild')

            # apply changes to fans and grids
            for ray_data in ray_data_items:
                # if oa_ray is not None:
                #     ray_data.image_pt_2d = oa_ray.t_abr
                ray_data.update_data(build='rebuild')

            # update and plot results
            for fig in figs:
                fig.clf()
                fig.plot()

        return slider_update

    alpha_tilt.observe(make_slider_callback(AttrChanger(mirror.decenter,
                                                        'euler', index=0)),
                       names='value')
    x_shift.observe(make_slider_callback(AttrChanger(mirror.decenter,
                                                     'dec', index=0)),
                    names='value')
    y_shift.observe(make_slider_callback(AttrChanger(mirror.decenter,
                                                     'dec', index=1)),
                    names='value')

    return alpha_tilt, x_shift, y_shift

class AttrChanger():
    """Changer built on an object/attribute pair. """
    def __init__(self, obj, attr, index=None):
        self.object = obj
        self.attr = attr
        self.index = index

    def get(self):
        value = getattr(self.object, self.attr, None)
        #print('AttrChanger.get:', self.attr, value, self.index)
        value = value if self.index is None else value[self.index]
        return value

    def set(self, value):
        if self.index is not None:
            seq = getattr(self.object, self.attr, None)
            seq[self.index] = value
            value = seq
        setattr(self.object, self.attr, value)
        #print('AttrChanger.set:', self.attr, value)

ext: py

max_stars_repo_name: ajeddeloh/ray-optics

max_issues_repo_name: ajeddeloh/ray-optics

max_forks_repo_name: ajeddeloh/ray-optics

reason: import numpy


text:
import numpy as np
from jacobi import jacobi
from pascalindex import pascalindex2d, pascalindex3d
from numpy.polynomial import Polynomial as poly

""" It is important to note that the koornwinder polynomials are not the shape functions (nodal basis functions). They are merely a 2D equivalent of the legendre orthogonal basis functions in 1D.
    In 1D: legendre polynomials formed the columns in the vandermonde matrix, so the columns were orthogonal. The nodal basis functions were created by solving VL=I, since the shape functions sampled at the nodes yields the identity matrix.
    The equivalent to skipping the legendre polynomials and using a monomial basis would be like skipping the koornwinder polynomials and forming the shape functions using a basis of {1, x1, x2, x1x2, x1^2, x2^2, ...}.
"""

def koornwinder(pts_eval, p):
    """
    Driver for the koornwinder shape function generation
    "Tensor products of Jacobi polynomials"
    """
    raise NotImplementedError('Use other functions instead')

    dim = pts_eval.shape[1]
    if dim == 1:
        return koornwinder1d(pts_eval, p)
    elif dim == 2:
        return koornwinder2d(pts_eval, p)
    elif dim == 3:
        raise NotImplementedError('Dimension 3 not implemented yet')
    else:
        raise ValueError('Dim must be 1, 2, or 3')

def koornwinder1d(pts_eval, porder):
    """
    X: Coordinates of the points wherethe polynomials are to be evaluated (NPOINTS)
    P: Maximum order of the shape functions. That is all polynomials of degree up to P, NPOLY=P+1

    Returns:
    F: Vandermonde matrix (NPOINTS,NPOLY)
    FX: Vandermonde matrix for the derivative of the Koornwinder polynomials w.r.t. x (NPOINTS,NPOLY)
    """

    if pts_eval.shape[1] > 1:
        print('The input list of x coordinates must be 1 dimensional')
        raise ValueError

    # Perform a transformation from the master element (triangle) to a larger master triangle centered at the origin with side legnth 2 - the element on the right of diagram on pg 25 of Interpolation 16.930 notes
    pts_eval = 2*pts_eval-1

    npoly = porder+1
    # Empty vandermonde matrix for the shape functions
    V = np.zeros((pts_eval.shape[0], npoly))
    # Empty vandermonde matrix for the shape function derivatives
    Vx = np.zeros_like(V)

    for n in np.arange(npoly):
        # Normalization factor such that the function integrates to 1
        poly = jacobi(n, 0, 0) * (2*n+1)**0.5
        dpoly = poly.deriv(1)

        V[:, n] = np.squeeze(poly(pts_eval))
        Vx[:, n] = np.squeeze(dpoly(pts_eval)[:, None])

    Vx *= 2

    return V, Vx

def koornwinder2d(pts_eval, porder):
    """
        KOORNWINDER2D Vandermonde matrix for Koornwinder polynomials in
                    the master triangle[0, 0]-[1, 0]-[0, 1]
        [F, FX, FY] = KOORNWINDER(X, P)

            X:         Coordinates of the points wherethe polynomials
                        are to be evaluated(npoints, dim)
            PORDER:    Maximum order of the polynomials consider. That
        is all polynomials of complete degree up to p,
                        npoly = (PORDER+1)*(PORDER+2)/2
            F:         Vandermonde matrix(npoints, npoly)
            FX:        Vandermonde matrix for the derivative of the Koornwinder
                        polynomials w.r.t. x(npoints, npoly)
            FY:        Vandermonde matrix for the derivative of the Koornwinder
                        polynomials w.r.t. y(npoints, npoly)
    """

    if pts_eval.shape[1] != 2:
        raise ValueError('The input list of x coordinates must be 2 dimensional')

    poly_indices = pascalindex2d(porder)
    npoly = poly_indices.shape[0]   # Equivalent to (k+1)*(k+2)/2

    # Empty vandermonde matrices for the shape functions and their derivatives
    V = np.zeros((pts_eval.shape[0], npoly))
    Vx = np.zeros_like(V)
    Vy = np.zeros_like(V)

    # Points come in as coordinates on the master triangle with vertices (0, 0), (1, 0), and (0, 1). They must be transformed first to the standard triangle centered at the origin, and then to the standard square sitting on the origin between (-1, 1)x(-1, 1)
    # Perform a transformation from the master element (triangle) to a larger master triangle centered at the origin with side legnth 2 - the element on the right of diagram on pg 25 of Interpolation 16.930 notes
    xi = 2*pts_eval-1   # Operates on both x and y - stretching and shifting triangle so that midpoint of hypotenuse is on the origin - in "standard" configuration

    ind_1 = np.argwhere(xi[:, 1] > 0.99999999)  # Or, we could have hard set the singularity while performing the transformation (see koornwinder3d)
    # id_1 if the corner point exists (for example): [[9]], shape = (1, 1)
    # id_1 if the corner point does not exist: [], shape = (0, 1)
    if ind_1.shape[0]:   # There will not always be a point at the corner, for example if koornwinder is being used to evaluate at the GQpts
        xi[ind_1, 1] = 0.99999999

    # Perform coordinate transformation from triangle to square - taken from slide 25 of 16.930 notes.
    eta = np.copy(xi)
    eta[:, 0] = 2*(1+xi[:, 0])/(1-xi[:, 1]) - 1
    #eta2 = xi2 in the transformation

    # Saving points at which to evaluate the derivative before resettting the corner point to 1
    eta_displaced = np.copy(eta)
    # Moving the point back to the corner in the actual eta vector
    eta[ind_1] = 1      # More rigorous testing needed into whether this or eta[ind_1] = 1 makes a difference in the 2D case

    jac = np.zeros_like(xi)
    jac[:, 0] = 2/(1-xi[:, 1])
    jac[:, 1] = 2*(1+xi[:, 0])/(1-xi[:, 1])**2

    for idx, (m, n) in enumerate(poly_indices):
        Px = jacobi(m, 0, 0)
        Py = jacobi(n, 2*m+1, 0)*(poly([1, -1])/2)**m
        dPx = Px.deriv(1)
        dPy = Py.deriv(1)

        # IMPORTANT: The points used to calculate the shape function values include the corner (singularity) because the polynomials are defined there...
        Px_val = Px(eta[:, 0])
        Py_val = Py(eta[:, 1])

        # However, the jacobian is singular at the corner, so we can't evaluate that point directly in our integral. Instead, we look at the point that is slightly displaced from the corner and evaluate the derivative of the function at that displaced point.
        # Note that since we have split up the function into multiple parts to differentiate via the product rule, we still need to evaluate both the derivative and re-evaluate the shape function at the displaced point. The only difference between eta and eta_displaced is that the eta coord of the corner is shifted down slightly.
        Px_val_disp = Px(eta_displaced[:, 0])
        Py_val_disp = Py(eta_displaced[:, 1])
        dPx_val = dPx(eta_displaced[:, 0])
        dPy_val = dPy(eta_displaced[:, 1])

        # Normalization factor such that the function integrates to 1
        norm = (2*(2*m+1)*(m+n+1))**0.5

        V[:, idx] = Px_val * Py_val*norm
        Vx[:, idx] = dPx_val*Py_val_disp*jac[:, 0]*norm
        Vy[:, idx] = norm*(dPx_val*Py_val_disp*jac[:, 1] + Px_val_disp*dPy_val)

    Vx *= 2     # Derivatives are magnified because the domain is compressed in the transformation back to the master triangle
    Vy *= 2

    return V, Vx, Vy

def koornwinder3d(pts_eval, porder):
    """
        KOORNWINDER3D Vandermonde matrix for Koornwinder polynomials in
                    the master tetrahedron [0, 0]-[1, 0]-[0, 1]-[1, 1]
        [F, FX, FY, FZ] = KOORNWINDER(X, P)

            X:         Coordinates of the points wherethe polynomials
                        are to be evaluated(npoints, dim)
            PORDER:    Maximum order of the polynomials consider. That
        is all polynomials of complete degree up to p,
                        npoly = (PORDER+1)*(PORDER+2)/2
            F:         Vandermonde matrix(npoints, npoly)
            FX:        Vandermonde matrix for the derivative of the Koornwinder
                        polynomials w.r.t. x(npoints, npoly)
            FY:        Vandermonde matrix for the derivative of the Koornwinder
                        polynomials w.r.t. y(npoints, npoly)
            FZ:        Vandermonde matrix for the derivative of the Koornwinder
                        polynomials w.r.t. z(npoints, npoly)
    """

    if pts_eval.shape[1] != 3:
        print('The input list of x coordinates must be 3 dimensional')
        raise ValueError

    poly_indices = pascalindex3d(porder)
    npoly = poly_indices.shape[0]   # Equivalent to (k+1)*(k+2)*(k+3)/6
    numpts = pts_eval.shape[0]

    # Empty vandermonde matrices for the shape functions and their derivatives
    V = np.zeros((numpts, npoly))
    Vx = np.zeros_like(V)
    Vy = np.zeros_like(V)
    Vz = np.zeros_like(V)

    # Points come in as coordinates on the master tetrahedron with vertices (0, 0), (1, 0), (0, 1), and (1, 1). They must be transformed first to the standard triangle centered at the origin, and then to the standard cube sitting on the origin between (-1, 1)^3
    # Perform a transformation from the master element (triangle) to a larger master triangle centered at the origin with side legnth 2 - the element on the right of diagram on pg 25 of Interpolation 16.930 notes
    xi = 2*pts_eval-1   # Operates on both x and y - stretching and shifting triangle so that midpoint of hypotenuse is on the origin - in "standard" configuration

    """
    DANGER: Known issue with calculation of gradient of koornwinder polynomials at the corners of the triangle/tet due to singularity of transformation. Four cases:
    KP = koornwinder polynomials
    1) KPs sampled at corners
    2) KP derivatives sampled at corners
    3) KPs sampled at GQ pts (exclusively interior)
    4) KP derivatives sampled at GQ pts (exclusively interior)

    Case 1 is used for interpolating onto a higher mesh
    Case 2 is used for computing normal vectors and the gradients of a scalar field                 THIS IS THE PROBLEM - initially when this code was written, this use case wasn't envisioned and wasn't a problem because of the known issue with sampling the derivatives at the corners
    Case 3 is used during normal construction of the elemental matrices during matrix assembly
    Case 4 is used during normal construction of the elemental matrices during matrix assembly

    """

    # # FIX
    # ind_1 = np.argwhere(xi[:, 1] > 0.99999999)
    # # id_1 if the corner point exists: [[9]], shape = (1, 1)
    # # id_1 if the corner point does not exist: [], shape = (0, 1)
    # if ind_1.shape[0]:   # There will not always be a point at the corner, for example if koornwinder is being used to evaluate at the GQpts
    #     xi[ind_1, 1] = 0.99999999

    # # FIX
    # # Perform coordinate transformation from triangle to square - taken from slide 25 of 16.930 notes. Note that eta2 = xi2 in the transformation.
    # eta = np.copy(xi)
    # eta[:, 0] = -2*(1+xi[:,0])/(xi[:,1]+xi[:,2]) - 1
    # eta[:, 1] = 2*(1+xi[:, 1])/(1-xi[:, 2]) - 1
    # # eta3 is equal to xi3 in the transformation

    # # Saving points at which to evaluate the derivative before resettting the corner point to 1
    # eta_displaced = np.copy(eta)
    # # Moving the point back to the corner in the actual eta vector
    # eta[ind_1] = 1

    eta = np.zeros_like(xi)
    eta_displaced = np.zeros_like(xi)

    # "Manual transform" - the points are transformed in a case-by-case basis depending on whether or not they lie at corners
    xi_displaced = np.copy(xi)
    eps = 1e-8
    for i, pt in enumerate(xi):
        if np.isclose(abs((pt[1] + pt[2])), 0, atol=eps):   # This disagrees with Cuong's implementation because it checks for the sum being close to 0, instead of checking that it is exactly equal to 0, which doesn't account for floating point roundoff.
            eta[i,0] = -1
            xi_displaced[i,2] = -pt[1] - eps
        else:
            # Normal eta1 transform
            eta[i, 0] = -2*(1+pt[0])/(pt[1]+pt[2]) - 1
        if np.isclose(pt[2], 1, atol=eps):
            eta[i, 1] = -1
            xi_displaced[i,2] = 1-eps
        else:
            # Normal eta 2 transform
            eta[i, 1] = 2*(1+pt[1])/(1-pt[2]) - 1

        eta[i:,2] = pt[2]

    eta_displaced[:, 0] = -2*(1+xi_displaced[:,0])/(xi_displaced[:,1]+xi_displaced[:, 2]) - 1
    eta_displaced[:, 1] = 2*(1+xi_displaced[:, 1])/(1-xi_displaced[:, 2]) - 1
    eta_displaced[:, 2] = xi_displaced[:,2]

    jac = np.zeros((numpts, 3,3))
    jac[:, 0, 0] = -2/(xi_displaced[:, 1] + xi_displaced[:, 2])
    jac[:, 1, 0] = 2*(1+xi_displaced[:, 0])/(xi_displaced[:, 1] + xi_displaced[:, 2])**2
    jac[:, 2, 0] = 2*(1+xi_displaced[:, 0])/(xi_displaced[:, 1] + xi_displaced[:, 2])**2  # same as above
    jac[:, 0, 1] = 0
    jac[:, 1, 1] = 2/(1-xi_displaced[:, 2])
    jac[:, 2, 1] = 2*(1+xi_displaced[:, 1])/(1-xi_displaced[:, 2])**2
    jac[:, 0, 2] = 0
    jac[:, 1, 2] = 0
    jac[:, 2, 2] = 1

    for idx, (m, n, l) in enumerate(poly_indices):
        Px = jacobi(m, 0, 0)
        Py = jacobi(n, 2*m+1, 0)*(poly([1, -1])/2)**m
        Pz = jacobi(l, 2*m+2*n+2, 0)*(poly([1, -1])/2)**(m+n)
        dPx = Px.deriv(1)
        dPy = Py.deriv(1)
        dPz = Pz.deriv(1)

        # IMPORTANT: The points used to calculate the shape function values include the corner (singularity) because the polynomials are defined there...
        Px_val = Px(eta[:, 0])
        Py_val = Py(eta[:, 1])
        Pz_val = Pz(eta[:, 2])

        # Normalization factor such that the function integrates to 1
        norm = (2*(2*m+1) *2*(m+n+1) * (m+n+l+2))**0.5
        # print(norm)
        # Why not add in a factor of 2^(4i+2j+6)?

        # However, the jacobian is singular at the corner, so we can't evaluate that point directly in our integral. Instead, we look at the point that is slightly displaced from the corner and evaluate the derivative of the function at that displaced point.
        # Note that since we have split up the function into multiple parts to differentiate via the product rule, we still need to evaluate both the derivative and re-evaluate the shape function at the displaced point. The only difference between eta and eta_displaced is that the eta coord of the corner is shifted down slightly.
        Px_val_disp = Px(eta_displaced[:, 0])
        Py_val_disp = Py(eta_displaced[:, 1])
        Pz_val_disp = Pz(eta_displaced[:, 2])

        dPx_val = dPx(eta_displaced[:, 0])
        dPy_val = dPy(eta_displaced[:, 1])
        dPz_val = dPz(eta_displaced[:, 2])

        df_deta1 = dPx_val*Py_val_disp*Pz_val_disp*norm
        df_deta2 = Px_val_disp*dPy_val*Pz_val_disp*norm
        df_deta3 = Px_val_disp*Py_val_disp*dPz_val*norm

        for pt in np.arange(numpts):
            [df_dxi1, df_dxi2, df_dxi3] = jac[pt,:,:] @ np.array([df_deta1[pt], df_deta2[pt], df_deta3[pt]])

            Vx[pt, idx] = df_dxi1
            Vy[pt, idx] = df_dxi2
            Vz[pt, idx] = df_dxi3

        V[:, idx] = norm*Px_val*Py_val*Pz_val

    Vx *= 2
    Vy *= 2
    Vz *= 2

    return V, Vx, Vy, Vz

if __name__ == '__main__':
    pass
    # np.set_printoptions(suppress=True, linewidth=np.inf)
    # np.set_printoptions(linewidth=np.inf, precision=10)
    # import sys
    # sys.path.insert(0, '../mesh')
    # sys.path.insert(0, '../util')
    # from mkmshlocal import mkmshlocal
    # from import_util import load_mat

    # pl2d, _ = mkmshlocal(3)
    # pl2d = pl2d[:,1:]
    # loc1d_idx = np.squeeze(np.argwhere(pl2d[:, -1] == 0))
    # pl1d = pl2d[loc1d_idx, :-1]

    # Testing 1D
    # V, Vx = koornwinder(pl1d[:,1:], 3)
    # print(V)
    # print(Vx)

    # Testing 2D
    # V, Vx, Vy = koornwinder(pl2d, 3)
    # print(V)
    # print()
    # print(Vx)
    # print()
    # print(Vy)

    # V_mat = load_mat('V')
    # Vx_mat = load_mat('Vx')
    # Vy_mat = load_mat('Vy')

    # print(np.allclose(V, V_mat))
    # print(np.allclose(Vx, Vx_mat))
    # print(np.allclose(Vy, Vy_mat))

    # Testing 3D

    from scipy.io import loadmat
    gq_mat = loadmat('porder3ndim3GQ.mat')
    dg_mat = loadmat('porder3ndim3DG.mat')
    e_mat = loadmat('e.mat')['e']
    xc_mat = loadmat('xc.mat')['xc']
    xi_mat = loadmat('xi.mat')['x']

    from gaussquad2d import gaussquad1d, gaussquad2d, gaussquad3d
    from masternodes import masternodes
    np.set_printoptions(suppress=True, linewidth=np.inf, precision=12)

    porder=3
    ploc, _, _, _, _, _, _ = masternodes(porder, 3)
    gpts, _ = gaussquad3d(6)

    f, fx, fy, fz = koornwinder3d(ploc, porder)

    print(dg_mat['fx']-fx)
    # print()
    # print(np.allclose(dg_mat['fx'], fx, rtol=0, atol=1e-13))
    # print(ploc)
    # print(fx)

reason: import numpy,from numpy,from scipy


text:
import sys
import itertools

import gym
import matplotlib
import numpy as np

import sklearn.pipeline
import sklearn.preprocessing
from sklearn.linear_model import SGDRegressor
from sklearn.kernel_approximation import RBFSampler

from core import plotting

class Estimator():
    def __init__(self, env, phi, p):
        self._env = env
        self._phi = phi
        self._p = p
        self._n_action = env.action_space.n

        self.models = []
        for _ in range(self._n_action):
            model = SGDRegressor(learning_rate="constant")
            s0 = env.reset()
            a0 = 0 #arbitrary

            l = int(self._p*a0/self._n_action)
            r = int(self._p*(a0+1)/self._n_action)
            phi_sa = self._phi(s0, a0)[0, l:r]
            model.partial_fit([phi_sa], [0])
            self.models.append(model)

    def predict(self, s, a=None):
        if a is None:
            Q = []
            for m, a in zip(self.models, range(self._env.action_space.n)):
                l = int(self._p*a/self._n_action)
                r = int(self._p*(a+1)/self._n_action)
                features = self._phi(s, a)[0, l:r]
                Q.append(m.predict([features])[0])
            return np.array(Q)
        else:
            l = int(self._p*a/self._n_action)
            r = int(self._p*(a+1)/self._n_action)
            features = self._phi(s, a)[0, l:r]
            return self.models[a].predict([features])[0]

    def update(self, s, a, y):
        """
        Updates the estimator parameters for a given state and action towards
        the target y.
        """
        l = int(self._p*a/self._n_action)
        r = int(self._p*(a+1)/self._n_action)
        features = self._phi(s, a)[0, l:r]
        self.models[a].partial_fit([features], [y])

class EpsilonGreedyPolicy(object):
    def __init__(self, estimator, eps, nA):
        self._estimator = estimator
        self._eps = eps
        self._nA = nA

    def _action_probs(self, s):
        probs = np.ones(self._nA, dtype=float) * self._eps / self._nA
        q_values = self._estimator.predict(s)
        best_action = np.argmax(q_values)
        probs[best_action] += (1.0 - self._eps)
        return probs

    def choose_action(self, s):
        p = self._action_probs(s)
        a = np.random.choice(np.arange(len(p)), p=p)
        return a

    @property
    def coef_(self):
        """TODO: Docstring for function.

        Parameters
        ----------

        Returns
        -------
        coef : array
            dim_action by dim_feature array

        """
        return np.array([np.hstack((m.coef_, m.intercept_)) for m in self._estimator.models])

class LinearQ(object):
    """Docstring for LinearQ3. """

    def __init__(self, env, phi, p, action_list, n_episode, epsilon, epsilon_decay, gamma):
        """TODO: to be defined1.

        Parameters
        ----------
        env : TODO
        phi : array
        action_list : list
        n_episode : TODO
        epsilon : TODO
        epsilon : TODO
        plotting : TODO
        gamma: float
            discount factor

        """
        self._env = env
        self._estimator = Estimator(env=env, phi=phi, p=p)
        self._n_episode = n_episode
        self._epsilon = epsilon
        self._epsilon_decay = epsilon_decay
        self._gamma = gamma

    def solve(self, reward_fn=None):
        """TODO: Docstring for solve.

        Parameters
        ----------
        reward_fn : function
            reward function for IRL

        Returns
        -------
        TODO

        """
        # Keeps track of useful statistics
        stats = plotting.EpisodeStats(
            episode_lengths=np.zeros(self._n_episode),
            episode_rewards=np.zeros(self._n_episode))

        for i_episode in range(self._n_episode):

            pi = EpsilonGreedyPolicy(estimator=self._estimator,
                                     eps=self._epsilon * self._epsilon_decay**i_episode,
                                     nA=self._env.action_space.n)

            last_reward = stats.episode_rewards[i_episode - 1]
            sys.stdout.flush()

            state = self._env.reset()

            for t in itertools.count():

                action = pi.choose_action(state)

                next_state, reward, done, _ = self._env.step(action)

                stats.episode_rewards[i_episode] += reward
                stats.episode_lengths[i_episode] = t

                q_values_next = self._estimator.predict(next_state)

                if reward_fn is not None:
                    reward_irl = reward_fn(state, action)
                    td_target = reward_irl + self._gamma * np.max(q_values_next)
                else:
                    td_target = reward + self._gamma * np.max(q_values_next)

                self._estimator.update(state, action, td_target)

                print("\rStep {} @ Episode {}/{} ({})".format(t, i_episode + 1, self._n_episode, last_reward), end="")

                if done:
                    break

                state = next_state

        return pi, stats

if __name__ == "__main__":
    estimator = Estimator()

    stats = q_learning(env, estimator, 100, epsilon=0.0)

    plotting.plot_cost_to_go_mountain_car(env, estimator)
    plotting.plot_episode_stats(stats, smoothing_window=25)

max_stars_repo_name: dtak/batch-apprenticeship-learning

max_issues_repo_name: dtak/batch-apprenticeship-learning

max_forks_repo_name: dtak/batch-apprenticeship-learning


text:
# -*- coding: utf-8 -*-
"""Functions to make 3D plots with M/EEG data."""

# Authors: Alexandre Gramfort <alexandre.gramfort@inria.fr>
#          Denis Engemann <denis.engemann@gmail.com>
#          Martin Luessi <mluessi@nmr.mgh.harvard.edu>
#          Eric Larson <larson.eric.d@gmail.com>
#          Mainak Jas <mainak@neuro.hut.fi>
#          Mark Wronkiewicz <wronk.mark@gmail.com>
#
# License: Simplified BSD

import numpy as np
from mne.io.constants import FIFF
from mne.io.pick import _picks_to_idx
from mne.utils import _validate_type, fill_doc, verbose

verbose_dec = verbose
FIDUCIAL_ORDER = (FIFF.FIFFV_POINT_LPA, FIFF.FIFFV_POINT_NASION,
                  FIFF.FIFFV_POINT_RPA)

@fill_doc
def plot_sensors_connectivity(info, con, picks=None,
                              cbar_label='Connectivity'):
    """Visualize the sensor connectivity in 3D.

    Parameters
    ----------
    info : dict | None
        The measurement info.
    con : array, shape (n_channels, n_channels) | Connectivity
        The computed connectivity measure(s).
    %(picks_good_data)s
        Indices of selected channels.
    cbar_label : str
        Label for the colorbar.

    Returns
    -------
    fig : instance of Renderer
        The 3D figure.
    """
    _validate_type(info, "info")

    from mne.viz.backends.renderer import _get_renderer
    from mne_connectivity.base import BaseConnectivity

    if isinstance(con, BaseConnectivity):
        con = con.get_data()

    renderer = _get_renderer(size=(600, 600), bgcolor=(0.5, 0.5, 0.5))

    picks = _picks_to_idx(info, picks)
    if len(picks) != len(con):
        raise ValueError('The number of channels picked (%s) does not '
                         'correspond to the size of the connectivity data '
                         '(%s)' % (len(picks), len(con)))

    # Plot the sensor locations
    sens_loc = [info['chs'][k]['loc'][:3] for k in picks]
    sens_loc = np.array(sens_loc)

    renderer.sphere(np.c_[sens_loc[:, 0], sens_loc[:, 1], sens_loc[:, 2]],
                    color=(1, 1, 1), opacity=1, scale=0.005)

    # Get the strongest connections
    n_con = 20  # show up to 20 connections
    min_dist = 0.05  # exclude sensors that are less than 5cm apart
    threshold = np.sort(con, axis=None)[-n_con]
    ii, jj = np.where(con >= threshold)

    # Remove close connections
    con_nodes = list()
    con_val = list()
    for i, j in zip(ii, jj):
        if np.linalg.norm(sens_loc[i] - sens_loc[j]) > min_dist:
            con_nodes.append((i, j))
            con_val.append(con[i, j])

    con_val = np.array(con_val)

    # Show the connections as tubes between sensors
    vmax = np.max(con_val)
    vmin = np.min(con_val)
    for val, nodes in zip(con_val, con_nodes):
        x1, y1, z1 = sens_loc[nodes[0]]
        x2, y2, z2 = sens_loc[nodes[1]]
        tube = renderer.tube(origin=np.c_[x1, y1, z1],
                             destination=np.c_[x2, y2, z2],
                             scalars=np.c_[val, val],
                             vmin=vmin, vmax=vmax,
                             reverse_lut=True)

    renderer.scalarbar(source=tube, title=cbar_label)

    # Add the sensor names for the connections shown
    nodes_shown = list(set([n[0] for n in con_nodes] +
                           [n[1] for n in con_nodes]))

    for node in nodes_shown:
        x, y, z = sens_loc[node]
        renderer.text3d(x, y, z, text=info['ch_names'][picks[node]],
                        scale=0.005,
                        color=(0, 0, 0))

    renderer.set_camera(azimuth=-88.7, elevation=40.8,
                        distance=0.76,
                        focalpoint=np.array([-3.9e-4, -8.5e-3, -1e-2]))
    renderer.show()
    return renderer.scene()

max_stars_repo_name: BIAPT/mne-connectivity

max_issues_repo_name: BIAPT/mne-connectivity

max_forks_repo_name: BIAPT/mne-connectivity


text:
#%%
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import phd.thermo
import phd.viz
colors, palette = phd.viz.phd_style()
INDUCTION = False
MUTS = False
GROWTH = False

# Load the compiled datasets
data = pd.read_csv('../../src/data/other/Garcia2011_Brewster2014_RazoMejia2018_Chure2019.csv',
                   comment='#')
growth_data = pd.read_csv('../../src/data/ch4_growth/analyzed_foldchange.csv', comment="#")
entropy_samples = pd.read_csv('../../src/data/ch4_growth/pooled_entropic_parameter_samples.csv', comment='#')

# Summarize the growth data.
growth_data = growth_data.groupby(['atc_ngml', 'date', 'run_number',
                                    'carbon', 'temp'])[
                                    ['repressors', 'fold_change']].mean().reset_index()
growth_data = growth_data.groupby(['atc_ngml', 'carbon', 'temp']).agg(('mean', 'sem')).reset_index()
growth_data

# Compute the bohr parameter
dfs = []
for g, d in growth_data.groupby(['carbon', 'temp']):
    if g[1] != 37:
        _samples = entropy_samples[(entropy_samples['temp']==g[1])]
        epRA = _samples[_samples['parameter']=='epRA_star']['value'].mean()
        epAI = _samples[_samples['parameter']=='epAI_star']['value'].mean()
    else:
        epRA = -13.9
        epAI = 4.5

    bohr = phd.thermo.SimpleRepression(d['repressors']['mean'] * 1.16, epRA, ka=139,
                           ki=0.53, ep_ai=epAI, effector_conc=0).bohr_parameter()
    d['bohr_parameter'] = bohr
    dfs.append(d)
growth = pd.concat(dfs, sort=False)

# Set up the figure canvas
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
phd.viz.despine(ax)
bohr_range = np.linspace(-12, 12, 200)
scaling_fn = (1 + np.exp(-bohr_range))**-1
_ = ax.plot(bohr_range, scaling_fn, 'k-', label='scaling function')
ax.set_xlabel('free energy [$k_BT$]')
ax.set_ylabel('fold-change')
ax.set_ylim([-0.05, 1.25])
ax.set_xlim([-12, 12])

induction = data[data['author']=='razo-mejia']

rep_colors = {22: colors['red'],
              60: colors['brown'],
              124: colors['green'],
              260: colors['orange'],
              1220: colors['purple'],
              1740: colors['blue']}
op_glyphs = {'O1':'^', 'O2':'o', 'O3':'v'}
class_glyphs = {'DNA': 'h', 'IND':'p', 'DBL':'X'}
mut_colors = {'Y20I':colors['red'], 'Q21A':colors['purple'], 'Q21M':colors['light_purple'],
              'F164T':colors['green'], 'Q294R':colors['brown'], 'Q294V':colors['orange'],
              'Q294K':colors['light_grey'], 'Y20I-F164T':colors['black'],
              'Y20I-Q294V':colors['dark_brown'], 'Y20I-Q294K':colors['dark_red'],
              'Q21A-F164T':colors['dark_blue'], 'Q21A-Q294V':colors['dark_purple'],
              'Q21A-Q294K':colors['dark_green'], 'Q21M-Q294V':colors['orange'],
              'Q21M-F164T':colors['pale_green'], 'Q21M-Q294K':colors['light_red'],
              'Q21M-Q294V':colors['light_purple']}

temp_colors = {32:colors['blue'], 42:colors['red']}
carbon_colors = {'glycerol':colors['green'], 'acetate':colors['brown']}

name = 'scaling_fn'
if INDUCTION:
    name += '_INDUCTION'
    for o, g in op_glyphs.items():
        _ = ax.plot([], [], f'{g}', color=colors['black'],
                    markeredgecolor=colors['grey'], markeredgewidth=0.5,
                    label=f'operator {o}', ms=4)
    for r, c in rep_colors.items():
        _ = ax.plot([], [], '-', color=c, label=f'{int(r)} repressors',
        ms=4)
    for g, d in induction.groupby(['repressors', 'operator']):
        marker = op_glyphs[g[1]]
        _color = rep_colors[g[0]]
        _ = ax.errorbar(d['bohr_parameter'], d['mean'], d['sem'], color=_color,
                        markeredgecolor=colors['grey'], fmt=marker, linestyle='none',
                        ms=4, markeredgewidth=0.5, label='__nolegend__')

DNA_muts = ['Q21M', 'Y20I', 'Q21A']
IND_muts =['F64T', 'Q294K', 'Q294V']
if MUTS:
    name += '_MUTS'
    for c, g in class_glyphs.items():
        if c == 'DNA':
            label = 'DNA binding mutant'
            _color = colors['purple']
        elif c == 'IND':
            label = 'inducer binding mutant'
            _color = colors['green']
        else:
             label = 'double mutant'
             _color = colors['blue']

        _ = ax.plot([], [], f'{g}', color=_color, markeredgecolor=colors['grey'],
                    markeredgewidth=0.5, label=label, ms=4)
    for g, d in data[(data['author']=='chure') & (data['mutant']!='wt')].groupby(['mutant']):
        if (g in DNA_muts) & ('-' not in g):
            marker = class_glyphs['DNA']
            _color = colors['purple']
        if (g in IND_muts) & ('-' not in g):
            marker = class_glyphs['IND']
            _color = colors['green']
        elif '-' in g:
            marker = class_glyphs['DBL']
            _color = colors['blue']

        _ = ax.errorbar(d['bohr_parameter'], d['mean'], d['sem'], fmt=marker, ms=4,
                    color=_color, markeredgewidth=0.5, markeredgecolor=colors['grey'],
                    linestyle='none', label='__nolegend__')

if GROWTH:
    name += '_GROWTH'
    _ = ax.plot([], [], '>', color=colors['red'], label='grown at 42°C', ms=4,
                markeredgewidth=0.5, markeredgecolor=colors['grey'])
    _ = ax.plot([], [], '>', color=colors['blue'], label='grown at 32°C', ms=4,
                markeredgewidth=0.5, markeredgecolor=colors['grey'])
    _ = ax.plot([], [], '<', color=colors['green'], label='grown in glycerol', ms=4,
                markeredgewidth=0.5, markeredgecolor=colors['grey'])
    _ = ax.plot([], [], '<', color=colors['brown'], label='grown in acetate', ms=4,
                markeredgewidth=0.5, markeredgecolor=colors['grey'])
    for g, d in growth.groupby(['temp', 'carbon']):
        if (g[0] != 37) | (g[1]!='glucose'):
            if g[0] != 37:
                marker = '>'
                _color = temp_colors[g[0]]
            else:
                marker = '<'
                _color = carbon_colors[g[1]]
            _ = ax.errorbar(d['bohr_parameter'], d['fold_change']['mean'],
                        d['fold_change']['sem'], fmt=marker, ms=4,
                        color=_color, linestyle='none', markeredgewidth=0.5,
                        markeredgecolor=colors['grey'], label='__nolegend__',
                        zorder=1000)
ax.legend(loc='upper left')

plt.savefig(f'../figs/collapse_{name}.pdf')

# %%

# %%

max_stars_repo_name: gchure/phd

max_issues_repo_name: gchure/phd

max_forks_repo_name: gchure/phd


text:
import matplotlib
matplotlib.use('Agg')

from pySDC.helpers.stats_helper import filter_stats

import numpy as np

from matplotlib import rc
import matplotlib.pyplot as plt

# noinspection PyShadowingBuiltins
def show_residual_across_simulation(stats, fname='residuals.png'):
    """
    Helper routine to visualize the residuals across the simulation (one block of PFASST)

    Args:
        stats (dict): statistics object from a PFASST run
        fname (str): filename
    """

    # get residuals of the run
    extract_stats = filter_stats(stats, type='residual_post_iteration')

    # find boundaries for x-,y- and c-axis as well as arrays
    maxprocs = 0
    maxiter = 0
    minres = 0
    maxres = -99
    for k, v in extract_stats.items():
        maxprocs = max(maxprocs, getattr(k, 'process'))
        maxiter = max(maxiter, getattr(k, 'iter'))
        minres = min(minres, np.log10(v))
        maxres = max(maxres, np.log10(v))

    # grep residuals and put into array
    residual = np.zeros((maxiter, maxprocs + 1))
    residual[:] = -99
    for k, v in extract_stats.items():
        step = getattr(k, 'process')
        iter = getattr(k, 'iter')
        if iter is not -1:
            residual[iter - 1, step] = np.log10(v)

    # Set up latex stuff and fonts
    rc('font', **{"sans-serif": ["Arial"], "size": 30})
    rc('legend', fontsize='small')
    rc('xtick', labelsize='small')
    rc('ytick', labelsize='small')

    # create plot and save
    fig, ax = plt.subplots(figsize=(15, 10))

    cmap = plt.get_cmap('Reds')
    plt.pcolor(residual.T, cmap=cmap, vmin=minres, vmax=maxres)

    cax = plt.colorbar()
    cax.set_label('log10(residual)')

    ax.set_xlabel('iteration')
    ax.set_ylabel('process')

    ax.set_xticks(np.arange(maxiter) + 0.5, minor=False)
    ax.set_yticks(np.arange(maxprocs + 1) + 0.5, minor=False)
    ax.set_xticklabels(np.arange(maxiter) + 1, minor=False)
    ax.set_yticklabels(np.arange(maxprocs + 1), minor=False)

    plt.savefig(fname, rasterized=True, transparent=True, bbox_inches='tight')

max_stars_repo_name: janEbert/pySDC

max_issues_repo_name: janEbert/pySDC

max_forks_repo_name: janEbert/pySDC


text:
#!/usr/bin/env python
"""
Abstract Base Class for audio feature extraction. Defines an interface for
audio feature extraction algorithms.

Inheriting classes must override :py:meth:`~spiegelib.features.FeaturesBase.get_features` and implement feature
extraction in that function, which accepts an :ref:`AudioBuffer <audio_buffer>` and
must return a `np.ndarray` containing results.

Feature extraction is run in inheriting classes through the __call__ member function.

Example of FFT, which inherits from FeaturesBase::

    audio = AudioBuffer('./some_audio.wav')
    fft = FFT()

    # Now feature extraction is run by treating
    # the FFT instance, fft, like a function
    spectrum = fft(audio)
"""

from abc import ABC, abstractmethod
import functools
import os.path
import numpy as np
import joblib
from tqdm import trange

from spiegelib.features import StandardScaler

class FeaturesBase(ABC):
    """
    Args:
        sample_rate (int, optional): audio sample rate, defaults to 44100
        time_major (bool, optional): for audio feature extraction that uses time slices,
            used to indicate the orientation of features and time slices in the resulting
            matrix. True indicates that results should be returned with an orientation
            like (time_slices, features). False (default) refers to an orientation
            of (features, time_slices).
        scale (bool, optional): whether to scale the results of feature extraction. The scaler
            must be set before scaling called be applied. See :py:meth:`~spiegelib.features.FeaturesBase.fit_scaler` or
            :py:meth:`~spiegelib.features.FeaturesBase.set_scaler`.
        scale_axis (int, tuple, None, optional): indicates the axis to use for fitting scalers and
            applying data scaling. A value of None flattens the dataset and calculates
            scaling parameters on that. A value of an int or a tuple indicates the axis
            or axes to use. Defaults to (0,) which causes scaling variables to be calculated
            on each feature independently.
    """

    def __init__(self, sample_rate=44100, time_major=False, scale=False, scale_axis=(0,)):
        """
        Constructor
        """

        super().__init__()

        self.sample_rate = sample_rate
        self.time_major = time_major

        self.should_scale = scale
        self.scaler = None
        self.scale_axis = scale_axis
        self.ScalerClass = StandardScaler

        self.input_modifiers = []
        self.prescale_modifiers = []
        self.output_modifiers = []

        # Update this in inheriting classes if you need
        self.dtype = np.float32

    def __call__(self, audio, scale=None):
        """
        Run this feature extraction pipeline.

        Applies functions in this order: input modifiers > feature extraction >
        prescale modifiers > data scaling > output modifiers

        Args:
            audio (:ref:`AudioBuffer <audio_buffer>`): Audio to extract features from
            scale (bool, optional): If set, will override scale attribute set
                during construction.

        Returns:
            np.ndarray: results from audio feature extraction with modifiers and scaling.
        """

        # Input data modification
        for modifer in self.input_modifiers:
            audio = modifer(audio)

        # Run feature extraction
        features = self.get_features(audio)

        # Apply any prescaling data modification
        for modifier in self.prescale_modifiers:
            features = modifier(features)

        # Normalize features
        should_scale = scale if scale != None else self.should_scale
        if should_scale:
            assert self.has_scaler(), "Scaler must be set first."
            features = self.scale(features)

        # Apply any output data modification
        for modifier in self.output_modifiers:
            features = modifier(features)

        return features

    @abstractmethod
    def get_features(self, audio):
        """
        Must be implemented. Run audio feature extraction on audio provided as parameter.
        Must check the `time_major` attribute and return data in the correct
        orientation.

        Args:
            audio (:ref:`AudioBuffer <audio_buffer>`): Audio to extract features from

        Returns:
            np.ndarray: Results of audio feature extraction
        """
        pass

    def add_modifier(self, modifier, type):
        """
        Add a data modifier to the feature extraction pipeline.

        Input modifiers are applied to raw `AudioBuffers <audio_buffer>` prior
        to feature extraction. Prescale modifiers are applied to results of audio
        feature extraction and before data scaling (if applicable). Output modifiers
        are applied after data scaling.

        Args:
            modifier (lambda): data modifier function. Should accept an np.array,
                apply some modification to that, and return a np.array
            type (str): Where to add function into pipeline. Must be one of
                ('input', 'prescale', or 'output')
        """

        if type == 'input':
            self.input_modifiers.append(modifier)
        elif type == 'prescale':
            self.prescale_modifiers.append(modifier)
        elif type == 'output':
            self.output_modifiers.append(modifier)
        else:
            raise ValueError('type must be one of ("input", "prescale", or '
                             '"output"), received: %s' % type)

    def fit_scaler(self, data, transform=True):
        """
        Fit scaler to dataset for future transforms.

        Args:
            data (np.ndarray): data to train (fit) scaler on
            transform (bool, optional): if True will also apply scaling
                to the data

        Returns:
            np.ndarray, None: Scaled data if transform parameter is True, otherwise\
                returns None.
        """

        # Fit a new scaler to dataset
        self.scaler = self.ScalerClass()
        self.scaler.fit(data, self.scale_axis)

        # Scale the input dataset
        if transform:
            return self.scaler.transform(data)
        else:
            return None

    def has_scaler(self):
        """
        Returns:
            bool: whether or not the scaler has been set
        """
        return self.scaler != None

    def set_scaler(self, scaler):
        """
        Set a scaler for a dimension, this will be used to normalize that dimension

        Args:
            scaler (object): a Scaler object to use to scale data. Must be a
                :ref:`DataScalerBase <data_scaler_base>` type to use the :py:meth:`~spiegelib.features.FeaturesBase.fit_scaler`
                method. Otherwise, other scalers like sklearn scalers could potentially
                be fit and then passed in here -- just needs to have a transform method
                to apply scaling. (note, using other objects outside of :ref:`DataScalerBase <data_scaler_base>` type
                has not been tested!)
        """

        self.scaler = scaler

    def scale(self, data):
        """
        Scale features using pre-trained scaler

        Args:
            data (np.ndarray): data to be scaled

        Returns:
            np.ndarray: scaled data
        """

        assert self.has_scaler(), "Scaler must be set first"
        return self.scaler.transform(data)

    def load_scaler(self, location):
        """
        Load trained scaler from a pickled file.

        Args:
            location (str): Location of pickled scaler object
        """

        self.scaler = joblib.load(location)

    def save_scaler(self, location):
        """
        Save the trained scaler for these features as a pickle for later use

        Args:
            location (str): Location to save pickled scaler object
        """

        joblib.dump(self.scaler, location)

max_stars_repo_name: spiegel-lib/spiegel

max_issues_repo_name: spiegel-lib/spiegel

max_forks_repo_name: spiegel-lib/spiegel


text:
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 11 10:06:15 2019

@author: Nathan
"""

import numpy as np
import random
from random import sample, choices, randrange
from heuristic import greedy
from itertools import accumulate, combinations
from bisect import bisect_left

def aco(locs, size=52, factor=.68, decay=.56):
    """
    Arguments:
        locs (atlas): atlas type object
        size (int): how many ants to use
        factor (float): pheromone factor for the ants must be > 0
        decay (float): decay factor for phermomone between 0. and 1.
    Yields:
        list: best path for this iteration
    """

    # residual pheromone and current iteration pheromone
    pheromone = np.ones_like(locs.dist)
    delta = np.empty_like(pheromone)

    while True:
        # best path and reset current iteration pheromone
        min_dist, min_path = float('inf'), None
        delta.fill(0.0)

        for _ in range(size):
            # current loc, path traveled, locs to vist yet
            cur, path, todo = 0, [0], list(range(1,len(locs)))

            # find cycle based on probability
            for _ in range(len(todo)):
                attraction =[pheromone[cur][i]/locs.dist[cur][i] for i in todo]
                total = sum(attraction)
                attraction = list(accumulate(x / total for x in attraction))
                cur = todo.pop(bisect_left(attraction, random.random()))
                path.append(cur)

            # update best path
            dist = locs.distance(path)
            if dist < min_dist:
                min_dist, min_path = dist, path

            # update delta pheromones
            for a, b in zip(path, path[1:] + [0]):
                delta[a][b] += factor / dist

        # update pheromones for decay, and symmetric delta pheromones
        pheromone = pheromone * decay + delta + delta.transpose()

        yield min_path

def pso(locs, size=1000, a=.85, b=.6):
    """
    Arguments:
        locs (atlas): atlas type object
        size (int): how many particles to use
        a (float): what amount of personal best swaps to use
        b (float): what amount of global best swaps to use
    Yields:
        list: path found using pso algorithm after each epoch
    """
    length = len(locs)
    particles = [sample(range(length), length) for _ in range(size)]
    pbest, gbest = particles.copy(), None

    while True:
        # find gbest
        gbest = min(particles, key=locs.distance)

        for num, pos in enumerate(particles):
            # copies of pbest and gbest
            alpha, beta, deltas = pbest[num][:], gbest[:], list()
            # chance to flop pos from diff with personal best
            for i in range(length):
                if pos[i] != alpha[i]:
                    j = alpha.index(pos[i])
                    alpha[i], alpha[j] = alpha[j], alpha[i]
                    if random.random() < a:
                        deltas.append((i, j))
            # chance flop pos compared to diff with global best
            for i in range(length):
                if pos[i] != beta[i]:
                    j = beta.index(pos[i])
                    beta[i], beta[j] = beta[j], beta[i]
                    if random.random() < b:
                        deltas.append((i, j))
            # flip changes found
            for i, j in deltas:
                pos[i], pos[j] = pos[j], pos[i]
            # update pbest
            pbest[num] = min(pos, pbest[num], key=locs.distance)

        yield gbest

def genetic(locs, select=100, size=100):
    """
    Arguments:
        locs (atlas): atlas type object
        select (int): how many members to keep after each generation
        size (int): the size of the population in each generation
    Yields:
        list: path found using gentic algorithm with SCX breeding
    """
    # intial random paths
    pop = [sample(range(len(locs)), len(locs)) for _ in range(size)]
    length = len(locs)

    def breed(par_a, par_b):
        # tsp path, last used node, set of not used locations
        last = randrange(length)
        path, legit = [last], set(range(length)) - {last}

        for i in range(length - 1):
            # find shorter first legitamate in either parent
            a = next((x for x in par_a[i:] if x in legit), *sample(legit, 1))
            b = next((x for x in par_b[i:] if x in legit), *sample(legit, 1))
            last = a if locs.dist[last][a] < locs.dist[last][b] else b

            # add to path and remove from legitamate
            legit.remove(last)
            path.append(last)

        # mutate
        a, b = sample(range(length), 2)
        path[a], path[b] = path[b], path[a]

        return path

    while True:
        # sort by fitness and yield best result
        pop = sorted(pop, key=locs.distance)[:select]
        yield pop[0]
        # breed those selected
        pop = [breed(*choices(pop, k=2)) for _ in range(size)]

def two_opt(locs, nn=False):
    """
    Arguments:
        locs (atlas): atlas type object
        nn (bool): initialize path as nearest neighbor greedy path
    Yields:
        list: the current path based on 2-opt heuristic
    """
    indexes = range(len(locs))
    path = greedy(locs) if nn else list(indexes)
    d = lambda x, y: locs.dist[path[x], path[y]]
    while True:
        change = True
        # generate all substrings of path
        for i, j in combinations(indexes, 2):
            yield path
            # flip substring if it shortens path
            if d(i-1, i) + d(j-1,j) > d(i-1,j-1) + d(i, j):
                path[i:j] = reversed(path[i:j])
                change = False

        # short circuit if no changes
        if change: return

def three_opt(locs, nn=False):
    """
    Arguments:
        locs (atlas): atlas type object
        nn (bool): initialize path as nearest neighbor greedy path
    Yields:
        list: the current path based on 3-opt heuristic
    """
    indexes = range(len(locs))
    path = greedy(locs) if nn else list(indexes)
    distance = lambda x, y: locs.dist[path[x], path[y]]
    pre = locs.distance(path)
    while True:
        # generate all two adjacent substrings of path
        for i, j, k in combinations(indexes, 3):
            yield path
            dist = distance(i-1,i) + distance(j-1,j) + distance(k-1,k)

            # modify those substrings to be shortest
            if dist > distance(i-1, j-1) + distance(i, j) + distance(k-1, k):
                path[i:j] = reversed(path[i:j])
            elif dist > distance(i-1, i) + distance(j-1, k-1) + distance(j, k):
                path[j:k] = reversed(path[j:k])
            elif dist > distance(k, i) + distance(j-1, j) + distance(k-1, i-1):
                path[i:k] = reversed(path[i:k])
            elif dist > distance(i-1, j) + distance(k-1, i) + distance(j-1, k):
                path[i:k] = path[j:k] + path[i:j]

        # short circuit if no changed made
        post = locs.distance(path)
        if pre == post:
            return
        else:
            pre = post

def sim_annealing(locs, steps=1000, preverse=0.5, temp=0.5, factor=0.999999):
    """
    Arguments:
        locs (atlas): atlas type object
        steps (int): how many steps of changes to do every loop
        preverse (float): prefer two opt over three opt factor
        temp (num): starting temperature
        factor (float): how much to degrade temperature each loop
    Yields:
        list: the current path based on simulated annealing algorithm
    """
    path = list(range(len(locs)))
    d = lambda x, y: locs.dist[path[x], path[y]]
    while True:
        changes = 0
        for i in range(steps):
            i, j, k = sorted(choices(path, k=3))

            # chance to do a two opt change
            if preverse > random.random():
                de = d(i - 1, j - 1) + d(i, j) - d(i - 1, i) - d(j - 1, j)

                if 0 > de or np.exp(-de/temp) > random.random():
                    path[i:j] = reversed(path[i:j])
                    changes += 1
            # chance to do a three opt change
            else:
                de = d(i, k-1) + d(j-1, k) + d(i-1, j) - \
                     d(i-1, i) - d(j-1, j) - d(k-1, k)

                if 0 > de or np.exp(-de/temp) > random.random():
                    path = path[i:j] + path[k:] + path[:i] + path[j:k]
                    changes += 1

        # update temperature return if it has stopped changing
        if changes == 0: return
        temp *= factor
        yield path


text:
import math

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def interpcurve(N, pX, pY):
    # Equally spaced in arclength
    N = np.transpose(np.linspace(0, 1, N))

    # How many points will be uniformly interpolated?
    nt = N.size

    # Number of points on the curve
    n = pX.size
    pxy = np.array((pX, pY)).T
    p1 = pxy[0, :]
    pend = pxy[-1, :]
    last_segment = np.linalg.norm(np.subtract(p1, pend))
    epsilon = 10 * np.finfo(float).eps

    # If the two end points are not close enough lets close the curve
    if last_segment > epsilon * np.linalg.norm(np.amax(abs(pxy), axis=0)):
        pxy = np.vstack((pxy, p1))
        nt = nt + 1
    else:
        print('Contour already closed')
    pt = np.zeros((nt, 2))

    # Compute the chordal arclength of each segment.
    chordlen = (np.sum(np.diff(pxy, axis=0)**2, axis=1))**(1 / 2)
    # Normalize the arclengths to a unit total
    chordlen = chordlen / np.sum(chordlen)
    # Cumulative arclength
    cumarc = np.append(0, np.cumsum(chordlen))
    tbins = np.digitize(N, cumarc)  # bin index in which each N is in

    # Catch any problems at the ends
    tbins[np.where(tbins <= 0 | (N <= 0))] = 1
    tbins[np.where(tbins >= n | (N >= 1))] = n - 1

    s = np.divide((N - cumarc[tbins]), chordlen[tbins - 1])
    pt = pxy[tbins, :] + np.multiply((pxy[tbins, :] - pxy[tbins - 1, :]),
                                     (np.vstack([s] * 2)).T)

    return pt

def haversine_dist(init_point, final_point):
    """Gives the Haversine distance between the initial and final point

    Parameters
    ----------
    init_point : array
        The initial point
    final_point : array
        The final point

    Returns
    -------
    float
        The Haversine distance
    """

    # Assuming the input is in degrees
    init_rad = init_point * math.pi / 180
    final_rad = final_point * math.pi / 180

    d_latitude = final_rad[0] - init_rad[0]
    d_longitude = final_rad[1] - init_rad[1]
    x = (d_longitude) * math.cos((final_rad[0] - init_rad[0]) / 2)
    y = d_latitude

    return [x * 6356.752e3, y * 6356.752e3]

def get_xy_position(config):
    """Get the x and y position of all the buildings
    """
    read_path = config['map_data_path'] + 'latitude_longitude.csv'
    data = np.genfromtxt(read_path, delimiter=',', skip_header=True)

    init_point = data[0]
    xy_pos = np.ones((data.shape[0], 3))
    for i, item in enumerate(data):
        temp_x = item.copy()
        xy_pos[i, 0:2] = haversine_dist(init_point, temp_x)

    # Rotation angle
    t = math.pi / 2 + math.atan2(114.2746 - 184.137, -67.67721 + 82.2549)
    rot = np.asarray([[math.cos(t), -math.sin(t), 0],
                      [math.sin(t), math.cos(t), 0], [0, 0, 1]])
    xy_pos = np.matmul(xy_pos, rot)

    plt.scatter(xy_pos[:, 0], xy_pos[:, 1])
    plt.show()
    df = pd.DataFrame(xy_pos, columns=['x', 'y'])
    filepath = config['map_save_path'] + 'co_ordinates.xlsx'
    df.to_excel(filepath, index=False)

    return None

max_stars_repo_name: JosephDistefano/offset-human-interface

max_issues_repo_name: JosephDistefano/offset-human-interface

max_forks_repo_name: JosephDistefano/offset-human-interface


text:
import numpy as np

from Algorithms.Preprocessing import Preprocessing

class PatchPreprocessing(Preprocessing):
    def __init__(self, data, labels, n_train_samples, n_validation_samples, use_ground_truth, patch_radius):
        super().__init__(data, labels, n_train_samples, n_validation_samples, use_ground_truth)

        self.patch_radius = patch_radius

    """
    input: (targets, labels)
    - targets: (y, x)-pairs representing pixels in the test set
    - labels: class labels for the respective samples in targets

    output: (test_data, labels)
    - test_data: preprocessed test data
    - labels: class labels for the respective samples in test_data
    """
    def model_get_test_batch(self, targets, labels):
        test_data = self.model_preprocess_pre()

        for (y, x) in targets:
            test_data = self.model_preprocess_one(test_data, x, y)

        return self.model_preprocess_post(test_data, labels, shuffle=False)

    """
    input: (x, y)
    - x: x coordinate of the pixel whose validity is being tested
    - y: y coordinate of the pixels whose validity is being tested

    output: (valid)
    - valid: whether this pixel can be preprocessed correctly (could be False for pixels on the edge)
    """
    def is_valid_training_target(self, x, y):
        if y < self.patch_radius or y > self.data.shape[0] - self.patch_radius - 1:
            return False

        if x < self.patch_radius or x > self.data.shape[1] - self.patch_radius - 1:
            return False

        return True

    """
    input: ()

    output: (data)
    - data: the structure that the preprocessed data will take, for example an array of three empty arrays for three datastreams
    """
    def model_preprocess_pre(self):
        return []

    """
    input: (data, x, y)
    - data: the current data array, new sample should be appended
    - x: x coordinate of the target pixel
    - y: y coordinate of the target pixel

    output: (data)
    - data: the new data array, now including the sample with coordinates (y, x)
    """
    def model_preprocess_one(self, data, x, y):
        data.append(self.data[
            y - self.patch_radius:y + self.patch_radius + 1,
            x - self.patch_radius:x + self.patch_radius + 1,
            :
        ])

        return data

    """
    input: (targets, data_train, data_validation)
    - targets: (y, x)-pairs valid for being selected for train and validation set
    - data_train: the structure returned by model_preprocess_pre for the training data, new samples should be appended
    - data_validation: the structure returned by model_preprocess_pre for the validation data, new samples should be appended

    output: (data_train, data_validation)
    - data_train: the filled data structure with training data
    - data_validation: the filled data structure with validation data
    """
    def model_preprocess(self, targets, data_train, data_validation):
        j = 0
        while j < self.n_train_samples + self.n_validation_samples:
            (y, x) = targets[min(j, len(targets) - 1)]

            if j < min(self.n_train_samples, len(targets) / 2.0):
                data_train = self.model_preprocess_one(data_train, x, y)
            elif j >= min(self.n_train_samples, len(targets) / 2.0) and j < self.n_train_samples:
                # repick an earlier training sample
                chosen_sample = np.random.randint(len(data_train))
                data_train.append(data_train[chosen_sample])
            elif j < self.n_train_samples + min(self.n_validation_samples, len(targets) / 2.0):
                data_validation = self.model_preprocess_one(data_validation, x, y)
            else:
                # repick an earlier validation sample
                chosen_sample = np.random.randint(len(data_validation))
                data_validation.append(data_validation[chosen_sample])

            j += 1

        return (data_train, data_validation)

    """
    input: (data, labels, shuffle)
    - data: the preprocessed data structure
    - labels: the corresponding labels for the data structure
    - shuffle: whether the data (and labels) should be shuffled

    output: (data, labels)
    - data: the postprocessed data structure
    - labels: the postprocessed labels, optionally shuffled in the same way as the data
    """
    def model_preprocess_post(self, data, labels, shuffle=False):
        if labels is not None:
            labels = np.array(labels)

        data = np.array(data)

        if shuffle:
            indices = np.arange(len(labels))
            np.random.shuffle(indices)

            if labels is not None:
                labels = labels[indices]

            data = data[indices]

        return (data, labels)

    """
    input: ()
    output: (shape)
    - shape: the shape of all the data after preprocessing
    """
    def get_all_data_shape(self):
        return [self.data.shape[0] - 2 * self.patch_radius, self.data.shape[1] - 2 * self.patch_radius]

max_stars_repo_name: zenodhaene/HyperspectralClassification

max_issues_repo_name: zenodhaene/HyperspectralClassification

max_forks_repo_name: zenodhaene/HyperspectralClassification


text:
import numpy as np
import os as os

from YGRW.data_interp import JumpDistFromAngle, AngleFromAngle

TEST_DIR = os.path.dirname(__file__)
DATA_DIR = os.path.join(TEST_DIR, "..", "data")

def test_open_data():
    """
    Test integrity of loading / unloading data
    Returns
    -------

    """

    isotropic_data_path = os.path.join(DATA_DIR, "jump_distances_isotropic.csv")
    data = np.loadtxt(isotropic_data_path, skiprows=1, usecols=[1, 2], delimiter=",")
    assert np.array_equal(data.shape, (512, 2))

    angle_corr_data_path = os.path.join(DATA_DIR, "angle_correlation.csv")
    data = np.loadtxt(angle_corr_data_path, skiprows=1, usecols=[1, 2], delimiter=",")
    assert np.array_equal(data.shape, (512, 2))

    angle_weighted_path = os.path.join(DATA_DIR, "jump_distances_angleweighted.csv")
    use_range = list(range(38))
    data = np.loadtxt(
        angle_weighted_path, skiprows=1, usecols=use_range[1:], delimiter=","
    )
    assert np.array_equal(data.shape, (500, 37))

    angle_weighted_path = os.path.join(DATA_DIR, "successive_angle_probability.csv")
    use_range = list(range(38))
    data = np.loadtxt(
        angle_weighted_path, skiprows=1, usecols=use_range[1:], delimiter=","
    )
    print(data.shape)
    assert np.array_equal(data.shape, (512, 37))

def test_jump_dist_from_angle():

    jdfa = JumpDistFromAngle()

    result = jdfa.distribution_from_angle(2.5)
    assert len(result) == 500

    import matplotlib.pyplot as plt

    import matplotlib.cm as cm
    from matplotlib.colors import Normalize

    map = cm.get_cmap("RdBu")
    norm = Normalize(vmin=0, vmax=180)

    for i, theta in enumerate(np.arange(2.5, 180, 5)):
        data_vals = jdfa.angle_data[:, i]
        assert np.allclose(
            data_vals / np.sum(data_vals), jdfa.distribution_from_angle(theta)
        )

    """
    for i, theta in enumerate(np.arange(0, 180, 1)):

        vals = jdfa.distribution_from_angle(theta)
        plt.plot(jdfa.jump_values, vals/ np.sum(5*vals), alpha=.1,
                label=str(theta),color=map(norm(theta)))
        plt.title("Step Size | Angle of New Step")
    plt.show()
    """

def test_angle_from_angle():

    afa = AngleFromAngle()

    result = afa.distribution_from_angle(2.5)
    assert len(result) == 512

    for i, theta in enumerate(np.arange(2.5, 180, 5)):

        data_vals = afa.angle_data[:, i]
        vals = afa.distribution_from_angle(theta)
        assert np.allclose(data_vals / np.sum(data_vals), vals)

    """
    import matplotlib.pyplot as plt
    import matplotlib.cm as cm
    from matplotlib.colors import Normalize

    map = cm.get_cmap("RdBu")
    norm = Normalize(vmin=0, vmax=180)
    for i, theta in enumerate(np.arange(2.5, 180, 1)):

        vals = afa.distribution_from_angle(theta)
        plt.plot(afa.next_angles, vals/ np.sum(5*vals), alpha=.1,
                 label=str(theta),color=map(norm(theta)))
        plt.title("Angle Step | Angle Step")
        #plt.legend()
    plt.show()
    """

max_stars_repo_name: MCnu/YGRW


text:
"""
The data creation script. Assumes a certain folder structure for the attacks (best viewed in the `paths` config (
`attacks_path_template`.
"""
import logging
import os
from collections import namedtuple

import hydra
import numpy as np
import torch
import zarr
from omegaconf import DictConfig, OmegaConf
from torch.utils.data import DataLoader

from daain.backbones.esp_net.espnet import CustomESPNet
from daain.config_schema import PATH_TO_CONFIGS
from daain.data import PerturbedDataset
from daain.data.activations_dataset import compute_activations, create_mask
from daain.utils.pytorch_utils import enforce_reproducibility
from daain.utils.zarr_utils import open_group

# def _get_transforms(dataset_dict):
#    return T.Compose([getattr(T, n)(**kwargs) for n, kwargs in dataset_dict["transforms"].items()])

def _get_loader(dataset_dict, attacks_path_template):
    """Source"""
    root_path_to_attacked_images = attacks_path_template.format(
        data_split=dataset_dict["split"], dataset=dataset_dict["dataset"]
    )

    data = PerturbedDataset(
        attacks=dataset_dict["attacks"], root_path=root_path_to_attacked_images, transforms=dataset_dict["transforms"]
    )
    return DataLoader(dataset=data, batch_size=1, shuffle=False)

def _get_activations_group(dataset_dict, activations_path_template, force=False):
    """Sink"""
    activations_path = activations_path_template.format(
        data_split=dataset_dict["split"], dataset=dataset_dict["dataset"]
    )
    return zarr.group(store=zarr.DirectoryStore(activations_path), overwrite=force)

def _gen_masks(
    mask_name, rejection_boundaries, activation_space_shape, min_distance, activations_group, masks_zarr, force=False
):
    group = open_group(activations_group, mask_name)

    group_raw = open_group(group, "raw")
    group_gridified = open_group(group, f"gridified")

    if force or f"{mask_name}_raw" not in masks_zarr:
        logging.info(f"using new mask for {mask_name}")
        t = create_mask(min_distance, rejection_boundaries, activation_space_shape)
        if rejection_boundaries is None:
            masks_zarr[f"{mask_name}_raw"] = t[0]
            masks_zarr[f"{mask_name}_gridified"] = t[1]
        else:
            masks_zarr[f"{mask_name}_raw"] = t
    else:
        logging.info(f"reusing existing mask {mask_name}")

    if rejection_boundaries is None:
        return zip(
            [np.array(masks_zarr[f"{mask_name}_{s}"]) for s in ["raw", "gridified"]], [group_raw, group_gridified]
        )
    else:
        return zip([np.array(masks_zarr[f"{mask_name}_{s}"]) for s in ["raw",]], [group_raw,])

def _gen_all_sampling_masks(activations_group, masks_zarr, mask_descriptions, force=False):
    """Used to combine all masks such that we load the image and thus activations only once. Then sample from all
    recorded activations in one go."""
    for mask in mask_descriptions:
        for selection, zarr_group in _gen_masks(
            mask_name=mask.mask_name,
            rejection_boundaries=mask.rejection_boundaries,
            activation_space_shape=mask.activation_space_shape,
            min_distance=mask.min_distance,
            activations_group=activations_group,
            masks_zarr=masks_zarr,
            force=force,
        ):
            yield selection, zarr_group, mask.recording_type

def generate_activations(
    model,
    dataset_descriptions,
    masks_zarr,
    mask_descriptions,
    attacks_path_template,
    activations_path_template,
    force=False,
    fast_dev_run=False,
):
    for dataset_dict in dataset_descriptions:
        loader = _get_loader(dataset_dict, attacks_path_template=attacks_path_template)
        activations_group = _get_activations_group(
            dataset_dict, activations_path_template=activations_path_template, force=force
        )

        logging.info(f"computing activations for {dataset_dict['split']}, {dataset_dict['dataset']}")
        selections = list(
            _gen_all_sampling_masks(
                activations_group=activations_group, masks_zarr=masks_zarr, mask_descriptions=mask_descriptions
            )
        )
        compute_activations(
            model,
            loader,
            attacks=dataset_dict["attacks"],
            selections=selections,
            force=force,
            n_data_points=10 if fast_dev_run else -1,
        )

logger = logging.getLogger(__name__)

# TODO this is causing some issues with typing a list of lists
# cs = ConfigStore.instance()
#
# cs.store(group="backbone", name="esp_dropout", node=ESPDropoutModelTrainingConfigSchema)
# cs.store(group="dataset_paths", name="cityscapes", node=DatasetPaths)
# cs.store(group="dataset_paths", name="bdd100k", node=DatasetPaths)
# cs.store(group="activation_spaces", name="esp_net_256_512", node=ActivationSpace)
# cs.store(name="data_creation", node=DataCreationSchema)

@hydra.main(config_path=PATH_TO_CONFIGS, config_name="data_creation")
# def main(cfg: DataCreationSchema):
def main(cfg: DictConfig):
    logging.info(OmegaConf.to_yaml(cfg))

    # for the poisson disc sampling algorithm
    enforce_reproducibility(cfg.random_seed)
    root_dir = cfg.paths["derived"]
    masks_zarr = zarr.group(store=zarr.DirectoryStore(os.path.join(cfg.paths["activations"], "masks.zarr")))
    activations_path_template = os.path.join(root_dir, cfg.paths["activation_path_template"])
    attacks_path_template = os.path.join(root_dir, cfg.paths["attacks_path_template"])

    MaskDescription = namedtuple(
        "MaskDescription",
        ["mask_name", "min_distance", "recording_type", "rejection_boundaries", "activation_space_shape"],
    )
    masks = [
        MaskDescription(
            **{
                **kwargs,
                "activation_space_shape": cfg.activation_spaces.activation_space_shapes[kwargs["recording_type"]],
            }
        )
        for kwargs in cfg.masks
    ]

    backbone_model = CustomESPNet(
        device=torch.device(cfg.device), p=cfg.backbone.model.p, q=cfg.backbone.model.q, num_classes=cfg.num_classes
    ).load_pretrained_weights(cfg.paths.model_weights)

    generate_activations(
        backbone_model,
        masks_zarr=masks_zarr,
        dataset_descriptions=cfg.dataset_descriptions,
        attacks_path_template=attacks_path_template,
        activations_path_template=activations_path_template,
        mask_descriptions=masks,
    )

if __name__ == "__main__":
    main()

max_stars_repo_name: merantix/mxlabs-daain

max_issues_repo_name: merantix/mxlabs-daain

max_forks_repo_name: merantix/mxlabs-daain


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed May 29 16:53:19 2019

@author: gparkes
"""

from numba import jit
import numpy as np
import matplotlib.pyplot as plt

# task 1
@jit(nopython=True)
def calc_positional_e(spins, J):
    """
    In this strategy, we use pure Python and beef up our code with JIT or Cython.
    """
    N = spins.shape[0]
    E = np.zeros_like(spins)
    # for loop
    for i in range(N):
        for j in range(N):
            if i == 0:
                spin_l = spins[N-1,j]
                spin_r = spins[i+1,j]
            elif i == N-1:
                spin_l = spins[i-1,j]
                spin_r = spins[0,j]
            else:
                spin_l = spins[i,j-1]
                spin_r = spins[i,j+1]

            if j == 0:
                spin_u = spins[i,N-1]
                spin_d = spins[i,j+1]
            elif j == N-1:
                spin_u = spins[i,j-1]
                spin_d = spins[i,0]
            else:
                spin_u = spins[i,j-1]
                spin_d = spins[i,j+1]

            E[i,j] = spin_l + spin_r + spin_u + spin_d
    return E

# task 2
@jit
def delta_e(spins, i, j, J):
    N = spins.shape[0]
    # check boundaries first!
    if i == 0:
        spin_l = spins[N-1,j]
        spin_r = spins[i+1,j]
    elif i == N-1:
        spin_l = spins[i-1,j]
        spin_r = spins[0,j]
    else:
        spin_l = spins[i-1,j]
        spin_r = spins[i+1,j]

    if j == 0:
        spin_u = spins[i,N-1]
        spin_d = spins[i,j+1]
    elif j == N-1:
        spin_u = spins[i,j-1]
        spin_d = spins[i,0]
    else:
        spin_u = spins[i,j-1]
        spin_d = spins[i,j+1]

    return 2. * J * spins[i,j] * (spin_l+spin_r+spin_u+spin_d)

# task 3
def total_mag_e(spins, J):
    return -J*np.sum(calc_positional_e(spins, J))

# task 4
def metropolis_hastings(N, n_steps, beta):
    spins = np.random.choice([-1,1], size=(N,N))
    J = 1
    E_series = np.ones(n_steps)
    E_series[0] = total_mag_e(spins, J)

    for s in range(1, n_steps):
        randx = np.random.randint(N)
        randy = np.random.randint(N)
        # calculate de
        dE = delta_e(spins, randx, randy, J)
        # if we decrease the energy or meet boltzmann requirements, flip the spin
        if dE < 0. or np.exp(-beta*dE) > np.random.rand():
            spins[randx, randy] *= -1
            E_series[s] = E_series[s-1] + dE
        else:
            E_series[s] = E_series[s-1]
    return spins, E_series

# task 5
def task_5():
    A = metropolis_hastings(40, int(5*10**5), 0.1)
    B = metropolis_hastings(40, int(5*10**5), 1.0)

    fig,ax=plt.subplots(ncols=3, figsize=(16,4))

    ax[0].plot(A[1], label=r"$\beta=0.1$")
    ax[0].plot(B[1], label=r"$\beta=1.0$")
    ax[0].legend()
    ax[1].imshow(A[0], cmap="gray")
    ax[2].imshow(B[0], cmap='gray')

    for a in [ax[1], ax[2]]:
        a.get_xaxis().set_visible(False)
        a.get_yaxis().set_visible(False)

    ax[0].set_xlabel("steps")
    ax[0].set_ylabel(r"Energy $E$")
    plt.show()

# task 6
def task_6():
    bvals = 50
    betas = np.linspace(0.1, 0.6, bvals)
    m_ser = np.zeros(bvals)

    def average_magnetization(spins, N):
        return (1 / N**2) * np.sum(spins)

    for b in range(bvals):
        # print("Running beta=%.4f" % betas[b])
        C, E = metropolis_hastings(20, 500000, betas[b])
        m_ser[b] = average_magnetization(C[0], 20)

    # plot
    fig = plt.figure(figsize=(13,5))
    plt.plot(betas, m_ser, "g*")
    plt.xlabel(r"$\beta$")
    plt.ylabel(r"$M$")
    plt.show()

# task 7

"""
Uncomment this code below to run in Jupyter notebook
"""
# %prun metropolis_hastings(40, 500000, 0.5)

# %timeit metropolis_hastings(40, 500000, 0.5)

# after running cython block...

# %timeit metropolis_hastings2(40, 500000, 0.5)

def task_7():
    C1, E1 = metropolis_hastings(40, 500000, .6)
    C2, E2 = metropolis_hastings2(40, 500000, .6)
    C3, E3 = metropolis_hastings(40, 500000, .1)
    C4, E4 = metropolis_hastings2(40, 500000, .1)
    fig,ax=plt.subplots(ncols=2, figsize=(16,4))

    ax[0].plot(E1, label=r"python")
    ax[0].plot(E2, label="cython")
    ax[0].legend()
    ax[1].plot(E3, label=r"python")
    ax[1].plot(E4, label="cython")
    ax[1].legend()

max_stars_repo_name: gregparkes/PythonTeaching

max_issues_repo_name: gregparkes/PythonTeaching

max_forks_repo_name: gregparkes/PythonTeaching

reason: import numpy,from numba