text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
% Copyright 2020 Niccolò Marcucci <niccolo.marcucci@polito.it>
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
"""

import meep as mp
import numpy as np
# import h5py as h5
from scipy.io import loadmat, savemat
from math import pi, cos, sin, tan, atan, atan2, sqrt
from matplotlib import pyplot as plt, cm, colors, widgets

from meep_objects.python_utils import *

def anisotropic_material (index, anisotropy1, anisotropy2=0, rot_angle_3=0,
                          rot_angle_1=0, rot_angle_2=0 ):
    """
    Creates and anysotropyc material with two extraordinary axis and one
    ordinary axis (the third one). The axis are assumed to be orthogonal.
    If all rotation angles are equal to 0 first axis would be x, second would
    be y and third z.

    Parameters
    ----------
    index : Real Positive
        Base refractive index.
    anisotropy1 : Real Positive
        Relative variation of the refractive index along the first axis with
        respect to the Base value, expressed as a percentage.
    anisotropy2 : Real Positive, optional
        Same as anisotropy1. The default is 0.
    rot_angle_3 : Real Positive, optional
        Rotation AROUND the third axis. The default is 0.
    rot_angle_1 : Real Positive, optional
        Rotation AROUND the third axis. The default is 0.
    rot_angle_2 : Real Positive, optional
        Rotation AROUND the third axis. The default is 0.

    Returns
    -------
    meep.Medium()
        DESCRIPTION.

    """

    if anisotropy1 == 0 and anisotropy2 == 0:
        return mp.Medium(index=index)

    eps = index**2

    eps_e1 = eps*(1+anisotropy1/100)
    eps_e2 = eps*(1+anisotropy2/100)
    eps_matrix = np.diag(np.array([eps_e1, eps_e2, eps]))

    # first rotation around z axis
    rot_matrix = np.array([[cos(rot_angle_3), -sin(rot_angle_3) , 0],
                           [sin(rot_angle_3),  cos(rot_angle_3) , 0],
                           [0                  ,  0             , 1]])
    eps_matrix = np.matmul(np.matmul(rot_matrix,eps_matrix),
                           np.linalg.inv(rot_matrix))

    # second rotation around x axis
    rot_matrix = np.array([[1 , 0                  , 0             ],
                           [0 , cos(rot_angle_1), -sin(rot_angle_1)],
                           [0 , sin(rot_angle_1),  cos(rot_angle_1)]])
    eps_matrix = np.matmul(np.matmul(rot_matrix,eps_matrix),
                           np.linalg.inv(rot_matrix))

    # second rotation around y axis
    rot_matrix = np.array([[ cos(rot_angle_2), 0 , sin(rot_angle_2)],
                           [0                , 1 , 0               ],
                           [-sin(rot_angle_2), 0 , cos(rot_angle_2)]])
    eps_matrix = np.matmul(np.matmul(rot_matrix,eps_matrix),
                           np.linalg.inv(rot_matrix))

    return mp.Medium(epsilon_diag = np.diag(eps_matrix),
                     epsilon_offdiag = np.array([eps_matrix[0,1],
                                                 eps_matrix[0,2],
                                                 eps_matrix[1,2]]))

def circular_DBR_cavity (medium_back=mp.Medium(epsilon=1),
                         medium_groove=mp.Medium(epsilon=2),
                         D=0.4, DBR_period=0.2, FF=0.5, N_rings=10,
                         thickness = 0, orientation = mp.Vector3(0,0,1)):
    """
    Circular DBR cavity created as a sequence of concentric cylinders

    Parameters
    ----------
    medium_back : mp.Medium(),
        Background medium. The default is Vacuum.
    medium_groove : mp.Medium(),
        Medium of the rings. The default is mp.Medium(epsilon=2).
    D : Real Positive
        Diameter of the centre of the cavity. The default is 0.4.
    DBR_period : Real Positive
        Radial period of the DBR. The default is 0.2.
    FF : Real, optional
        Fill factor defined as a number between 0 and 1. The default is 0.5.
    N_rings : Integer
        Number of DBR rings. The default is 10.

    Returns
    -------
    List
        Returns a list of (meep) geometric objects.

    """
    device = []
    rings = []
    for i in range(1,N_rings+1):
        c1 = mp.Cylinder(radius = D/2+(N_rings-i+FF)*DBR_period,
                         height = thickness,
                         axis = orientation,
                         material = medium_groove)
        c2 = mp.Cylinder(radius = D/2+(N_rings-i)*DBR_period,
                         height = thickness,
                         axis = orientation,
                         material = medium_back)
        rings.append(c1)
        rings.append(c2)
    device.extend(rings)
    return device

def dipole_source(f, df=0, source_pos=mp.Vector3(0,0,0),
                  theta=0, phi=0, amplitude = 1):
    """
    Supports only cylindrical coordinates for now

    Parameters
    ----------
    f : Real
        Centre frequency of the source Gaussian spectrum.
    df : Real
        Width of the source Gaussian spectrum. If set to 0 the source will not
        be Gaussian but a continous sinusoidal source.
        The default is 0.
    source_pos : meep.Vector3().
        Source posistion. The default is meep.Vector3(0,0,0).
    theta : Real, optional
        Inclination of the dipole moment with respect to z axis, in radiants.
        The default is 0.
    phi : Real, optional
        Azimutal orientation of the dipole moment, in radiants.
        e.g. theta = pi/2, phi = 0 implies the dipole moment oriented along y.
        The default is 0.
    amplitude: Complex, optional
        Complex amplitude multiplying the current source. Can impart a phase
        shift to the real part of the overall current.
        The default is 1.

    Returns
    -------
    list
        A list containing the three sources (meep objects) required for
        orienting a source arbitrarily in space.

    """

    if df == 0 :
        source_type = mp.ContinuousSource(f,width=0.1 )
    else :
        source_type = mp.GaussianSource(f,df)

    source_x = mp.Source(source_type,
                         component=mp.Ex,
                         amplitude=cos(phi + pi/2)*sin(theta)*amplitude,
                         center=source_pos,)

    source_y = mp.Source(source_type,
                         component=mp.Ey,
                         amplitude=sin(phi + pi/2)*sin(theta)*amplitude,
                         center=source_pos)

    source_z = mp.Source(source_type,
                         component=mp.Ez,
                         amplitude = cos(theta)*amplitude,
                         center=source_pos,)
    return [source_x,source_y, source_z]

def plane_wave_source(f, df, k, center, size, inc_plane_norm, amplitude=1):
    """
    Plane wave source.
    """
    uk = k.unit()
    inc_plane_norm = inc_plane_norm.unit()

    if uk.dot(inc_plane_norm) != 0.:
        # The dot product might never return 0 if inc_plane_norm i not one of
        # the axis due to numerical error.
        raise ValueError("the wavevector has to be orthogonal to the incidence plane")

    uE = uk.cross(inc_plane_norm)
    E = uE * amplitude

    def pw_amp(E, k, x0, component):
        def _pw_amp(x):
            if component == mp.Ex :
                ampl = E.x
            elif component == mp.Ey :
                ampl = E.y
            elif component == mp.Ez :
                ampl = E.z

            return ampl * np.exp( 1j*k.dot(x-x0) )

        return _pw_amp

    plane_wave = [mp.Source(mp.ContinuousSource(f,fwidth=0.1,is_integrated=True) if df==0 else mp.GaussianSource(f,fwidth=df,is_integrated=True),
                            component = mp.Ex,
                            center = center,
                            size = size,
                            amp_func = pw_amp(E, k, center, mp.Ex)),
                  mp.Source(mp.ContinuousSource(f,fwidth=0.1,is_integrated=True) if df==0 else mp.GaussianSource(f,fwidth=df,is_integrated=True),
                            component = mp.Ey,
                            center = center,
                            size = size,
                            amp_func = pw_amp(E, k, center, mp.Ey)),
                  mp.Source(mp.ContinuousSource(f,fwidth=0.1,is_integrated=True) if df==0 else mp.GaussianSource(f,fwidth=df,is_integrated=True),
                            component = mp.Ez,
                            center = center,
                            size = size,
                            amp_func = pw_amp(E, k, center, mp.Ez))]
    return plane_wave

def dielectric_multilayer(design_file, substrate_thickness, x_width,
                          y_width=1, unit = 'um', used_layer_info={}, exclude_last_layer=False,
                          buried=False, axis = mp.Z) :
    """
    Dielectric multilayer stack
    """

    data = loadmat(design_file)

    N = np.size(data['idx_layers'])
    idx_layers = data['idx_layers'].reshape(N)
    d_layers   = data['d_layers'].reshape(N)
    # print(used_layer_info)
    # print(d_layers)
    d_layers[ used_layer_info['used_layer'] ] = used_layer_info['thickness']*1e-6
    idx_layers[ used_layer_info['used_layer'] ] = used_layer_info['refractive index']
    used_layer = used_layer_info['used_layer']

    if unit == 'nm' :
        d_layers *= 1e9
    elif unit == 'um' :
        d_layers *= 1e6
    else :
        raise ValueError("Unit not supported")

    d_layers[ 0] = substrate_thickness
    d_layers[-1] = 0

    multilayer = []

    for i, d in enumerate(d_layers) :
        if d == 0 :
            continue

        z = np.sum( d_layers[:i] ) + d/2 - np.sum(d_layers[:used_layer]) - d_layers[used_layer]/2
        z = np.round(z, 4)
        if not (i == N-2 and exclude_last_layer) :
            if axis == mp.Z :
                size = [x_width, y_width, d]
                centre = [0, 0, z]
            elif axis == mp.Y:
                size = [x_width, d, y_width]
                centre = [0, z, 0]
            elif axis == mp.Y:
                size = [d, y_width, x_width]
                centre = [z, 0, 0]

            multilayer.append(mp.Block(material = mp.Medium(index = np.real(idx_layers[i])),
                                       size     = mp.Vector3(*size),
                                       center   = mp.Vector3(*centre)))
            multilayer[i].name  = f'Layer_{i}'
            multilayer[i].group = 'Multilayer'

    thickness = np.sum(d_layers[1:])

    design_specs = {"d_layers"  :d_layers,
                    "idx_layers":idx_layers }
    return multilayer, thickness, design_specs

def my_sum(array, i, j):
    # equivalent for sum(array[i:j]), valid also for j<i
    # j is excluded from the sum, as per stadard python
    i, j = [array.size + k if k < 0 else k for k in [i, j] ]
    if i <= j:
        return +np.sum(array[i:j])
    else:
        return -np.sum(array[j:i])

def elliptic_DBR_cavity(medium_groove=mp.Medium(epsilon=2),
                        D=0.4, d=0.3, DBR_period=0.2, FF=0.5, N_rings=10,
                        thickness=0,
                        center=mp.Vector3(), axial_rotation=0):
    """
    Elliptic DBR cavity created as a sequence of concentric cylinders.

    Parameters
    ----------
    medium_groove : mp.Medium(),
        Medium of the rings. The default is mp.Medium(epsilon=2).
    D : Real Positive
        Major axis diametre of the centre of the cavity. The default is 0.4.
    d : Real Positive
        Miinor axis diameter of the centre of the cavity. The default is 0.3.
    DBR_period : Real Positive
        Radial period of the DBR. The default is 0.2.
    FF : Real, optional
        Fill factor defined as a number between 0 and 1. The default is 0.5.
    N_rings : Integer
        Number of DBR rings. The default is 10.
    axial_rotation : radians
        In plane rotation of the ellipse.

    Returns
    -------
    List
        Returns a list of (meep) geometric objects.

    """
    rings = []
    polygons = grating_veritices(DBR_period, D/2, d/2, N_rings, n_arms = 0,  FF = FF)
    for ring in polygons:
        vertices = [mp.Vector3(v[0]*cos(axial_rotation)-v[1]*sin(axial_rotation),
                               v[0]*sin(axial_rotation)+v[1]*cos(axial_rotation),
                               0)   for v in np.transpose(ring)]
        centroid =  sum(vertices, mp.Vector3(0)) * (1.0 / len(vertices))

        c1 = mp.Prism(vertices = vertices,
                      height = thickness,
                      axis = mp.Vector3(z=1),
                      center = center + centroid,
                      material = medium_groove)
        rings.append(c1)
    return rings

def linear_DBR_cavity(medium_groove=mp.Medium(epsilon=2),
                      D=0.4, DBR_period=0.2, FF=0.5, N_periods=10,
                      width=1, thickness=0, axis=mp.Z,
                      axial_rotation=0, center=mp.Vector3()):
    """
    Linear DBR cavity created as a sequence of rectangles.

    Parameters
    ----------
    medium_groove : mp.Medium(),
        Medium of the rings. The default is mp.Medium(epsilon=2).
    D : Real Positive
        spacer at the centre of the cavity. The default is 0.4.
    DBR_period : Real Positive
        period of the DBR. The default is 0.2.
    FF : Real, optional
        Fill factor defined as a number between 0 and 1. The default is 0.5.
    N_periods : Integer
        Number of DBR lines on each side of the spacer. The default is 10.
    width : Real, optional
        Length of each groove. Default is 10um
    axial_rotation : radians
        In plane rotation of the DBR.

    Returns
    -------
    List
        Returns a list of (meep) geometric objects.

    """
    device = []
    for i in range(1,N_periods+1):
        groove_size = mp.Vector3(FF*DBR_period, width, thickness)
        groove_centre = mp.Vector3(D/2 + FF*DBR_period/2 + (N_periods-i)*DBR_period, 0, 0)

        phi = axial_rotation
        e1 = mp.Vector3(cos(phi), sin(phi), 0)
        e2 = mp.Vector3(cos(phi+pi/2), sin(phi+pi/2), 0)
        e3 = mp.Vector3(z=1)
        groove_centre_right = groove_centre.rotate(mp.Vector3(z=1), phi)
        groove_centre_left  = groove_centre.rotate(mp.Vector3(z=1), phi + pi)

        if axis == mp.Y :
            e1 = e1.rotate(mp.Vector3(x=1), -pi/2)
            e2 = e2.rotate(mp.Vector3(x=1), -pi/2)
            e3 = e3.rotate(mp.Vector3(x=1), -pi/2)
            groove_centre_right = groove_centre_right.rotate(mp.Vector3(x=1), -pi/2) + center
            groove_centre_left = groove_centre_left.rotate(mp.Vector3(x=1), -pi/2) + center
        elif axis == mp.X :
            e1 = e1.rotate(mp.Vector3(y=1), -pi/2)
            e2 = e2.rotate(mp.Vector3(y=1), -pi/2)
            e3 = e3.rotate(mp.Vector3(y=1), -pi/2)
            groove_centre_right = groove_centre_right.rotate(mp.Vector3(y=1), -pi/2) + center
            groove_centre_left = groove_centre_left.rotate(mp.Vector3(y=1), -pi/2) + center
        elif axis == mp.Z :
            groove_centre_right += center
            groove_centre_left += center
        else:
            raise ValueError("Wrong axis value.")

        b1 = mp.Block(e1 = e1,
                      e2 = e2,
                      e3 = e3,
                      size = groove_size,
                      center = groove_centre_right,
                      material = medium_groove)
        b2 = mp.Block(e1 = e1,
                      e2 = e2,
                      e3 = e3,
                      size = groove_size,
                      center = groove_centre_left,
                      material = medium_groove)

        device.append(b1)
        device.append(b2)
    return device

def metasurface_radial_grating(medium_groove=mp.Medium(epsilon=2),
                               D=2, metasurface_period=0.2, scatter_length = 0.4,
                               scatter_width=0.5, scatter_tilt=pi/4,
                               scatter_type = 'radial', N_rings=10, N_arms=0,
                               thickness=0, orientation=mp.Vector3(0, 0, 1)):
    """
    Complex metasurface-like grating.

    Parameters
    ----------
    medium_groove : mp.Medium(),
        Medium of the rings. The default is mp.Medium(epsilon=2).
    D : Real Positive
        Diameter of the centre of the metasurface. The default is 2.
    metasurface_period : Real Positive
        Radial period of the DBR. The default is 0.2.
    scatter_width : TYPE, optional
        DESCRIPTION. The default is 0.5.
    scatter_tilt : TYPE, optional
        DESCRIPTION. The default is pi/4.
    N_rings : Integer
        Number of DBR rings. The default is 10.
    N_arms : TYPE, optional
        DESCRIPTION. The default is 0.
    thickness : TYPE, optional
        DESCRIPTION. The default is 0.
    orientation : TYPE, optional
        DESCRIPTION. The default is mp.Vector3(0, 0, 1).

    Raises
    ------
    ValueError
        DESCRIPTION.

    Returns
    -------
    None.

    """
    metasurface = []
    start_radius = D / 2
    if N_arms != 0:
        print('To be Implemented!')
        # for j in range(np.abs(N_arms)):
        #     b1 = mp.Block(e1=mp.Vector3(1,0,0).rotate(mp.Vector3(0,0,1),rotation),
        #      e2=mp.Vector3(0,1,0).rotate(mp.Vector3(0,0,1),rotation),
        #      size=groove_size,
        #      center=groove_centre.rotate(mp.Vector3(0,0,1),rotation),
        #      material = medium_groove)

    else:
        if N_rings != 0:
            for n in range(N_rings):
                radius = D/2 + n * metasurface_period

                if scatter_type == "filled":
                    N_scatters = round(2*pi * radius / scatter_length)
                elif scatter_type == "radial":
                    N_scatters = round(2*pi * start_radius / scatter_length)
                elif scatter_type == "radial_filled":
                    N_scatters = round(2*pi * start_radius / scatter_length)
                else:
                    raise ValueError()

                for k in range(N_scatters):
                    location_tilt = 2*pi / N_scatters * k

                    tilt = (pi/2 + scatter_tilt + location_tilt + n*pi/6)

                    metasurface.append(mp.Block(
                            e1=mp.Vector3(cos(tilt), sin(tilt), 0),
                            e2=mp.Vector3(cos(tilt+pi/2), sin(tilt+pi/2), 0),
                            size=mp.Vector3(scatter_length*0.75, scatter_width, thickness),
                            center=mp.Vector3(
                                radius * cos(location_tilt),
                                radius * sin(location_tilt),
                                0),
                            material=medium_groove))

                    metasurface[-1].name  = f'Scatter_{n}_{k}'
                    metasurface[-1].group = 'Metasurface'

    return metasurface

def linear_pol_splitting_grating(  medium_groove=mp.Medium(epsilon=2),
                                metasurface_period=0.4, scatter_length = 0.4,
                                scatter_width=0.1, scatter_tilt=pi/3,
                                scatter_shape = '',
                                N_periods_x=9, N_periods_y=1,
                                thickness=1, center=mp.Vector3(0, 0, 0)):
    """
    Circular polarization sensitive grating. Similar tu the metasurface-like
    grating, but with specific purpose.
    """
    sc_width = scatter_width
    sc_length = scatter_length

    # following lists are for v-shaped scatters
    rel_width = sc_width / sc_length
    scatter_vertices_x = [0, .5, .5-rel_width/0.866, 0, -.5+rel_width*0.866, -.5]
    scatter_vertices_y = [-.5, .5, .5, -.5+rel_width/0.5, .5, .5]
    scatter_vertices = np.array([scatter_vertices_x, scatter_vertices_y]).transpose()

    metasurface = []

    if N_periods_x != 0:
        for n in range(N_periods_x):
            for k in range(N_periods_y):
                tilt = n * scatter_tilt + pi/2

                if scatter_shape == 'V' or scatter_shape == 'v' :
                    tilt = tilt - pi
                    vertices = [mp.Vector3(v[0],v[1],0).rotate(mp.Z,tilt)*sc_length for v in scatter_vertices]
                    centroid =  sum(vertices, mp.Vector3(0)) * (1.0 / len(vertices))

                    scatter =  mp.Prism(vertices = vertices,
                                        height = thickness,
                                        center=mp.Vector3(n * metasurface_period - (N_periods_x-1)/2,
                                                          k * metasurface_period - (N_periods_y-1)/2,
                                                          0) + center + centroid,
                                        axis = mp.Vector3(z=1),
                                        material = medium_groove)
                else:
                    scatter = mp.Block(e1=mp.Vector3(cos(tilt), sin(tilt), 0),
                                       e2=mp.Vector3(cos(tilt+pi/2), sin(tilt+pi/2), 0),
                                       size=mp.Vector3(sc_length, sc_width, thickness),
                                       center=mp.Vector3( (n - (N_periods_x-1)/2) * metasurface_period,
                                                          (k - (N_periods_y-1)/2) * metasurface_period,
                                                         0) + center,
                                       material=medium_groove)

                metasurface.append(scatter)
                metasurface[-1].name  = f'Scatter_{n}_{k}'
                metasurface[-1].group = 'Metasurface'

    return metasurface

def pol_splitting_grating(  medium_groove=mp.Medium(epsilon=2),
                            D=2, metasurface_period=0.4, scatter_length = 0.4,
                            scatter_width=0.1, scatter_tilt=pi/3,
                            scatter_shape = '', scatter_disposition='radial',
                            topology='circular',
                            n_rings=9, n_arms=0, lambda_bsw=0.5,
                            thickness=1, center=mp.Vector3(0, 0, 0)):
    """
    Circular polarization sensitive grating. Similar tu the metasurface-like
    grating, but with specific purpose.
    """
    sc_width = scatter_width
    sc_length = scatter_length

    # following lists are for v-shaped scatters
    rel_width = sc_width / sc_length
    scatter_vertices_x = [0, .5, .5-rel_width/0.866, 0, -.5+rel_width*0.866, -.5]
    scatter_vertices_y = [-.5, .5, .5, -.5+rel_width/0.5, .5, .5]
    scatter_vertices = np.array([scatter_vertices_x, scatter_vertices_y]).transpose()

    metasurface = []

    if n_rings != 0:
        for n in range(n_rings):
            r0 = D/2 + n * metasurface_period

            # the following if statement is just for setting the actual
            # beginning of the scatter at D/2
            if scatter_shape == 'V' or scatter_shape == 'v':
                r0 += metasurface_period/2
            else:
                r0 += sc_width/2

            if scatter_disposition == 'filled':
                L = 2*pi*r0 + pi*metasurface_period*n_arms      # length of on turn in archimede's spiral
                N_scatters = round(L / metasurface_period)
            if scatter_disposition == 'radial':
                r1 = r0 - n * metasurface_period
                L = 2*pi*r1 + pi*metasurface_period*n_arms
                N_scatters = round(L / metasurface_period)

            if n_arms != 0:
                N_scatters -= np.mod(N_scatters,n_arms)     # make N_scatters divisible by n_arms
                theta = np.linspace(0, 2*pi/n_arms, int(N_scatters/n_arms))
                theta = np.tile(theta, n_arms)
            else:
                theta = np.zeros((N_scatters,1))

            for k in range(N_scatters):
                location_tilt = 2*pi / N_scatters * k

                if topology == 'spiral':
                    radius = r0 + (lambda_bsw * n_arms * theta[k]/2/pi)
                    tilt = pi/2 + location_tilt + n*scatter_tilt
                elif topology == 'circular':
                    radius = r0
                    tilt = pi/2 + location_tilt + n*scatter_tilt + theta[k]
                else:
                    raise ValueError('Topology can be either "spiral" or "circular"')

                if scatter_shape == 'V' or scatter_shape == 'v' :
                    tilt = tilt - pi
                    vertices = [mp.Vector3(v[0],v[1],0).rotate(mp.Z,tilt)*sc_length for v in scatter_vertices]
                    centroid =  sum(vertices, mp.Vector3(0)) * (1.0 / len(vertices))

                    scatter =  mp.Prism(vertices = vertices,
                                        height = thickness,
                                        center=mp.Vector3(radius * cos(location_tilt),
                                                          radius * sin(location_tilt),
                                                          0) + center + centroid,
                                        axis = mp.Vector3(z=1),
                                        material = medium_groove)
                else:
                    scatter = mp.Block(e1=mp.Vector3(cos(tilt), sin(tilt), 0),
                                       e2=mp.Vector3(cos(tilt+pi/2), sin(tilt+pi/2), 0),
                                       size=mp.Vector3(scatter_length*0.75, scatter_width, thickness),
                                       center=mp.Vector3(radius * cos(location_tilt),
                                                         radius * sin(location_tilt),
                                                         0) + center,
                                       material=medium_groove)

                metasurface.append(scatter)
                metasurface[-1].name  = f'Scatter_{n}_{k}'
                metasurface[-1].group = 'Metasurface'

    return metasurface

def spiral_grating(medium_groove=mp.Medium(epsilon=2),
                   D=0.4, d=None, DBR_period=0.2, FF=0.5, N_rings=10,
                   N_arms=2, thickness=0, center=mp.Vector3(0, 0, 0)):
    """
    Elliptic DBR cavity created as a sequence of concentric cylinders.

    Parameters
    ----------
    medium_groove : mp.Medium(),
        Medium of the rings. The default is mp.Medium(epsilon=2).
    D : Real Positive
        Major axis diametre of the centre of the cavity.
        The default is 0.4.
    d : Real Positive
        Minor axis diameter of the centre of the cavity.
        The default is 0.3.
    DBR_period : Real Positive
        Radial period of the DBR.
        The default is 0.2.
    FF : Real, optional
        Fill factor defined as a number between 0 and 1.
        The default is 0.5.
    N_rings : Integer
        Number of DBR rings.
        The default is 10.
    N_arms : Integer
        Number of arms of the spiaral. a.k.a. topological charge.
        Default is 2
    orientation : meeep Vector3()
        Orientation of the prism. Useful for 1D symulations.
        Default is along z-axis

    Returns
    -------
    List
        Returns a list of (meep) geometric objects.

    """
    device = []

    if d == None:
        d = D

    polygons = grating_veritices(DBR_period, D/2, d/2, N_rings, n_arms = N_arms,  FF = FF)

    for polygon in polygons:
        vertices = [mp.Vector3(v[0],v[1],0) for v in np.transpose(polygon)]
        centroid =  sum(vertices, mp.Vector3(0)) * (1.0 / len(vertices))

        c1 = mp.Prism(vertices = vertices,
                      height = thickness,
                      axis = mp.Vector3(0,0,1),
                      center = center + centroid,
                      material = medium_groove)
        device.append(c1)

    return device

def grating_veritices(period, start_radius1,
                      start_radius2=0, N_periods=10, n_arms=0,  FF=0.5,
                      spacer='empty') :
    """
    Function for generating the list of vertices por the circular, spiral
    and elliptic gratings.
    """

    if start_radius2 == 0:
        start_radius2 = start_radius1
    a = start_radius1
    b = start_radius2

    vert_list = []

    if n_arms != 0 :
        half_res = np.sum( [ (max([a,b]) + period*i) * 2*pi for i in range(N_periods)]) / period
        half_res = max(int(half_res), 32*int(N_periods/n_arms))
        half_res -= np.mod(half_res, int(N_periods/n_arms)) + 1

        res = 2 * half_res

        theta= np.linspace(0, 2*pi*N_periods/n_arms, half_res)

        # Each arm of the spiral is defined by a polygon
        for j in range(np.abs(n_arms)) :

            # from where to start the spiral, might even be elliptic
            if a == b == 0 :
                start_radius = np.zeros((half_res,))
            else:
                start_radius = a * b / np.sqrt(
                                          (b * np.cos(theta + 2*pi*j/n_arms) )**2 +
                                          (a * np.sin(theta + 2*pi*j/n_arms) )**2 )
            # parametrize the radius
            radius = start_radius+period*theta/(2*pi/n_arms)

            vertices = np.zeros((2, res))
            vertices[0, 0:half_res] = radius*np.cos(theta+2*pi*j/n_arms)
            vertices[1, 0:half_res] = radius*np.sin(theta+2*pi*j/n_arms)
            vertices[0, half_res:res] = np.flip( (radius+period*FF) * np.cos(theta+2*pi*j/n_arms) )
            vertices[1, half_res:res] = np.flip( (radius+period*FF) * np.sin(theta+2*pi*j/n_arms) )

            vert_list.append(vertices)

            # close the spiral on the origin if the initial radius is null
            if min(a,b) == 0 :
                centre_res = 15

                r_centre = np.linspace(period*FF, 0, centre_res)

                theta_centre = np.linspace(0, -pi/n_arms, centre_res)

                vert = np.zeros( (2,centre_res) )
                vert[0, :] = r_centre * np.cos(theta_centre + 2*pi*j/n_arms)
                vert[1, :] = r_centre * np.sin(theta_centre + 2*pi*j/n_arms)

                vert_list.append(vert)

    else :
        if spacer == 'full':
            # make a circle at the centre
            extra_r = period * (1 - FF)

            theta = np.linspace(0, 2*pi*(1-1/32), 32)
            radius = a * b / np.sqrt( (b*np.cos(theta))**2 +
                                      (a*np.sin(theta))**2 )

            vertices = np.zeros((2,32))
            vertices[0,:] = radius * np.cos(theta)
            vertices[1,:] = radius * np.sin(theta)

            vert_list.append(vertices)
        elif spacer == 'empty':
            extra_r = 0
            print('spacer is empty')
        else:
            raise ValueError("Invalid spacer value. Either 'empty' or 'full'")

        for j in range(N_periods):
            half_res = int( (max([a,b]) + period*j) * 2*pi / period)
            half_res = max(int(half_res), 32)
            half_res -= np.mod(half_res,4) +1

            res = 2 * half_res

            theta = np.linspace(0, 2*pi, half_res)

            # first circle radius, can be elliptic
            start_radius = a * b / np.sqrt( (b*np.cos(theta))**2 +
                                            (a*np.sin(theta))**2 )
            radius = start_radius + period*j + extra_r
            vertices = np.zeros((2,res))
            vertices[0,0:half_res] = radius * np.cos(theta)
            vertices[1,0:half_res] = radius * np.sin(theta)
            vertices[0,half_res:res] = np.flip( (radius+period*FF) * np.cos(theta))
            vertices[1,half_res:res] = np.flip( (radius+period*FF) * np.sin(theta))
            vertices[1,half_res] -= .001    # this is necessary for solving a bug
                                            # see https://github.com/NanoComp/libctl/issues/61
            vert_list.append(vertices)

    return vert_list

def create_openscad(sim, scale_factor=1):
    try:
        import openpyscad as ops
    except ModuleNotFoundError:
        print("WARNING openpyscad is not installed in this environment")
    else:
        try:
            scad_name = f"{sim.name}.scad"
        except:
            print("WARNING Simulation doesn't have a name")
            scad_name = "meep_sim.scad"

        with open(scad_name,"w") as file:
            file.write(f"//Simulation {scad_name}\n")

        full_scad = ops.Union()
        for obj in sim.geometry:
            centroid = mp.Vector3()
            if obj.__class__ == mp.Block:
                cube = ops.Cube([obj.size.x*scale_factor,
                                 obj.size.y*scale_factor,
                                 obj.size.z*scale_factor], center = True)
                tiltz = np.arctan2(obj.e1.y, obj.e1.x)
                # print(obj.e1)
                if tiltz != 0:
                    cube = cube.rotate([0, 0, tiltz/np.pi*180])
                # tiltx = np.arctan2(obj.e1.z, obj.e1.x)
                # if tiltx != 0:
                #     cube = cube.rotate([0, -tiltx/np.pi*180, 0])
                # tilty = np.arctan2(obj.e1.z, obj.e1.y)
                # if tilty != 0:
                #     cube = cube.rotate([tilty/np.pi*180, 0, 0])

                index = np.round(np.sqrt(obj.material.epsilon_diag.x),2)
                if index == 2.53:
                    color = [1, 0, 0]
                elif index == 1.65:
                    color = [0, 1, 0]
                elif index == 1.46:
                    color = [0, .6, .6]
                elif index == 1.48:
                    color = [1, 1, 0]
                elif index == 2.08:
                    color = [.5, 0, 0]
                    continue
                else:
                    color = [0, .3, .33]
                cube = cube.color(color)

                scad = cube

            elif obj.__class__ == mp.Cylinder:
                cyl = ops.Cylinder(h=obj.height*scale_factor,
                                    r=obj.radius*scale_factor, center = True)

                scad = cyl

            elif obj.__class__ == mp.Prism:
                vertices = [ [v.x*scale_factor, v.y*scale_factor ] for v in obj.vertices ]
                vertices[int(len(vertices)/2)][1] += .001 # undo bug correction in function "grating_veritices"
                centroid =  sum(obj.vertices, mp.Vector3(0)) * (scale_factor / len(vertices))
                base = ops.Polygon(vertices)
                scad = ops.Linear_Extrude(
                    height = obj.height*scale_factor,
                    center = True)
                scad.append(base)

            scad = scad.translate([obj.center.x*scale_factor - centroid.x,
                                   obj.center.y*scale_factor - centroid.y,
                                   obj.center.z*scale_factor])

            if full_scad == None:
                full_scad = scad
            else:
                full_scad.append(scad)

        full_scad.write(scad_name)#, mode='a')
        print("Openscad created")

        # sim_domain = ops.Cube([(sim.cell_size.x - sim.PML_width) * scale_factor,
        #                         (sim.cell_size.y - sim.PML_width) * scale_factor,
        #                         (sim.cell_size.z - sim.PML_width) * scale_factor], center = True)
        # sim_domain = sim_domain.color([.5, .5, .5, .5])
        # sim_domain = sim_domain.translate([sim.geometry_center.x*scale_factor,
        #                                     sim.geometry_center.y*scale_factor,
        #                                     sim.geometry_center.z*scale_factor])


text:
import os
import time

import numpy as np

import tensorflow as tf
import tensorflow.keras.backend as K
from tensorflow.compat.v1.train import Saver

def sanity_check(layers, layers_to_exploit):
    """
    Basic sanity check for layers and gradients to exploit based on model layers
    """
    if layers_to_exploit and len(layers_to_exploit):
        assert np.max(layers_to_exploit) <= len(layers),\
            "layer index greater than the last layer"

def time_taken(self, start_time, end_time):
    """
    Calculates difference between 2 times
    """
    delta = end_time - start_time
    hours = int(delta / 3600)
    delta -= hours * 3600
    minutes = int(delta / 60)
    delta -= minutes * 60
    seconds = delta
    return hours, minutes, np.int(seconds)

class attack_utils():
    """
    Utilities required for conducting membership inference attack
    """

    def __init__(self, directory_name='latest'):
        self.root_dir = os.path.abspath(os.path.join(
                                        os.path.dirname(__file__),
                                        "..", ".."))
        self.log_dir = os.path.join(self.root_dir, "logs")
        self.aprefix = os.path.join(self.log_dir,
                                    directory_name,
                                    "attack",
                                    "model_checkpoints")
        self.dataset_directory = os.path.join(self.root_dir, "datasets")

        if not os.path.exists(self.aprefix):
            os.makedirs(self.aprefix)
        if not os.path.exists(self.dataset_directory):
            os.makedirs(self.dataset_directory)

    def get_gradshape(self, variables, layerindex):
        """
        Returns the shape of gradient matrices
        Args:
        -----
        model: model to attack
        """
        g = (layerindex-1)*2
        gradshape = variables[g].shape
        return gradshape

    def get_gradient_norm(self, gradients):
        """
        Returns the norm of the gradients of loss value
        with respect to the parameters

        Args:
        -----
        gradients: Array of gradients of a batch
        """
        gradient_norms = []
        for gradient in gradients:
            summed_squares = [K.sum(K.square(g)) for g in gradient]
            norm = K.sqrt(sum(summed_squares))
            gradient_norms.append(norm)
        return gradient_norms

    def get_entropy(self, model, features, output_classes):
        """
        Calculates the prediction uncertainty
        """
        entropyarr = []
        for feature in features:
            feature = tf.reshape(feature, (1, len(feature.numpy())))
            predictions = model(feature)
            predictions = tf.nn.softmax(predictions)
            mterm = tf.reduce_sum(input_tensor=tf.multiply(predictions,
                                                           np.log(predictions)))
            entropy = (-1/np.log(output_classes)) * mterm
            entropyarr.append(entropy)
        return entropyarr

    def split(self, x):
        """
        Splits the array into number of elements equal to the
        size of the array. This is required for per example
        computation.
        """
        split_x = tf.split(x, len(x.numpy()))
        return split_x

    def get_savers(self, attackmodel):
        """
        Creates prefixes for storing classification and inference
        model
        """
        # Prefix for storing attack model checkpoints
        prefix = os.path.join(self.aprefix, "ckpt")
        # Saver for storing checkpoints
        attacksaver = Saver(attackmodel.variables)
        return prefix, attacksaver

    def createOHE(self, num_output_classes):
        """
        creates one hot encoding matrix of all the vectors
        in a given range of 0 to number of output classes.
        """
        return tf.one_hot(tf.range(0, num_output_classes),
                          num_output_classes,
                          dtype=tf.float32)

    def one_hot_encoding(self, labels, ohencoding):
        """
        Creates a one hot encoding of the labels used for
        inference model's sub neural network

        Args:
        ------
        zero_index: `True` implies labels start from 0
        """
        labels = tf.cast(labels, tf.int64).numpy()
        return tf.stack(list(map(lambda x: ohencoding[x], labels)))

    def intersection(self, to_remove, remove_from, batch_size):
        """
        Finds the intersection between `to_remove` and `remove_from`
        and removes this intersection from `remove_from`
        """
        to_remove = to_remove.unbatch()
        remove_from = remove_from.unbatch()

        m1, m2 = dict(), dict()
        for example in to_remove:
            hashval = hash(bytes(np.array(example)))
            m1[hashval] = example
        for example in remove_from:
            hashval = hash(bytes(np.array(example)))
            m2[hashval] = example

        # Removing the intersection
        extracted = {key: value for key,
                     value in m2.items() if key not in m1.keys()}
        dataset = extracted.values()
        features, labels = [], []
        for d in dataset:
            features.append(d[0])
            labels.append(d[1])
        finaldataset = tf.compat.v1.data.Dataset.from_tensor_slices(
            (features, labels))
        return finaldataset.batch(batch_size=batch_size)

ext: py

reason: import numpy


text:
#!/usr/bin/python
"""
Plot subsystems (x, vx), (y, vy), (z, vz) by listening to the PositionStateEstimator.

Authors: David Fridovich-Keil ( dfk@eecs.berkeley.edu )
"""

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

import rospy
from crazyflie_msgs.msg import PositionStateStamped

#matplotlib.rcParams['text.usetex'] = True

# Set up planner reference point in global scope.
planner_x = 0.0
planner_y = 0.0
planner_z = 0.0

received_reference = False

# Set up plots in global scope.
plt.ion()
fig, ax = plt.subplots(3)
plots = []

# X-subsystem.
plots.append(ax[0].scatter([], []))
ax[0].set_title("X subsystem")
ax[0].set_xlabel("V_{x} (m/s)")
ax[0].set_ylabel("x (m)")

#ax[0].set_xlabel(r"$V_{x}$ (m/s)")
#ax[0].set_ylabel(r"$x$ (m)")
ax[0].set_xlim(-2, 2)
ax[0].set_ylim(-2, 2)

# Y-subsystem.
plots.append(ax[1].scatter([], []))
ax[1].set_title("Y subsystem")
ax[1].set_xlabel("V_{y} (m/s)")
ax[1].set_ylabel("y (m)")

#ax[1].set_xlabel(r"$V_{y}$ (m/s)")
#ax[1].set_ylabel(r"$y$ (m)")
ax[1].set_xlim(-2, 2)
ax[1].set_ylim(-2, 2)

# Z-subsystem.
plots.append(ax[2].scatter([], []))
ax[2].set_title("Z subsystem")
ax[2].set_xlabel("V_{z} (m/s)")
ax[2].set_ylabel("z (m)")

#ax[2].set_xlabel(r"$V_{z}$ (m/s)")
#ax[2].set_ylabel(r"$z$ (m)")
ax[2].set_xlim(-2, 2)
ax[2].set_ylim(-2, 2)

# Callback to plot a new state in each of the three subsystem plots.
def StateCallback(msg):
    global planner_x, planner_y, planner_z
    global points, ax
    global received_reference

    if not received_reference:
        return

    MAX_POINTS = 5000

    # Unpack msg.
    x = msg.state.x - planner_x
    vx = msg.state.x_dot
    y = msg.state.y - planner_y
    vy = msg.state.y_dot
    z = msg.state.z - planner_z
    vz = msg.state.z_dot

    # X subsystem.
    points = plots[0].get_offsets()
    points = np.append(points, np.array([vx, x]))

    if len(points) > MAX_POINTS:
        points = points[-MAX_POINTS:]

    plots[0].set_offsets(points)
    ax[0].set_xlim(points[:, 0].min() - 0.1, points[:, 0].max() + 0.1)
    ax[0].set_ylim(points[:, 1].min() - 0.1, points[:, 1].max() + 0.1)

    # Y subsystem.
    points = plots[1].get_offsets()
    points = np.append(points, np.array([vy, y]))

    if len(points) > MAX_POINTS:
        points = points[-MAX_POINTS:]

    plots[1].set_offsets(points)
    ax[1].set_xlim(points[:, 0].min() - 0.1, points[:, 0].max() + 0.1)
    ax[1].set_ylim(points[:, 1].min() - 0.1, points[:, 1].max() + 0.1)

    # Z subsystem.
    points = plots[2].get_offsets()
    points = np.append(points, np.array([vz, z]))

    if len(points) > MAX_POINTS:
        points = points[-MAX_POINTS:]

    plots[2].set_offsets(points)
    ax[2].set_xlim(points[:, 0].min() - 0.1, points[:, 0].max() + 0.1)
    ax[2].set_ylim(points[:, 1].min() - 0.1, points[:, 1].max() + 0.1)

# Callback to update the planner state.
def ReferenceCallback(msg):
    global planner_x, planner_y, planner_z
    global received_reference

    planner_x = msg.state.x
    planner_y = msg.state.y
    planner_z = msg.state.z

    received_reference = True

# Main.
if __name__ == "__main__":
    rospy.init_node("state_plotter")
    rospy.Subscriber("/state/position", PositionStateStamped, StateCallback)
    rospy.Subscriber("/ref/planner", PositionStateStamped, ReferenceCallback)

    while not rospy.is_shutdown():
        # Update the figure.
        fig.canvas.draw()

        # Go to sleep.
        rospy.sleep(1.0)


text:
import numpy as np
import pandas as pd
import math

import ee
ee.Initialize()

import dirs

Fusion_Table_Field_1 = ee.FeatureCollection('ft:1rrcNOppcOV0na_6JHEmP8ycdQFP6ywHbheFV3wIi')
Fusion_Table_Field_2 = ee.FeatureCollection('ft:19rJUyf1HEqGFo1jMpkJ1-fsJV7WUlNOltVDMoLX1')
Fusion_Table_Field_3 = ee.FeatureCollection('ft:1R8vN0XlfE8wt22qfDgKNQ2nOnxQdp3bNifl01TgZ')
Fusion_Table_Field_4 = ee.FeatureCollection('ft:1gT2Ga-d3P8cKo0oUV9NPo25uFETfj9elqb_d9Tyr')
Fusion_Table_Field_5 = ee.FeatureCollection('ft:1i-qZRKlRPavbNPRIBCGLlHWqEr8uO5ninOudb_Ng')
Fusion_Table_Field_6 = ee.FeatureCollection('ft:1aLFSaWGfWonPwBcZ4IJPq-475zpzleS3Yh2zD09Q')
Fusion_Table_Field_7 = ee.FeatureCollection('ft:1pycN-8TvbR955LmbfETeAPGtQx9eJhRF5Dfaz2OM')
Fusion_Table_Field_8 = ee.FeatureCollection('ft:1uKwBKKhJtgWPfv3EtYssRNGWdUD6Rs5NZI5wuU63')
Fusion_Table_Field_9 = ee.FeatureCollection('ft:1zGUGjN3b_SVAyBil8styHJPfxf8jYEwz2tRWKJEd')
Fusion_Table_Field_10 = ee.FeatureCollection('ft:1zcgTzJ7uOo6LVxdI4pRTZRalKYMm3aFOVBOPB-UB')
Fusion_Table_Field_11 = ee.FeatureCollection('ft:1Z0ABnpqOqdLaXKLvdZEfuObaqMCfG98_YJprpM4h')
Fusion_Table_Field_12 = ee.FeatureCollection('ft:1hhASXotpHuL6Jw1YUoF-yVT81h035oLm-UNA2VbJ')
Fusion_Table_Field_13 = ee.FeatureCollection('ft:1tvvmMplT-rqFy2qr27g2aQp3rNA5QB7UA-OVoZkG')
Fusion_Table_Field_14 = ee.FeatureCollection('ft:1If-KPk4E8aHBht3CFLxgG_R0tZ4SUg6xLldWKZ6k')
Fusion_Table_Field_15 = ee.FeatureCollection('ft:1_un2eXsAg23SQNUlKeEdweIgw2igd5lxhW95sfUH')
Fusion_Table_Field_16 = ee.FeatureCollection('ft:1eXW_SExNomJ9yWfkYUrHDTo0pzZw5J6RApeZHdvx')
Fusion_Table_Field_17 = ee.FeatureCollection('ft:1zdnJlHfba4xoBa6E6hnywBIGz4Q4l5DI90DMK65A')
Fusion_Table_Field_18 = ee.FeatureCollection('ft:1DQih_EnaIq6MyTYxV5HTpAhRm3zPegRTRsCUVRln')
Fusion_Table_Field_19 = ee.FeatureCollection('ft:1XAedrGrZaLJYq15HPcWivtYlMmsg9_Sa3x-9zgXe')
Fusion_Table_Field_20 = ee.FeatureCollection('ft:18--k354WJBBQeeR9WmZNW330LC7F6hvNcFSOKojo')
Fusion_Table_Field_21 = ee.FeatureCollection('ft:1w3vVi440gNuOJMdlHJ3f_-ik3WUmzrChOFCLnje2')
Fusion_Table_Field_22 = ee.FeatureCollection('ft:1VcjBYp1ecNRXr0RXRU8nwqDM5Ja_hh8mUKLBi_xW')
crs_7 = 'EPSG:32610'
crs_8 = 'EPSG:32610'
crs_9 = 'EPSG:32610'
crs_10 = 'EPSG:32610'
crs_11 = 'EPSG:32610'
crs_12 = 'EPSG:32610'
crs_13 = 'EPSG:32610'
crs_14 = 'EPSG:32610'
crs_15 = 'EPSG:32610'
crs_16 = 'EPSG:32610'
crs_19 = 'EPSG:32610'
crs_20 = 'EPSG:32610'
crsTransform_7 = [30, 0, 596385, 0, -30, 4257315]
crsTransform_8 = [30, 0, 596385, 0, -30, 4257315]
crsTransform_9 = [30, 0, 596385, 0, -30, 4257315]
crsTransform_10 = [30, 0, 596385, 0, -30, 4257315]
crsTransform_11 = [30, 0, 596385, 0, -30, 4257315]
crsTransform_12 = [30, 0, 596385, 0, -30, 4257315]
crsTransform_13 = [30, 0, 596385, 0, -30, 4257315]
crsTransform_14 = [30, 0, 596385, 0, -30, 4257315]
crsTransform_15 = [30, 0, 596385, 0, -30, 4257315]
crsTransform_16 = [30, 0, 596385, 0, -30, 4257315]
crsTransform_19 = [30, 0, 596385, 0, -30, 4257315]
crsTransform_20 = [30, 0, 596385, 0, -30, 4257315]
crs_1 = 'EPSG:32611'
crs_2 = 'EPSG:32611'
crs_3 = 'EPSG:32611'
crs_4 = 'EPSG:32611'
crs_5 = 'EPSG:32611'
crs_6 = 'EPSG:32611'
crs_17 = 'EPSG:32611'
crs_18 = 'EPSG:32611'
crs_21 = 'EPSG:32611'
crs_22 = 'EPSG:32611'
crsTransform_1 = [30, 0, 155985, 0, -30, 4101915]
crsTransform_2 = [30, 0, 155985, 0, -30, 4101915]
crsTransform_3 = [30, 0, 155985, 0, -30, 4101915]
crsTransform_4 = [30, 0, 155985, 0, -30, 4101915]
crsTransform_5 = [30, 0, 155985, 0, -30, 4101915]
crsTransform_6 = [30, 0, 155985, 0, -30, 4101915]
crsTransform_17 = [30, 0, 155985, 0, -30, 4101915]
crsTransform_18 = [30, 0, 155985, 0, -30, 4101915]
crsTransform_21 = [30, 0, 155985, 0, -30, 4101915]
crsTransform_22 = [30, 0, 155985, 0, -30, 4101915]

def export_image_from_field(field_index, description = 'CMU_Research', description_temp = 'AVE Temp', start_date = '2007-01-01', end_date = '2013-12-31'):
    extrema_i = pd.read_csv(dirs.extreme_degrees)['extrema_field_' + str(field_index)]
    table = eval('ee_analysis.Fusion_Table_Field_' + str(field_index))
    out_name = 'FieldA'+str(field_index)
    out_name_temp='mean_temperature_field'+str(field_index)
    crs=eval('ee_analysis.crs_' + str(field_index))
    crsTransform=eval('ee_analysis.crsTransform_' + str(field_index))

    rect, rect_coordinates = create_rectangle(extrema_i)
    image = get_image_7years(table, rect, start_date, end_date)
    print(rect.getInfo())
    task_1 = export_image(image, description, out_name, rect_coordinates, crs, crsTransform)
    task_2 = export_temp_csv(start_date, end_date, description_temp, out_name_temp, rect)
    return task_1, task_2

def create_rectangle(extrema_i):#rect of same size as the one from ee code in gee
# check
    def delta_lat(meters):
        magic_number = 111111.0#ee.Number(math.pi).multiply(6371000.785).divide(180)  # transformation number
        dl = ee.Number(meters).divide(magic_number)
        d2 = ee.Number(meters).divide(magic_number)
        return dl,d2

# check (working with min max need two)
    def delta_lon(meters, min_lat, max_lat):
        magic_number = 111111.0#ee.Number(math.pi).multiply(6371000.785).divide(180)
        min_lat_rad = ee.Number(min_lat).divide(180).multiply(math.pi)
        max_lat_rad = ee.Number(max_lat).divide(180).multiply(math.pi)
        d_min = ee.Number(meters).multiply(ee.Number(1)).divide(min_lat_rad.cos()).divide(magic_number)
        d_max = ee.Number(meters).multiply(ee.Number(1)).divide(max_lat_rad.cos()).divide(magic_number)
        return d_min, d_max

    buffer_size_meters = 600
    d_lat_up, d_lat_down = delta_lat(buffer_size_meters)
    d_min_lon, d_max_lon = delta_lon(buffer_size_meters, extrema_i[1], extrema_i[3])
    lower_lat = ee.Number(extrema_i[1]).subtract(d_lat_down)
    lower_lon = ee.Number(extrema_i[0]).subtract(d_min_lon)  # how to call
    upper_lat = ee.Number(extrema_i[3]).add(d_lat_up)
    upper_lon = ee.Number(extrema_i[2]).add(d_max_lon)
    rect = ee.Geometry.Rectangle(lower_lon, lower_lat, upper_lon, upper_lat)
    rect_coordinates = rect.coordinates()

    return rect, rect_coordinates

def get_image_7years(table, rect, start_date, end_date):

    def get_image_collections():
        landsat_7_collection = ee.ImageCollection('LANDSAT/LE7_SR')#check
        crop_layer_collection = ee.ImageCollection(['USDA/NASS/CDL/2007b','USDA/NASS/CDL/2008','USDA/NASS/CDL/2009','USDA/NASS/CDL/2010','USDA/NASS/CDL/2011','USDA/NASS/CDL/2012','USDA/NASS/CDL/2013'])#check
        national_elevation = ee.Image('USGS/NED')#check
        return landsat_7_collection, crop_layer_collection, national_elevation

    def get_precip_monthly_average(start_date, end_date):

        def set_month_year_int(image):
            y = image.date().get('year').int16()
            m = image.date().get('month').int16()
            return image.set({'month': m, 'year': y})

        def get_monthly_mean(y):
            months = ee.List.sequence(1, 12)
            def go_over_months(m):
                monthly_image = ee.Image(precipitation_collection_daily.filterMetadata('year', 'equals', y)
                                         .filterMetadata('month', 'equals', m)
                                         .select('precipitation').sum()########## SUM() ?!!!!!#######
                                         .set('year', y)
                                         .set('month', m)
                                         .set('date', ee.Date.fromYMD(y, m, 1)))
                return monthly_image

            list_of_monthly_images_one_year = ee.List(months.map(go_over_months))
            return list_of_monthly_images_one_year

        precipitation_collection_daily = ee.ImageCollection('UCSB-CHG/CHIRPS/PENTAD').filterDate(start_date,end_date)
        precipitation_collection_daily = ee.ImageCollection(precipitation_collection_daily.map(set_month_year_int))

        years = ee.List.sequence(ee.Date(start_date).get('year'), ee.Date(end_date).get('year'))
        precipitation_collection = ee.ImageCollection(years.map(get_monthly_mean).flatten())
        return precipitation_collection

    def stack_ls7(image_1, image_2):
        millis_to_day = 1000 * 60 * 60 * 24
        date_code = image_1.date().millis().long().divide(millis_to_day).int16()#days
        year = image_1.date().get('year').int16()
        return ee.Image(image_2).addBands(ee.Image(image_1)
            .select("B1","B2","B3","B4","B5","B7")#
            .addBands(ee.Image(image_1).select("cfmask").toInt16()
            .addBands(ee.Image.constant(year)).toInt16()
            .addBands(ee.Image.constant(date_code)).toInt16()))

    def stack_crop_to_landsat(image_1, image_2):#check
        date_code = image_1.date().get('year').int16()
        helper = ee.String('cropland_')
        new_band = helper.cat(ee.Number.format(date_code))
        old_names = image_1.bandNames()
        new_names = old_names.set(-1, new_band)
        return ee.Image(image_2).addBands(ee.Image(image_1)
                                      .rename(new_names).select(new_band).int16())

    def stack_precip_to_landsat(image_1, image_2):#check
        year = ee.Date(image_1.get('date')).get('year').int16()
        month = ee.Date(image_1.get('date')).get('month').int16()
        helper = ee.String('precipitation_')
        new_band_1 = helper.cat(ee.Number.format(year))
        new_band_2 = new_band_1.cat(ee.Number.format(month))
        old_names = image_1.bandNames()
        new_names = old_names.set(-1, new_band_2)
        return ee.Image(image_2).addBands(ee.Image(image_1)
                                      .rename(new_names).select(new_band_2).int16())#maybe float()?

    landsat_7_collection, crop_layer_collection, national_elevation = get_image_collections()

    landsat_7_collection = landsat_7_collection.filterDate(start_date, end_date).filterBounds(rect)#left as it was
    out_image = ee.Image(landsat_7_collection.iterate(stack_ls7, ee.Image(-1000))).slice(1)

    crop_layer_collection = crop_layer_collection.filterBounds(rect)  # hopefully working keep in mind/worked in GEE
    out_image = ee.Image(crop_layer_collection.iterate(stack_crop_to_landsat, out_image))

    precipitation_collection = get_precip_monthly_average(start_date, end_date)
    out_image = ee.Image(precipitation_collection.iterate(stack_precip_to_landsat, out_image))

    out_image = out_image.addBands(national_elevation).int16()#left as it was
    out_image = out_image.addBands(ee.Terrain.aspect(national_elevation)).int16()
    out_image = out_image.addBands(ee.Terrain.slope(national_elevation)).int16()

    fkey = ee.Number(table.first().get('FieldID'))
    out_image = out_image.addBands(ee.Image.constant(fkey).int16())

    return out_image

def export_temp_csv(start_date, end_date, description_temp, out_name_temp, rect, scale=30):
    def get_temp_monthly_average(start_date, end_date):

        def temp_set_month_year_int(image):
            y = image.date().get('year').int16()
            m = image.date().get('month').int16()
            return image.set({'month': m, 'year': y})

        def temp_get_monthly_mean(y):
            months = ee.List.sequence(1, 12)
            def temp_go_over_months(m):
                monthly_image = ee.Image(temperature_collection_daily.filterMetadata('year', 'equals', y)
                                         .filterMetadata('month', 'equals', m)
                                         .select('air').mean()
                                         .set('year', y)
                                         .set('month', m)
                                         .set('date', ee.Date.fromYMD(y, m, 1)))
                return monthly_image

            list_of_monthly_images_one_year = ee.List(months.map(temp_go_over_months))
            return list_of_monthly_images_one_year

        temperature_collection_daily = ee.ImageCollection('NCEP_RE/surface_temp').filterDate(start_date,end_date)
        temperature_collection_daily = ee.ImageCollection(temperature_collection_daily.map(temp_set_month_year_int))

        years = ee.List.sequence(ee.Date(start_date).get('year'), ee.Date(end_date).get('year'))
        temperature_collection = ee.ImageCollection(years.map(temp_get_monthly_mean).flatten())
        return temperature_collection

    def temp_to_feature_collection(image):
        feature = ee.Feature(None, ee.Dictionary(image.reduceRegion(ee.Reducer.mean(),
                                            geometry = rect, scale = ee.Number(scale))
                    .set('year', image.get('year')))
                    .set('month', image.get('month')))
        return feature

    def export_table(feature_collection, description_temp, out_name_temp):  # left alone

        task = ee.batch.Export.table.toDrive(collection = feature_collection, folder = description_temp,
                                             description = out_name_temp, fileFormat = 'CSV')
        task.start()
        return task

    temperature_collection = get_temp_monthly_average(start_date, end_date)
    temp_feature_collection = ee.FeatureCollection(temperature_collection.map(temp_to_feature_collection))
    task = export_table(temp_feature_collection, description_temp, out_name_temp)
    return task

def export_image(image, description, out_name, rect_coordinates, crs, crsTransform):#left alone
   task = ee.batch.Export.image.toDrive(image = image, folder = description, description = out_name,
                                        crs=crs, crsTransform=crsTransform, region=rect_coordinates.getInfo()[0])
   task.start()
   return task

#running the code:
#extrema_1 = pd.read_csv(dirs.extreme_degrees)['extrema_field_' + str(7)]

#task = export_image_from_field(extrema_i = extrema_1, table = Fusion_Table_Field_7, description = 'CMU_Research', out_name = 'Field7_test',
##                                description_temp = 'average_temp', out_name_temp='mean_temperature',
#                                start_date = '2007-01-01', end_date = '2013-12-31',crs=crs_7, crsTransform=crsTransform_7)


text:
# classifies doodle drawings of $CATEGORIES.
# usage: call the the script with a trasnparent drawing as an argument
import os
import math
import sys
import keras
import numpy as np
from PIL import Image

model_path = os.path.join(os.getcwd(), "assets/doodle-model.h5")
model = keras.models.load_model(model_path)
TARGET_DIM = 28 # target image dimension for the neural network

CATEGORIES = [
    "Baseball",
    "Bowtie",
    "Clock",
    "Hand",
    "Hat",
    "Lightning",
    "Lollipop",
    "Mountain",
    "Pizza",
    "Star"
]

def load_image(path):
    # TODO load as black and white (0 or 1)
    img = Image.open(path) #Can be many different formats.
    img.load()
    print(f'loaded image of size {img.size}')
    print(f'image trasnparency status {img.getbbox()}')
    return img

def trim_image(img):
    img=img.crop(img.getbbox())
    # pad to square
    longer_side = max(img.size)
    horizontal_padding = (longer_side - img.size[0]) / 2
    vertical_padding = (longer_side - img.size[1]) / 2
    img = img.crop(
        (
            -math.ceil(horizontal_padding),
            -math.ceil(vertical_padding),
            img.size[0] + math.floor(horizontal_padding),
            img.size[1] + math.floor(vertical_padding)
        )
    )
    return img

# convert transparent pixels to white
def rasterize(img):
    canvas = Image.new('RGBA', img.size, (255,255,255,255)) # Empty canvas colour (r,g,b,a)
    canvas.paste(img, mask=img) # Paste the image onto the canvas, using it's alpha channel as mask
    return canvas

def resize_image(image):
    assert image.size[0] == image.size[1], f"image should be square {image.size}"
    wpercent = (TARGET_DIM/float(image.size[0]))
    hsize = int((float(image.size[1])*float(wpercent)))
    image = image.resize((TARGET_DIM,hsize), Image.ANTIALIAS)
    return image

def save_image(img, path):
    img.save(path)

# convert grayscale to black and white
def black_white(pixels):
    pixels[pixels < 255] = 0
    return pixels

def inverse(pixels):
    pixels *= -1
    pixels += 255
    return pixels

if __name__ == '__main__':
    if (len(sys.argv) < 2):
        print('pass in an image')
        sys.exit()
    path = sys.argv[1]
    image = load_image(path)
    image = trim_image(image)
    image = rasterize(image)
    image = resize_image(image)
    image = image.convert('L') # make it grayscale
    # save_image(image, 'cropped_drawing.png')

    image_data = np.asarray(image, np.float32)
    image_data = image_data.reshape(28*28) # redundant?
    image_data = inverse(image_data)

    image_data = image_data.reshape(1, 28, 28, 1)
    image_data /= 255 #
    print('image shape', image_data.shape)

    probs = model.predict(image_data)[0]
    predictions = dict(zip(CATEGORIES, probs))
    predictions = sorted(predictions.items(), key=lambda x:x[1], reverse=True)
    best_guess = predictions[0]
    print(predictions)
    print(f'{best_guess[0]},{best_guess[1]}')

max_issues_repo_name: brollb/sketch-adventures

max_forks_repo_name: brollb/sketch-adventures


text:
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import accuracy_score, confusion_matrix, classification_report

from osgeo import gdal

def readImg(img):
    # function to read .tif image files

    # Read heatmap
    image = gdal.Open(img)
    # Get band of heatmap, it's gray scale image!
    img_band = image.GetRasterBand(1)
    # Read the image as array
    image = img_band.ReadAsArray()
    # Normalize the pixel values in the range 0-1 acc. to max. normalization
    image = (image - image.min()) / (image.max() - image.min())
    return image.astype('float32')

def predict_eval(SavedModel, img_range, testPath, fn_input, fn_target, out_path, scale, nr = 1):
    # Predict and compare with ground truth

    image_arr = readImg(testPath + fn_input)
    image_tar = readImg(testPath + fn_target)

    if len(img_range) == 4: # If range was set
        xmin, xmax, ymin, ymax = img_range
        image_arr = image_arr[xmin:xmax, ymin:ymax]
        image_tar = image_tar[xmin:xmax, ymin:ymax]

    print("- " + str(scale) +"k", 'Example: ', image_arr.shape)
    evals_orig = evaluate(image_arr, image_tar)

    image_arr = rescaleImg(image_arr)
    image_tar = rescaleImg(image_tar)

    print("+ " + str(scale) +"k", 'Prediction: ', image_arr.shape)
    pred = model_predict(SavedModel, image_arr, num_runs = nr)
    evals_pred = evaluate(pred > 0.5, image_tar)

    save_prediction(pred, out_path, fn_input[:-4], '_' + str(scale) + '_' + str(nr) + '_out.png')
    return [evals_orig, evals_pred]

def predict_only(SavedModel, img_range, testPath, fn_input, out_path, scale, nr = 1):
    # Predict only

    image_arr = readImg(testPath + fn_input)

    if len(img_range) == 4: # If range was set
        xmin, xmax, ymin, ymax = img_range
        image_arr = image_arr[xmin:xmax, ymin:ymax]

    image_arr = rescaleImg(image_arr)

    print("+ " + str(scale) +"k", 'Prediction: ', image_arr.shape)
    pred = model_predict(SavedModel, image_arr, num_runs = nr)

    save_prediction(pred, out_path, fn_input[:-4], '_' + str(scale) + '_' + str(nr) + '_out.png')

def model_predict(newModel, input_image, num_runs):
    # Predict with model for n times
    m,n = input_image.shape
    for i in range(num_runs):
        input_image = np.reshape(input_image, (1, m, n, 1))
        input_image = newModel.predict([input_image])
        input_image = np.reshape(input_image,(m, n)) > 0.5
    return input_image

def rescaleImg(image_arr):
    # Rescale the image to 8n x 8n, because 3x downsampling

    if image_arr.shape[0] % 8 != 0:
        n = image_arr.shape[0] % 8
        new_x = image_arr.shape[0] - n
    else:
        new_x = image_arr.shape[0]

    if image_arr.shape[1] % 8 != 0:
        n = image_arr.shape[1] % 8
        new_y = image_arr.shape[1] - n
    else:
        new_y = image_arr.shape[1]

    image_arr = image_arr[:new_x, :new_y]

    return image_arr

def IoUcheck(img_input, img_output):
    # Pixelwise IoU score

    logic_and = np.sum(np.logical_and(img_output, img_input))
    logic_or = np.sum(np.logical_or(img_output, img_input))
    return logic_and/logic_or

def evaluate(image_arrA, image_arrB):
    # Calculate the evaluation metrics

    target_names = ['0', '1']

    y_true = image_arrB.flatten().astype(bool)
    y_pred = image_arrA.flatten().astype(bool)

    Accuracy = accuracy_score(y_true, y_pred)
    IoU = IoUcheck(y_true, y_pred)
    conf = confusion_matrix(y_true, y_pred)
    report = classification_report(y_true, y_pred, target_names=target_names, digits=8)

    print('Acc & IoU:', Accuracy, IoU)
    print('Error:', 1 - Accuracy)
    print(conf)
    print(report)

    return Accuracy, IoU, conf, report

def save_prediction(img, out_path, out_filename, subfix):

    fig = plt.figure(figsize=(img.shape[1] / 1000, img.shape[0] / 1000), dpi=100, frameon=False)
    ax = plt.Axes(fig, [0., 0., 1., 1.])
    ax.set_axis_off()
    fig.add_axes(ax)

    plt.imshow(img, cmap='gray')
    fig.savefig(out_path + out_filename + subfix, dpi=1000)

max_stars_repo_name: yuzzfeng/MapGen

max_issues_repo_name: yuzzfeng/MapGen

max_forks_repo_name: yuzzfeng/MapGen


text:
# Copyright (c) 2015, Andrew Delong and Babak Alipanahi All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this
# list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation and/or
# other materials provided with the distribution.
#
# 3. Neither the name of the copyright holder nor the names of its contributors
# may be used to endorse or promote products derived from this software without
# specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
# ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# Author's note:
#     This file was distributed as part of the Nature Biotechnology
#     supplementary software release for DeepBind. Users of DeepBind
#     are encouraged to instead use the latest source code and binaries
#     for scoring sequences at
#        http://tools.genes.toronto.edu/deepbind/
#
import sys
import re
import os
import os.path
import argparse
import numpy as np
import numpy.random as npr
import gzip
import csv
import cPickle
import deepbind_util as util
import deepity
import scipy
import shutil

# Warning: The multiprocessing module may be used indirectly, so do not put any
# unintentional statements outside of main().

def main():
    util.enable_reversecomplement()

    args   = loadargs()
    models = loadmodels(args)
    tfgroups = load_tfgroups(args)
    util.globals.flags.push("normalize_targets", True)

    for tfgroup in tfgroups:
        trdata = None
        if len(tfgroup["ids"]) == 0:
            print "No TFs to train on microarray %s"%tfgroup["train_fold"]
            continue

        if "calib" in args.steps:
            trdata = load_pbmdata(trdata, tfgroup["ids"], tfgroup["train_fold"], args, remove_probe_bias=True)
            util.calibrate(models, trdata, args.calibdir, nfold=args.nfold, ncalib=args.ncalib, allfolds=True)

        if "train" in args.steps:
            trdata = load_pbmdata(trdata, tfgroup["ids"], tfgroup["train_fold"], args, remove_probe_bias=True)
            util.train(models, trdata, args.calibdir, args.finaldir, nfold=1,          ntrial=args.ntrial, metric_key="pearson.r")

    for tfgroup in tfgroups:
        tedata = None

        newids = []
        for id in tfgroup["ids"]:
            if os.path.exists(args.outdir+"/final/"+id+"/model.pkl"):
                newids.append(id)
            else:
                print "WARNING: did not find model for %s, skipping" % id
        tfgroup["ids"] = newids

        if len(tfgroup["ids"]) == 0:
            print "No TFs to test on microarray %s"%tfgroup["train_fold"]
            continue

        if "test" in args.steps:
            tedata = load_pbmdata(tedata, tfgroup["ids"], tfgroup["test_fold"], args, remove_probe_bias=False)
            save_test_performance(tedata, tfgroup["ids"], tfgroup["test_fold"], args)

        if "report" in args.steps:
            tedata = load_pbmdata(tedata, tfgroup["ids"], tfgroup["test_fold"], args, remove_probe_bias=False)
            util.save_featuremaps(tedata, args.finaldir, args.reportdir)

    if "report" in args.steps:
        all_tfids = sum([tfgroup["ids"] for tfgroup in tfgroups],[])
        save_report(args.finaldir, args.reportdir, all_tfids, index_metric="pearson")
        save_pbm_performance_table(args, all_tfids)

    if "chip" in args.steps:
        all_tfids = sum([tfgroup["ids"] for tfgroup in tfgroups],[])
        save_chip_performance_table(args, all_tfids)

#########################################################################

def load_tfgroups(args):
    with open("../data/dream5/pbm/tfids.txt") as f:
        tfids = sorted([line.rstrip("\r\n") for line in f.readlines()])

    # Narrow down the list of ids based on which chunk we've been asked to compute
    tfids  = util.getchunktargets(args, tfids)

    # Group the TF ids into those that are trained on A versus those that are trained on B
    tfgroups = [ {"ids" : list(set(tfids).intersection(["C_%d"%(i+1) for i in range(20)]
                                                      + ["TF_%d"%(i+1) for i in range( 0,33)])), "train_fold" : "A", "test_fold" : "B" },
                 {"ids" : list(set(tfids).intersection( ["TF_%d"%(i+1) for i in range(33,66)])), "train_fold" : "B", "test_fold" : "A" }]

    return tfgroups

#################################

def loadargs():
    args = argparse.ArgumentParser(description="Generate the DREAM5 PBM and/or CHIP experiments.")
    args = util.parseargs("dream5", args)
    return args

#################################

def loadmodels(args):
    models = util.loadmodels(args, "cfg/regression/maxpool")
    for cfg in models.itervalues():
        cfg["model"].conv_seq[0].fsize = 24  # Override default max motif length
    return models

#################################

def load_probe_biases():

    if not os.path.exists("../data/dream5/pbm/probe_biases.npz"):

        # First find which rows belong to which microarray (A or B)
        with gzip.open("../data/dream5/pbm/sequences.tsv.gz") as f:
            f.readline()
            foldids = [line[0] for line in f]
            rowidx  = { "A" : [i for i in range(len(foldids)) if foldids[i] == "A"],
                        "B" : [i for i in range(len(foldids)) if foldids[i] == "B"] }

        # Then load the targets and split the rows into rows for A and rows for B
        with gzip.open("../data/dream5/pbm/targets.tsv.gz") as f:
            reader = csv.reader(f, delimiter='\t')
            colnames = reader.next()
            targets = np.array([[float(x) for x in row] for row in reader])

        biases = {}
        for foldid in ("A","B"):
            # For this microarray (A or B), calculate a multiplicative bias
            # for each one of its individual probes (rows).
            microarray_measurements = targets[rowidx[foldid]]
            bias = []
            for i in range(len(microarray_measurements)):
                probe_measurements = microarray_measurements[i,:].ravel()
                if np.all(np.isnan(probe_measurements)):
                    bias.append(np.nan)
                else:
                    bias.append(np.median(probe_measurements[~np.isnan(probe_measurements)]))
            biases[foldid] = np.array(bias)

        np.savez("../data/dream5/pbm/probe_biases.npz", **biases)

    with np.load("../data/dream5/pbm/probe_biases.npz") as f:
        biases = { key : f[key] for key in f.keys() }
    return biases

#################################

def load_pbmdata(trdata, tfids, fold, args, remove_probe_bias=False):
    if trdata is not None:
        return trdata  # Training data was already loaded in earlier step.
    maxrows = 10000 if args.quick else None

    # Determine which targetnames we're responsible for
    targetnames = gzip.open("../data/dream5/pbm/targets.tsv.gz").readline().rstrip().split("\t")
    targetcols = [i for i in range(len(targetnames)) if targetnames[i] in tfids]

    data = util.datasource.fromtxt("../data/dream5/pbm/sequences.tsv.gz", None, "../data/dream5/pbm/targets.tsv.gz", targetcols=targetcols, foldfilter=fold,  maxrows=maxrows)

    if remove_probe_bias:
        # Remove per-probe multiplicative bias
        biases = load_probe_biases()[fold].reshape((-1,1))
        data.targets /= biases
        data.Y /= biases

    return data

#################################

def save_test_performance(data, tfids, fold, args):

    print "Predicting...", tfids
    # Generate predictions for each model
    predictions = util.predict(data, args.finaldir, args.reportdir)
    print "done"

    for tfname in predictions:
        print tfname
        Z = predictions[tfname]

        # Restore the original (unpreprocessed) scale and normalization of the predictions
        # Then add the opposite microarray's estimated probe value to each probe
        with open(args.finaldir + "/" + tfname +"/preprocessors.pkl") as f:
            normalizer = cPickle.load(f)['targets'][0]
        Z *= normalizer.scales[0]
        Z += normalizer.biases[0]

        # Add per-probe multiplicative bias to predictions
        probe_biases = load_probe_biases()[fold].reshape((-1,1))
        if args.quick:
            probe_biases = probe_biases[:len(Z),:]
        M = ~np.isnan(probe_biases)
        Z[M] *= probe_biases[M]

        #with open("predictions/%s_%s.tsv" % (tfname, {"A":"B","B":"A"}[foldid]), 'w') as f:
        #    f.writelines([str(x)+"\n" for x in Z.ravel()])
        z = Z.ravel()
        y = data.Y[:,data.targetnames.index(tfname)].ravel()
        rowidx = data.rowidx.ravel()
        mask = ~np.isnan(y)
        zscore4 = np.std(y[mask])*4+np.mean(y[mask])
        util._update_metrics(args.finaldir, tfname, "test", rowidx, z, y, aucthresh=(zscore4,zscore4))

#################################

def save_pbm_performance_table(args, tfids):
    metric_keys = ["pearson.r","spearman.r","pearson.p","spearman.p","auc","auc.mean"]
    with open(args.reportdir+"/performance_pbm.txt", "w") as f:
        f.write("protein\t%s\n" % "\t".join(metric_keys))
        for tfid in ["TF_%d" % i for i in range(1,67)]:
            f.write(tfid)
            if tfid in tfids:
                metrics = util.load_metrics(args.finaldir+"/"+tfid+"/metrics.txt")
                for key in metric_keys:
                    f.write("\t%s" % metrics["test"].get(key,np.nan))
            f.write("\n")

    return

#################################

def save_chip_performance_table(args, tfids):
    with open(args.reportdir+"/performance_chipseq.txt", "w") as f:
        f.write("Background\tFile\tauc.mean\tauc.std\n")
        chipseqids = ["TF_%d"%i for i in [23,25,31,40,44]]
        for tfid in chipseqids:
            if tfid not in tfids:
                continue
            for windowsize in [51,100]:
                for background in ["dinuc", "full_genomic", "genomic"]:
                    seqfile = "../data/dream5/chipseq/%s_CHIP_%d_%s.seq" % (tfid, windowsize, background)
                    invivodata = util.datasource.fromtxt(seqfile+"[0]", None, seqfile+"[1]", sequencenames=["seq"], targetnames=[tfid])
                    invivodata.targetnames = [tfid]

                    predictions = util.predict(invivodata, args.finaldir, args.reportdir, include=[tfid])
                    z = predictions[tfid].ravel()
                    y = invivodata.Y.ravel()
                    metrics = util.calc_metrics(z, y)
                    s = "%s\t%s\t%.3f\t%.3f\n" % (background, seqfile, metrics['auc.mean'], metrics['auc.std'])
                    print s,
                    f.write(s)

#################################

def save_report(modeldir, reportdir, tfids, index_metric="auc"):

    util.makepath(reportdir)
    index_html = open(reportdir+"/index.html",'w')
    index_html.write("<title>Training report</title>\n")
    index_html.write("
\n")
    index_html.write("
Name	train %s	test %s
\n" % (index_metric, index_metric))

    for tfid in ["TF_%d" % i for i in range(1,67)]:
        if tfid not in tfids:
            continue
        print tfid
        util.makepath(reportdir+"/"+tfid)

        # Load PFMs, convert them to logo images, and dump them to the report directory
        with open(reportdir+"/%s.pfms.pkl"%tfid) as f:
            _ = cPickle.load(f)
            pfms = _["pfms"]
            ics  = _["ic"]
            counts = _["counts"]
        pfm_order = np.argsort(-ics)
        logos = []
        for j in pfm_order:
            pfm = pfms[j]
            ic  = ics[j]
            if ic <= 0:
                continue

            pfm_rev = np.fliplr(np.flipud(pfm))
            logo_fwd = deepity.tape2logo.tape2logo(pfm.T,     height=50, letterwidth=10, bufferzoom=4, vmax=1.0, style="seqlogo", rna=False)
            logo_rev = deepity.tape2logo.tape2logo(pfm_rev.T, height=50, letterwidth=10, bufferzoom=4, vmax=1.0, style="seqlogo", rna=False)
            logo_filename = "%s/%s/pfm%02d" % (reportdir, tfid, len(logos))
            scipy.misc.imsave(logo_filename+"_fwd.png", logo_fwd)
            scipy.misc.imsave(logo_filename+"_rev.png", logo_rev)
            logos.append((j, os.path.basename(logo_filename)+"_fwd.png", ic, counts[j]))

        # Load train/test metrics so we can print them in the HTML report
        metrics_file = "%s/%s/metrics.txt" % (modeldir, tfid)
        metrics = deepity.load_metrics(metrics_file)

        # Add row for this TF in main index.html
        index_html.write("
\n")
        index_html.write("%s	\n" % (tfid, tfid))
        if index_metric=="auc":
            index_html.write("%s        ± %s	\n" % (metrics["train"]["auc.mean"], metrics["train"]["auc.std"]))
            index_html.write("%s ± %s	\n" % (metrics["test"]["auc.mean"],  metrics["test"]["auc.std"]))
        elif index_metric=="pearson":
            index_html.write("%s        (p=%s)	\n" % (metrics["train"]["pearson.r"], metrics["train"]["pearson.p"]))
            index_html.write("%s (p=%s)	\n" % (metrics["test"]["pearson.r"],  metrics["test"]["pearson.p"]))
        index_html.write("
\n")

        # Build page showing filters and sequence logos for this specific TF
        tf_html = open(reportdir+"/"+tfid+"/index.html", 'w')
        tf_html.write("<title>Training report - %s</title>\n" % tfid)
        tf_html.write("
%s
\n" % tfid)
        # tf_html.write("gradient maps)
\n"%(tfid,tfid))
        tf_html.write("""

        """)

        with open(metrics_file) as f:
            metrics_text = f.read()
        tf_html.write("
%s
\n"% metrics_text)
        if os.path.exists(modeldir+"/%s/predict.scatter-test.png" % tfid):
            tf_html.write("

TEST predictions

\n" % tfid)
        if os.path.exists(modeldir+"/%s/predict.scatter-train.png" % tfid):
            tf_html.write("

TRAINING predictions

\n" % tfid)

        # Then generate a table for the complete model
        tfdir = modeldir+"/"+tfid
        tf_html.write("

Feature Logos
\n")
        tf_html.write("<input type=\"button\" value=\"toggle strand\" onclick=\"toggle_strand();\"/>
")
        tf_html.write("
\n")
        for filter_index, logo, ic, count in logos:
            tf_html.write("
%d		%.1f bits,	 %d activations
" % (filter_index, logo, ic, count))
        tf_html.write("

\n")

        # Now show the actual model.
        '''
        shutil.copy(tfdir_final[0]+"/fold0.report/filters_.conv_seq(0).color.png", basedir+"report/%s/final_filters.png" % (tfid))
        shutil.copy(tfdir_final[0]+"/fold0.report/filters_.conv_seq(0).logo.png",  basedir+"report/%s/final_logos.png" % (tfid))
        shutil.copy(tfdir_final[0]+"/fold0.report/filters_.conv_seq(1).png",       basedir+"report/%s/final_biases.png"  % (tfid))
        shutil.copy(tfdir_final[0]+"/fold0.report/filters_.combiner.png",          basedir+"report/%s/final_weights.png" % (tfid))
        os.system("convert %sreport/%s/final_filters.png -rotate -90 %sreport/%s/final_filters_rot.png" % (basedir,tfid,basedir,tfid))
        '''

        tf_html.write("

Actual DeepBind model
\n")
        if os.path.exists(tfdir + "/model.conv_seq(1).color.png"):
            tf_html.write("Filters:
\n")
            tf_html.write("
\n" % tfid)

        if os.path.exists(tfdir + "/model.combiner.png"):
            tf_html.write("
Combiner layer:
\n")
            tf_html.write("
\n" % tfid)

        tf_html.write("\n")

    index_html.write("
\n")
    index_html.write("\n")
    index_html.close()

if __name__=="__main__":
    #util.disable_multiprocessing()
    main()

max_stars_repo_name: gifford-lab/deepbind-docker

max_issues_repo_name: gifford-lab/deepbind-docker

max_forks_repo_name: gifford-lab/deepbind-docker

reason: import numpy,import scipy