text:
"""Visualizer for rllab-trained experiments."""

import argparse
import joblib
from matplotlib import pyplot as plt
import numpy as np
import os

from flow.core.util import emission_to_csv

from rllab.sampler.utils import rollout

def visualizer_rllab(args):
    # extract the flow environment
    data = joblib.load(args.file)
    policy = data['policy']
    env = data['env']

    # FIXME(ev, ak) only one of these should be needed
    # unwrapped_env = env._wrapped_env._wrapped_env.env.unwrapped
    # unwrapped_env = env.wrapped_env.env.env.unwrapped

    # if this doesn't work, try the one above it
    unwrapped_env = env._wrapped_env.env.unwrapped

    # Recreate experiment params
    tot_cars = unwrapped_env.vehicles.num_vehicles
    rl_cars = unwrapped_env.vehicles.num_rl_vehicles
    max_path_length = int(env.horizon)
    flat_obs = env._wrapped_env.observation_space.flat_dim
    obs_vars = unwrapped_env.obs_var_labels or []
    num_obs_var = flat_obs / tot_cars

    # Set sumo to make a video
    sumo_params = unwrapped_env.sumo_params
    sumo_params.emission_path = './test_time_rollout/'
    if args.no_render:
        sumo_params.render = False
    else:
        sumo_params.render = True
    unwrapped_env.restart_sumo(
        sumo_params=sumo_params, render=sumo_params.render)

    # Load data into arrays
    all_obs = np.zeros((args.num_rollouts, max_path_length, flat_obs))
    all_rewards = np.zeros((args.num_rollouts, max_path_length))
    rew = []
    for j in range(args.num_rollouts):
        # run a single rollout of the experiment
        path = rollout(env=env, agent=policy)

        # collect the observations and rewards from the rollout
        new_obs = path['observations']
        all_obs[j, :new_obs.shape[0], :new_obs.shape[1]] = new_obs
        new_rewards = path['rewards']
        all_rewards[j, :len(new_rewards)] = new_rewards

        # print the cumulative reward of the most recent rollout
        print('Round {}, return: {}'.format(j, sum(new_rewards)))
        rew.append(sum(new_rewards))

    # print the average cumulative reward across rollouts
    print('Average, std return: {}, {}'.format(np.mean(rew), np.std(rew)))

    # ensure that a reward_plots folder exists in the directory, and if not,
    # create one
    if not os.path.exists('plots') and not os.environ.get('TEST_FLAG', 1):
        os.makedirs('plots')

    # create an array of time
    sim_step = unwrapped_env.sumo_params.sim_step
    t = np.arange(max_path_length) * sim_step

    for obs_var_idx in range(int(num_obs_var)):
        if len(obs_vars) < obs_var_idx + 1:
            obs_var = 'Observation {0}'.format(obs_var_idx)
        else:
            obs_var = obs_vars[obs_var_idx]

        # plot mean value for observation for each vehicle across rollouts
        plt.figure()
        for car in range(tot_cars):
            center = np.mean(
                all_obs[:, :, tot_cars * obs_var_idx + car], axis=0)
            plt.plot(
                range(max_path_length),
                center,
                lw=2.0,
                label='Veh {}'.format(car))
        plt.ylabel(obs_var, fontsize=15)
        plt.xlabel('time (s)', fontsize=15)
        plt.title(
            '{2}, Autonomous Penetration: {0}/{1}'.format(
                rl_cars, tot_cars, obs_var),
            fontsize=16)
        plt.legend(loc=0)

        # save the plot in the "plots" directory unless we're testing
        if not os.environ.get('TEST_FLAG', 1):
            plt.savefig(
                'plots/{0}_{1}.png'.format(args.plotname, obs_var),
                bbox='tight')

        # plot mean values for the observations across all vehicles and all
        # rollouts
        car_mean = np.mean(
            np.mean(
                all_obs[:, :, tot_cars * obs_var_idx:tot_cars *
                        (obs_var_idx + 1)], axis=0), axis=1)
        plt.figure()
        plt.plot(t, car_mean)
        plt.ylabel(obs_var, fontsize=15)
        plt.xlabel('time (s)', fontsize=15)
        plt.title(
            'Mean {2}, Autonomous Penetration: {0}/{1}'.format(
                rl_cars, tot_cars, obs_var),
            fontsize=16)

        # save the plot in the "plots" directory
        if not os.environ.get('TEST_FLAG', 1):
            plt.savefig(
                'plots/{0}_{1}_mean.png'.format(args.plotname, obs_var),
                bbox='tight')

    # Make a figure for the mean rewards over the course of the rollout
    mean_reward = np.mean(all_rewards, axis=0)

    plt.figure()
    plt.plot(t, mean_reward, lw=2.0)
    plt.ylabel('reward', fontsize=15)
    plt.xlabel('time (s)', fontsize=15)
    plt.title(
        'Reward, Autonomous Penetration: {0}/{1}'.format(rl_cars, tot_cars),
        fontsize=16)

    # save the rewards plot in the "reward_plots" directory
    if not os.environ.get('TEST_FLAG', 1):
        plt.savefig('plots/{0}_reward.png'.format(args.plotname), bbox='tight')

    # if prompted, convert the emission file into a csv file
    if args.emission_to_csv:
        dir_path = os.path.dirname(os.path.realpath(__file__))
        emission_filename = '{0}-emission.xml'.format(
            unwrapped_env.scenario.name)

        emission_path = \
            '{0}/test_time_rollout/{1}'.format(dir_path, emission_filename)

        emission_to_csv(emission_path)

def create_parser():
    parser = argparse.ArgumentParser()
    parser.add_argument('file', type=str, help='path to the snapshot file')
    parser.add_argument(
        '--num_rollouts',
        type=int,
        default=100,
        help='Number of rollouts we will average over')
    parser.add_argument(
        '--no_render',
        action='store_true',
        help='Whether to render the result')
    parser.add_argument(
        '--plotname',
        type=str,
        default='traffic_plot',
        help='Prefix for all generated plots')
    parser.add_argument(
        '--emission_to_csv',
        action='store_true',
        help='Specifies whether to convert the emission file '
             'created by sumo into a csv file')
    return parser

if __name__ == '__main__':

    parser = create_parser()
    args = parser.parse_args()
    visualizer_rllab(args)

max_stars_repo_name: nskh/flow

max_issues_repo_name: nskh/flow

max_forks_repo_name: nskh/flow

reason: import numpy


text:
import numpy as np

def dist3D_segment_to_segment(p0,p1,p2,p3): # S1=(p0,p1), S2=(p2,3)
    # c++ from http://geomalgorithms.com/a07-_distance.html#dist3D_Segment_to_Segment%28%29
    SMALL_NUM =  0.00000001
    dot = lambda u,v: np.dot(u,v)
    norm = lambda v: np.sqrt(dot(v,v))
    d = lambda u,v: norm(u-v)

    u = p1-p0 # Vector   u = S1.P1 - S1.P0
    v = p3-p2 # Vector   v = S2.P1 - S2.P0
    w = p0-p2 # Vector   w = S1.P0 - S2.P0
    a = dot(u,u)         # always >= 0
    b = dot(u,v)
    c = dot(v,v)         # always >= 0
    d = dot(u,w)
    e = dot(v,w)
    D = a*c - b*b        # always >= 0
    sc, sN, sD = D,D,D       # sc = sN / sD, default sD = D >= 0
    tc, tN, tD = D,D,D       # tc = tN / tD, default tD = D >= 0

    # compute the line parameters of the two closest points
    if (D < SMALL_NUM):  # the lines are almost parallel
        sN = 0.0         # force using point P0 on segment S1
        sD = 1.0         # to prevent possible division by 0.0 later
        tN = e
        tD = c

    else:                 # get the closest points on the infinite lines
        sN = (b*e - c*d)
        tN = (a*e - b*d)
        if (sN < 0.0):        # sc < 0 => the s=0 edge is visible
            sN = 0.0
            tN = e
            tD = c

        elif (sN > sD):  # sc > 1  => the s=1 edge is visible
            sN = sD
            tN = e + b
            tD = c

    if (tN < 0.0):            # tc < 0 => the t=0 edge is visible
        tN = 0.0
        # recompute sc for this edge
        if (-d < 0.0):
            sN = 0.0
        elif (-d > a):
            sN = sD
        else:
            sN = -d
            sD = a

    elif (tN > tD):      # tc > 1  => the t=1 edge is visible
        tN = tD
        # recompute sc for this edge
        if ((-d + b) < 0.0):
            sN = 0
        elif ((-d + b) > a):
            sN = sD
        else:
            sN = (-d +  b)
            sD = a

    # finally do the division to get sc and tc
    sc,tc = 0.0,0.0
    # sc = (abs(sN) < SMALL_NUM ? 0.0 : sN / sD)
    if np.abs(sN) < SMALL_NUM:
        sc = 0.0
    else:
        sc = sN / sD

    #tc = (abs(tN) < SMALL_NUM ? 0.0 : tN / tD)
    if np.abs(tN) < SMALL_NUM:
        tc = 0.0
    else:
        tc = tN/tD

    # get the difference of the two closest points
    dP = w + (sc * u) - (tc * v)  # =  S1(sc) - S2(tc)

    return norm(dP)   # return the closest distance

ext: py


text:
"""
Glue outputs into a single binary file

Jinbo Wang
<jinbow@gmail.com>
Scripps Institution of Oceanography
August 26, 2015 """

import glob,os,sys
import numpy as np

#change this to the path of the data output
folder='/home/jwang0/temp_project/offline_particle/sose100/DP/run/'

#you can glue one case or multi-cases at the same time, 'DIMES_0001' etc. are casenames
casenames=['DIMES_0001','DIMES_0002','DIMES_0003']
casenames=['DIMES_0006']

#specify particle numbers and case numbers (the NPP value in the namelist) here:
npts=5000
npp=73

#you don't need to change the following.
varn=['XYZ','MLD','TSG','GRAD']
nn = [3,1,4,4]
for cn in casenames:
    for i in range(1,npp+1):
        casename='%s_%04i'%(cn,i)
        for iv,var in enumerate(varn):
            fns=sorted(glob.glob(folder+casename+'.%s.*.data'%var))
            print casename,var,"total %i files"%len(fns)
            f=open('filenames_%s'%casename,'w')
            n=nn[iv]
            d=np.fromfile(fns[0],'>f4')
            print d.shape,n
            d=d.reshape(n,npts)
            nxy,nopt=d.shape
            del d
            t0=fns[0].split('.')[-2]
            t1=fns[-1].split('.')[-2]
            print t0,t1
            dds=np.memmap(casename+'.%s.%s.%s.data'%(var,t0,t1),dtype='>f4',shape=(len(fns),n,nopt),mode='write')
            for i in range(len(fns)):
                d=np.fromfile(fns[i],'>f4').reshape(n,-1)
                f.writelines(fns[i])
                dds[i,...]=d
                del d
            f.close()
            del dds,

max_stars_repo_name: jinbow/Octupus

max_issues_repo_name: jinbow/Octupus

max_forks_repo_name: jinbow/Octupus


text:
"""
UAV Toolbox
"""

# Dependencies
import arcpy
import os
import cv2
import numpy as np
import scipy as sp
import sklearn

class Toolbox(object):
    def __init__(self):
        """
	Define the toolbox (the name of the toolbox is the name of the .pyt file).
	"""
        self.label = "Toolbox"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = [Tool, RGB2HSV]

class Tool(object):
    def __init__(self):
        """
	Define the tool (tool name is the name of the class).
	"""
        self.label = "Tool"
        self.description = ""
        self.canRunInBackground = False

    def getParameterInfo(self):
        """
	Define parameter definitions
	"""
        params = None
        return params

    def isLicensed(self):
        """
	Set whether tool is licensed to execute.
	"""
        return True

    def updateParameters(self, parameters):
        """
	Modify the values and properties of parameters before internal
        validation is performed.  This method is called whenever a parameter
        has been changed.
	"""
        return

    def updateMessages(self, parameters):
        """
	Modify the messages created by internal validation for each tool
        parameter.  This method is called after internal validation.
	"""
        return

    def execute(self, parameters, messages):
        """
	The source code of the tool.
	"""
        return

class RGB2HSV(object):
    def __init__(self):
        """
	Define the tool (tool name is the name of the class).
	"""
        self.label = "RGB to HSV"
        self.description = ""
        self.canRunInBackground = False

    def getParameterInfo(self):
        """
	Define parameter definitions
	"""

        # Input
        input_raster = arcpy.Parameter(
            displayName="Input RGB Raster",
            name="in_raster",
            datatype="DERasterDataset",
            parameterType="Required",
            direction="Input")

        # Output
        output_raster = arcpy.Parameter(
            displayName="Output HSV Raster",
            name="out_raster",
            datatype="DERasterDataset",
            parameterType="Required",
            direction="Output")

        output_raster.parameterDependencies = [input_raster.name]
        output_raster.schema.clone = True

        params = [input_raster, output_raster]
        return params

    def isLicensed(self):
        """
	Set whether tool is licensed to execute.
	"""
        return True

    def updateParameters(self, parameters):
        """
	Modify the values and properties of parameters before internal
        validation is performed.  This method is called whenever a parameter
        has been changed.
	"""
        return

    def updateMessages(self, parameters):
        """
	Modify the messages created by internal validation for each tool
        parameter.  This method is called after internal validation.
	"""
        return

    def execute(self, parameters, messages):
        """
	The source code of the tool.
	"""
        input_raster = parameters[0].valueAsText
        output_raster = parameters[1].valueAsText
        arcpy.AddMessage("Input Raster: %s" % input_raster)
        arcpy.AddMessage("Output Raster: %s" % output_raster)

        RGB_array = arcpy.RasterToNumPyArray(input_raster)
        R_array = RGB_array[0,:,:]
        G_array = RGB_array[1,:,:]
        B_array = RGB_array[2,:,:]
        RGB_array = np.dstack((R_array, G_array, B_array))
        HSV_array = cv2.cvtColor(RGB_array, cv2.COLOR_RGB2HSV)
        arcpy.AddMessage("Generated HSV")

        Ymin = float(str(arcpy.GetRasterProperties_management(input_raster, "BOTTOM")))
        Xmin = float(str(arcpy.GetRasterProperties_management(input_raster, "LEFT")))
        Xmax = float(str(arcpy.GetRasterProperties_management(input_raster, "RIGHT")))
        Ymax = float(str(arcpy.GetRasterProperties_management(input_raster, "TOP")))
        Xcell = float(str(arcpy.GetRasterProperties_management(input_raster, "CELLSIZEX")))
        Ycell = float(str(arcpy.GetRasterProperties_management(input_raster, "CELLSIZEY")))
        xymin = arcpy.Point(Xmin, Ymin)
        HSV_raster = arcpy.NumPyArrayToRaster(HSV_array, xymin, Xcell, Ycell)
        #HSV_raster.save(output_raster)
        #arcpy.CopyRaster_management(HSV_raster, output_raster)
        return

ext: pyt

max_stars_repo_name: trevstanhope/arc-cropland

max_issues_repo_name: trevstanhope/arc-cropland

max_forks_repo_name: trevstanhope/arc-cropland

reason: import numpy,import scipy


text:
import pytest
import numpy as np
import lenstronomy.LightModel.Profiles.ellipsoid as torus
from lenstronomy.LightModel.Profiles.ellipsoid import Ellipsoid

class TestTorus(object):
    """
    class to test Shapelets
    """
    def setup(self):
        pass

    def test_function(self):
        """

        :return:
        """
        output = torus.function(x=1, y=1, amp=1., sigma=2, center_x=0, center_y=0)
        assert output == 0.079577471545947673

class TestEllipsoid(object):

    def setup(self):
        pass

    def test_function(self):
        """

        :return:
        """
        ellipsoid = Ellipsoid()
        output = ellipsoid.function(x=1, y=1, amp=1., radius=1, e1=0, e2=0, center_x=0, center_y=0)
        assert output == 0
        output = ellipsoid.function(x=0.99, y=0, amp=1., radius=1, e1=0, e2=0, center_x=0, center_y=0)
        assert output == 1./np.pi

if __name__ == '__main__':
    pytest.main()


text:
"""
Swendson & Wang Replica Monte Carlo for general ising model.

Author: Colin Clement
Date: 2016-6-11

This implementation follows: arXiv:cond-mat/0407273v1
The only modification is a general implementation for any
graph with two-spin interactions.

"""

import numpy as np
import scipy.sparse as sps
from itertools import islice
from collections import defaultdict
from functools import wraps
from subprocess import check_output
from math import exp
from scipy.spatial.distance import squareform
try:
    import cPickle as pickle
except ImportError:
    import pickle
try:
    from numba import jit
    hasnumba = True
except ImportError as perror:
    hasnumba = False

class GeneralSpinModel(object):
    """Swendsen and Wang replica Monte Carlo for ising model with arbitrary
    couplings and interaction graphs"""
    def __init__(self, J, templist = [1.], observers = None, sparse = True, **kwargs):
        """
            Monte carlo simulation for generalized ising couplings.
            Parameters:
                J: coupling array of size (N_spins, N_spins), ideally sparse
                templist: list of temperatures to simulate, should be ordered.
                observers : list of observer functions
            kwargs:
                seed : Seed for initializing random number generator
                therm : Initially thermalize for integer number of steps provided
            Kwargs:
        """
        if len(J.shape)==1:
            self.J = squareform(J)
        elif (J-J.T).sum() > 0:
            raise(RuntimeError, "Couplings must be symmetric")
        else:
            self.J = J
        if sparse:
            self.J = sps.csr_matrix(self.J)
        self.N = J.shape[0]
        self.templist = templist
        self.Ntemps = len(templist)
        self.time = 0 #increments each MC step
        self.seed = kwargs.get('seed', 92089)
        self.therm = kwargs.get('therm', 0)
        self.observers = observers or []

        # Data structures of J matrix that are more efficient
        self.bonds = defaultdict(set)
        self.neighborsets = [neighbors(i, self.J) for i in range(self.N)]
        self.neighborlist = []
        self.neighborbonds = []
        self.bondmap = {}
        for i, neighs in enumerate(self.neighborsets):
            self.neighborlist += [np.array(list(neighs))]
            self.neighborbonds += [self.J[i, list(neighs)].toarray().flatten()]
            for n in neighs:
                if i < n: self.bonds[i].update(set([(i,n)]))
                else: self.bonds[i].update(set([(n,i)]))
                self.bondmap[(i,n)] = self.J[i,n]

        self.cumSpinPair = np.zeros((self.N, self.N))
        self.pairtime = 0.
        self._seeded_rand() #creates self._my_rng random number generator
        self.spinarray = 2*(self._my_rng.rand(len(self.templist), self.N)>0.5) - 1

        if self.therm: self.thermalize(self.therm)

    def __repr__(self):
        formatstring = "<{cls}(Ntemps={Ntemps}, N={N}, time = {time}, seed = {seed})>"
        return formatstring.format(cls=self.__class__.__name__, **self.__dict__)

    def __getstate__(self):
        """
        for pickling, observers are not pickleable
        """
        gooddict = {key:value for key,value in self.__dict__.iteritems() if not
                    key=='observers'}
        gooddict['observers'] = []
        return gooddict

    def __iter__(self):
        """ Make an instance of this class an iterator. See the next method."""
        return self

    def _seeded_rand(self):
        """Set up random number generator"""
        self._my_rng = np.random.RandomState()
        self._my_rng.seed(self.seed)

    def _increment_time(self):
        """Count time steps (total MC sweeps)"""
        self.time += 1
        self.notify_observers() #always notify, observer decides if time is right

    def thermalize(self, relax_time):
        """
        performs relax_time (int) monte carlo steps
        """
        for state in islice(self,relax_time):
            pass

    def energy(self, s, J = None):
        """Return the energy of the current state"""
        J = J or self.J
        return -s.dot(J.dot(s))/2.

    def _oneSweep(self, s, T):
        """
        One Monte Carlo step for each spin in s
        input:
            s : spin array of length self.N
            T : temperature float
        """
        randtest = self._my_rng.rand(self.N)
        for i, test in enumerate(randtest):
            dE = 2*s[i]*(self.J[i].dot(s))
            if dE < 0 or test <= exp(-dE/T): s[i] *= -1

    def clusterStep(self, s0, T0, s1, T1):
        """
        Find clusters of constant overlap between
        spin configurations s0 and s1, perform one sweep
        of replica cluster Monte Carlo.
        input:
            s0: spin configuration 0
            T0: Temperature of spin configuration 0
            s1: spin configuration 1
            T1: Temperature of spin configuration 1
        """
        tau = s0*s1
        cls = getClusters(tau, self.neighborsets)
        k_ab = (1./T0 - 1./T1)*self.getClusterCouplings(cls, s0)
        eta = np.ones(len(cls))
        randtest = self._my_rng.rand(len(cls))
        for i, test in enumerate(randtest):
            dE = 2*eta[i]*(k_ab[i].dot(eta))
            if dE < 0 or test <= exp(-dE):
                eta[i] *= -1
                cl = cls[i]
                s0[cl], s1[cl] = -s0[cl], -s1[cl]

    def getClusterCouplings(self, clusters, s0):
        """
        Get the effective cluster couplings as prescribed
        by Swendsen and Wang.
        input:
            clusters: list of lists of cluster indicies
            s0: reference spin configuration
        returns:
            sparse matrix of effective couplings
        """
        bdrys = [boundaries(cls, self.bonds) for cls in clusters]
        bondmap = self.bondmap
        rows, cols, data = [], [], []
        for i, bdry_i in enumerate(bdrys[:-1]):
            for jj, bdry_j in enumerate(bdrys[i+1:]):
                shared = bdry_i.intersection(bdry_j)
                if shared:
                    coupling = sum([s0[bd[0]]*s0[bd[1]]*bondmap[bd] for bd in shared])
                    rows += [i, jj+i+1]
                    cols += [jj+i+1, i] #symmetric matrix
                    data += [coupling, coupling]
        N_c = len(clusters)
        return sps.csr_matrix((data, (rows, cols)), shape=(N_c, N_c))

    def next(self):
        """ For the iterator, do one monte carlo step"""
        for nt, T in enumerate(self.templist):
            for i, rand in enumerate(self._my_rng.rand(self.N)):
                flipFast(i, self.spinarray[nt], T, rand,
                         self.neighborlist[i], self.neighborbonds[i])
        for nt in range(self.Ntemps-1):
            self.clusterStep(self.spinarray[nt], self.templist[nt],
                             self.spinarray[nt+1], self.templist[nt+1])
        self._increment_time()
        return self

    def notify_observers(self):
        """ Send state to observers """
        for obs in self.observers:
            obs.send(self)

#-----------------------------
# Utilities
#-----------------------------

def neighbors(ind, adj):
    """
    Given adjacency matrix, return the neighbors
    of node 'ind'
    """
    if type(adj) is sps.csr.csr_matrix:
        return set(adj[ind].indices)
    elif type(adj) is np.ndarray:
        conn = adj[ind]
        return set(np.arange(len(conn))[adj[ind]>0])

def boundaries(cls, bonds):
    return reduce(lambda s, t: s.symmetric_difference(t),
                  [bonds[i] for i in cls])

def getClusters(taus, neighborsets):
    """
    input:
        taus : values of overlap between two spin configurations.
               Array of N spin-values.
        neighborsets: list ordered by spin index, whose
        elements are the spins neighboring the spin at that index.
    output:
        list of cluster indices split according to the adjacency
        matrix to find contiguous regions of 'up' or 'down' spins
    """
    result = []
    spins = set(range(len(taus)))
    while spins:
        seed = spins.pop()
        newclust = {seed}
        queue = [seed]
        while queue:
            nxt = queue.pop(0)
            neighs = neighborsets[nxt].difference(newclust)
            for neigh in list(neighs):
                if not taus[neigh] == taus[seed]: neighs.remove(neigh)
            newclust.update(neighs)
            queue.extend(neighs)
            spins.difference_update(neighs)
        result.append(list(newclust))
    return result

@jit("void(int64, int64[:], float64, float64, int32[:], float64[:])",
     nopython=True)
def flipFast(index, s, T, rand, neighs, bonds):
    """
    Much more efficient version of _oneSweep MC
    input:
        index: spin to attempt flipping
        s: spin configuration
        T: temperature
        rand: a uniform random float
        neighs: list of neighboring sites to index
        bonds: list of bonds connections index to its neighbors
    """
    dE = 0.
    for nn, jj in zip(neighs, bonds):
        dE += s[nn]*jj
    dE *= 2*s[index]
    if dE < 0 or rand <= exp(-dE/T): s[index] *= -1

def coroutine(func):
    """ Decorator that initializes coroutines (return yield)"""
    @wraps(func)
    def warmed_up(*args, **kwargs):
        cr = func(*args, **kwargs)
        cr.send(None)
        return cr
    return warmed_up

def fileName(model, directory = 'time_series/', protocol = '.pkl',  appendnote = ''):
    """Creates a unique file name for a given model final string first 10
    characters from the current git commit hash"""
    hashstr = str(check_output("git log -n 1 | grep commit | sed s/commit\ //", shell=True)[:10])
    strformat = 'EA-seed_{seed}-size_{size}-T_{T}-' + appendnote + hashstr
    return directory + strformat.format(**model.__dict__) + protocol

def pickleLoader(pklfile):
    try:
        while True:
            yield pickle.load(pklfile)
    except EOFError:
        pass

def timeSeriesLoader(filename):
    series = []
    with open(filename, 'rb') as infile:
        for state in pickleLoader(infile):
            series += [state]
    return np.array(series)

def packSpins(spins):
    return np.packbits((spins+1)/2)

def unpackSpins(ints, N):
    return (2*np.unpackbits(ints).astype('int')-1)[:N]

#----------------------------
# Observers
#----------------------------

@coroutine
def saveTimeSeries(period = 1, directory = 'time_series/', protocol = '.npy', **kwargs):
    """Observer which saves spin state of model to a file. specify directory to save in,
    protocol: 'pkl' or 'h5 for pickling or hdf5"""
    #Make file upon recieving model with fileName()
    model = yield
    if protocol == '.pkl':
        name = fileName(model, directory, protocol, **kwargs)
        if 'filename' in kwargs:
            name = kwargs['filename']
        while True:
            model = yield
            if model.time % period == 0:
                with open(name,'ab') as outfile:
                    pickle.dump(model.spins.copy(), outfile,-1)
    elif protocol == '.npy':
        with open(name, 'wb') as outfile:
            while True:
                model = yield
                if model.time % period == 0:
                    np.save(outfile, packSpins(model.spins))
                    outfile.flush()

@coroutine
def onlinePairCorr(period=100, **kwargs):
    """Calculate the pair correlation function using the same online
    Welford-type algorithm as for the FIM. Only calculates cross moments.
    PairCorr is calculated upon pickling by dividing by self.pairtime and T**2.
    """
    model = yield
    model.pairtime = 0
    N = model.N
    model.meanSpinPair = np.zeros(N)
    model.cumSpinPair = np.zeros((N, N))
    while True:
        model = yield
        if model.time % period==0:
            model.pairtime += 1
            meanSpin_old = model.meanSpinPair.copy().astype('float')
            spins = model.spins.copy().astype('float')
            model.meanSpinPair += (spins - meanSpin_old)/(model.pairtime+0.0)
            model.cumSpinPair += np.outer(spins-model.meanSpinPair,
                                          spins-meanSpin_old)

#----------------------------
#   SuperCluster Calculations
#----------------------------

def simplifyNestedClustering(supercluster, cluster):
    """
    Unpacks a clustering of clusters. supercluster is a list of indices from
    cluster, cluster is a list of original indicies. The output is the indices
    of cluster organizes according to supercluster
    """
    flattened = []
    for scl in supercluster:
        flat_super = []
        for cl in scl:
            flat_super += cluster[cl]
        flattened += [flat_super]
    return flattened

def mapCij(clust, CIJ):
    """
    Takes a ClusterModel correlation function and, using the cluster assignment
    inherited from an REMC model, maps the correlation function (CIJ) onto
    the dimensions of the original system correlations
    """
    N = np.sum(map(len, clust))
    mapped_CIJ = np.zeros((N,N))
    for cl1, row in zip(clust, CIJ):
        for cl2, corr in zip(clust, row):
            for i_x in cl1:
                for i_y in cl2:
                    mapped_CIJ[i_x,i_y] = corr
    return mapped_CIJ

max_stars_repo_name: bcdaniels/coniii

max_issues_repo_name: bcdaniels/coniii

max_forks_repo_name: bcdaniels/coniii

reason: import numpy,import scipy,from scipy,from numba


text:
import numpy as np
import warnings
from sage import utils

class Imputer:
    '''Imputer base class.'''
    def __init__(self, model):
        self.model = utils.model_conversion(model)

    def __call__(self, x, S):
        raise NotImplementedError

class DefaultImputer(Imputer):
    '''Replace features with default values.'''
    def __init__(self, model, values):
        super().__init__(model)
        if values.ndim == 1:
            values = values[np.newaxis]
        elif values[0] != 1:
            raise ValueError('values shape must be (dim,) or (1, dim)')
        self.values = values
        self.values_repeat = values
        self.num_groups = values.shape[1]

    def __call__(self, x, S):
        # Prepare x.
        if len(x) != len(self.values_repeat):
            self.values_repeat = self.values.repeat(len(x), 0)

        # Replace specified indices.
        x_ = x.copy()
        x_[~S] = self.values_repeat[~S]

        # Make predictions.
        return self.model(x_)

class MarginalImputer(Imputer):
    '''Marginalizing out removed features with their marginal distribution.'''
    def __init__(self, model, data):
        super().__init__(model)
        self.data = data
        self.data_repeat = data
        self.samples = len(data)
        self.num_groups = data.shape[1]

        if len(data) > 1024:
            warnings.warn('using {} background samples may lead to slow '
                          'runtime, consider using <= 1024'.format(
                            len(data)), RuntimeWarning)

    def __call__(self, x, S):
        # Prepare x and S.
        n = len(x)
        x = x.repeat(self.samples, 0)
        S = S.repeat(self.samples, 0)

        # Prepare samples.
        if len(self.data_repeat) != self.samples * n:
            self.data_repeat = np.tile(self.data, (n, 1))

        # Replace specified indices.
        x_ = x.copy()
        x_[~S] = self.data_repeat[~S]

        # Make predictions.
        pred = self.model(x_)
        pred = pred.reshape(-1, self.samples, *pred.shape[1:])
        return np.mean(pred, axis=1)

max_stars_repo_name: shubhaguha/sage

reason: import numpy,from sage


text:
import torch
import json
import nltk
import pickle
import numpy as np
from tqdm import tqdm
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
import nltk
nltk.download("punkt")
from nltk.tokenize import word_tokenize
import gensim
from sklearn.svm import SVC
from sklearn.metrics import f1_score
from sys import exit

def reader(file):
    data = []
    with open(f'../data/{file}', 'r') as file:
        for line in file.readlines()[:]:
            data.append(json.loads(line))
    return data

dev, test, train = [reader(file) for file in ['music_reviews_dev.json',
                                              'music_reviews_test_masked.json',
                                              'music_reviews_train.json']]

def tokenize(file):
        X = []
        y = []
        for review in file:
            if 'reviewText' in review.keys():
                X.append(word_tokenize(review['reviewText']))
                y.append(1)
            else:
                y.append(-1)
                X.append("_")

        return X, y

test_X, test_y = tokenize(test)

#print("LOADING EMBEDDINGS")
model = gensim.models.KeyedVectors.load_word2vec_format('../data/twitter.bin', binary=True)
vocab = set(model.vocab.keys())
#print("EMBEDDINGS LOADED")

def embed(data):
    out = []
    for review in tqdm(data):
        tmp = []
        for word in review:
            tmp.append(model[word]) if word in vocab else tmp.append(model['<U>'])
        t = np.mean(np.array(tmp),axis=0)
        out.append(t)
    return out

#print("STARTED HERE")

test_X = embed(test_X)

pred = []
counter = 0
#print("PREDICTING")
clf = pickle.load(open("SVM","rb"))
for x_inst,y_inst in zip(test_X,test_y):
    if y_inst == -1:
        counter += 1
        pred.append(1)
    else:
        pred.append(clf.predict(x_inst.reshape(1,-1))[0])

#print("NO TEXT DATA COUNTER:",counter)
print("id,prediction")
for i in range(len(pred)):
    print(f"{test[i]['id']},{pred[i]}")

max_stars_repo_name: syrkis/neurolingo

max_issues_repo_name: syrkis/neurolingo

max_forks_repo_name: syrkis/neurolingo


text:
#!/bin/python

import numpy as np
from astropy.io import fits
import rca
from joblib import Parallel, delayed, parallel_backend, cpu_count

import click
@click.command()

@click.option(
    "--n_comp",
    default=4,
    type=int,
    help="RCA number of eigenPSFs.")
@click.option(
    "--upfact",
    default=1,
    type=int,
    help="Upsampling factor.")
@click.option(
    "--ksig",
    default=3.,
    type=float,
    help="Denoising parameter K.")
@click.option(
    "--run_id",
    default="rca",
    type=str,
    help="Id used for the saved models and validations.")
@click.option(
    "--psf_size",
    default=1.25,
    type=float,
    help="PSF size in the type of the psf_size_type parameter.")
@click.option(
    "--psf_size_type",
    default="sigma",
    type=str,
    help="PSF size type.")
@click.option(
    "--saving_dir",
    default="/n05data/tliaudat/wf_exps/outputs/rca/",
    type=str,
    help="Path to the saving directory. Should include the directories /models, /validation and /metrics.")
@click.option(
    "--input_data_dir",
    default="/n05data/tliaudat/wf_exps/datasets/rca_shifts/",
    type=str,
    help="Input dataset directory. Should have /train and /test directories.")

def main(**args):
    print(args)
    with parallel_backend("loky", inner_max_num_threads=1):
        results = Parallel(n_jobs=1)(
            delayed(rca_procedure)(**args) for i in range(1)
        )

def rca_procedure(**args):
    # Model parameters
    n_comp = args['n_comp']
    upfact = args['upfact']
    ksig = args['ksig']

    run_id = args['run_id']

    # Load data
    saving_base_path = args['saving_dir']
    input_base_path = args['input_data_dir']
    model_save_dir_path = saving_base_path + 'models/'
    val_save_dir_path = saving_base_path + 'validation/'
    metrics_save_dir_path = saving_base_path + 'metrics/'

    catalog_ids = [200, 500, 1000, 2000]  # [200]
    ccd_tot = 36  # 3

    rmse_list = []
    rel_rmse_list = []

    # Iterate over the catalogs
    for catalog_id in catalog_ids:
        # Create lists for exposure metrics
        psf_interp_list = []
        star_list = []

        # Iterate over the ccds
        for ccd_it in range(ccd_tot):

            ## Training
            # Train data
            train_file_name = 'train/train_stars-%07d-%02d.fits'%(catalog_id, ccd_it)
            train_catalog = fits.open(input_base_path + train_file_name)

            # Prepare input data
            # In RCA format (batch dim at the end)
            obs_stars = rca.utils.rca_format(
                train_catalog[1].data['VIGNET']
            )
            obs_pos = np.array([
                train_catalog[1].data['XWIN_IMAGE'],
                train_catalog[1].data['YWIN_IMAGE']
            ]).T

            # Create RCA instance
            rca_inst = rca.RCA(
                n_comp=n_comp,
                upfact=upfact,
                ksig=ksig,
                verbose=1
            )

            # If we have 2 stars or less we skip the ccd
            if obs_pos.shape[0] <= 2:
                continue

            # fit it to stars
            _, _ = rca_inst.fit(
                obs_stars,
                obs_pos,
                psf_size=args['psf_size'],
                psf_size_type=args['psf_size_type']
            )
            # Save model
            save_path = model_save_dir_path + run_id + '_fitted_model-%07d-%02d.npy'%(catalog_id, ccd_it)
            rca_inst.quicksave(save_path)

            ## Validation
            # Test data
            test_file_name = 'test/test_stars-%07d-%02d.fits'%(catalog_id, ccd_it)
            test_catalog = fits.open(input_base_path + test_file_name)

            test_stars = rca.utils.rca_format(
                test_catalog[1].data['VIGNET']
            )

            # Test positions
            test_pos = np.array([
                test_catalog[1].data['XWIN_IMAGE'],
                test_catalog[1].data['YWIN_IMAGE']
            ]).T

            # Recover test PSFs
            interp_psfs_upfact1 = rca_inst.validation_stars(test_stars, test_pos, upfact=1)
            interp_psfs = rca_inst.validation_stars(test_stars, test_pos, upfact=upfact)

            # Save validation PSFs and Ground truth stars
            val_dict = {
                'PSF_VIGNETS': interp_psfs,
                'PSF_UPFACT1': interp_psfs_upfact1,
                'GT_VIGNET': test_catalog[1].data['VIGNET'],
                'POS': test_pos,
            }
            val_save_name = run_id + '_validation-%07d-%02d.npy'%(catalog_id, ccd_it)
            np.save(val_save_dir_path + val_save_name, val_dict, allow_pickle=True)

            # Add images to lists
            psf_interp_list.append(np.copy(interp_psfs))
            star_list.append(np.copy(test_catalog[1].data['VIGNET']))

        # Calculate RMSE metric on all the CCDs
        exp_psfs = np.concatenate(psf_interp_list)
        exp_stars = np.concatenate(star_list)

        # Calculate residuals
        residuals = np.sqrt(np.mean((exp_psfs - exp_stars)**2, axis=(1,2)))
        GT_star_mean = np.sqrt(np.mean((exp_stars)**2, axis=(1,2)))

        # RMSE calculations
        rmse = np.mean(residuals)
        rel_rmse = 100. * np.mean(residuals/GT_star_mean)

        rmse_list.append(rmse)
        rel_rmse_list.append(rel_rmse)

    rmse_list = np.array(rmse_list)
    rel_rmse_list = np.array(rel_rmse_list)

    # Save the metrics
    metrics_dict = {
        'rmse': rmse_list,
        'rel_rmse': rel_rmse_list
    }
    metrics_save_name = run_id + '_metrics.npy'
    np.save(metrics_save_dir_path + metrics_save_name, metrics_dict, allow_pickle=True)

if __name__ == "__main__":
  main()

max_stars_repo_name: tobias-liaudat/wf-psf

max_issues_repo_name: tobias-liaudat/wf-psf

max_forks_repo_name: tobias-liaudat/wf-psf

reason: import numpy,from astropy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Written by Tong He and CBIG under MIT license:
https://github.com/ThomasYeoLab/CBIG/blob/master/LICENSE.md
"""

import os
import argparse
import numpy as np
import scipy.io as sio
from shutil import copyfile

def get_phes(folder, dataset):
    '''get the phenotypes of test meta-set

    Args:
        folder (str): folder path that stores data input files
        dataset (str): state the test meta-set dataset

    Returns:
        List: list of test meta-set phenotypes
    '''

    if dataset == 'ukbb':
        txt = os.path.join(folder, 'ukbb_test_final_phe_list.txt')
    elif dataset == 'HCP':
        txt = os.path.join(folder, 'HCP_diff_roi_final_phe_list.txt')
    elif dataset == 'exp':
        txt = os.path.join(folder, 'exp_test_final_phe_list.txt')
    else:
        raise NameError('wrong dataset name')
    phes = np.genfromtxt(txt, dtype=str)
    return phes

def get_split(dataset='ukbb'):
    '''save the split files from KRR classical and meta-matching

    Args:
        dataset (str): state the test meta-set dataset

    Returns:
        None
    '''

    if dataset == 'exp':
        rngs = 3
        n_subjects = 1000
        base_dir = os.path.join(
            os.environ['CBIG_CODE_DIR'],
            'stable_projects/predict_phenotypes/He2022_MM/examples/exp_output')
        input_dir = os.path.join(
            os.environ['CBIG_CODE_DIR'],
            'stable_projects/predict_phenotypes/He2022_MM/examples/exp_input')
        classical_dir = os.path.join(base_dir, 'output_KRR_classical_exp')
        mm_dir = os.path.join(base_dir, 'output_KRR_mm')
    elif dataset == 'unit_tests':
        rngs = 2
        n_subjects = 400
        base_dir = os.path.join(
            os.environ['CBIG_CODE_DIR'],
            'stable_projects/predict_phenotypes/He2022_MM/unit_tests/output')
        input_dir = os.path.join(os.environ['CBIG_CODE_DIR'],
                                 'stable_projects', 'predict_phenotypes',
                                 'He2022_MM', 'examples', 'unit_tests_input')
        classical_dir = os.path.join(base_dir, 'output_KRR_classical_exp')
        mm_dir = os.path.join(base_dir, 'output_KRR_mm')
        dataset = 'exp'
    else:
        rngs = 100
        base_dir = os.path.join(
            os.environ['CBIG_CODE_DIR'],
            'stable_projects/predict_phenotypes/He2022_MM/replication')
        input_dir = os.path.join(
            os.environ['CBIG_REPDATA_DIR'],
            'stable_projects/predict_phenotypes/He2022_MM')
        if dataset == 'ukbb':
            n_subjects = 10000
            classical_dir = os.path.join(base_dir, 'output_KRR_classical_ukbb')
            mm_dir = os.path.join(base_dir, 'output_KRR_mm')
        elif dataset == 'HCP':
            n_subjects = 1019
            classical_dir = os.path.join(base_dir, 'output_KRR_classical_HCP')
        else:
            raise NameError('wrong dataset name ' + dataset)

    output_dir = os.path.join(base_dir, 'output_intermediate')
    os.makedirs(output_dir, exist_ok=True)

    phes_tes = get_phes(input_dir, dataset)
    ks = [10, 20, 50, 100, 200]
    split_ind_dict = {}

    for _, phe in enumerate(phes_tes):
        split_ind_dict[phe] = np.zeros((rngs, len(ks), n_subjects))

    # get split file for k shot
    for _, phe in enumerate(phes_tes):
        for ik, k in enumerate(ks):
            for rng in range(1, rngs + 1):
                folder = os.path.join(
                    classical_dir, 'output_phe_' + phe, dataset + '_' + phe +
                    '_k_' + str(k) + '_rng_num_' + str(rng))
                split_mat = os.path.join(folder,
                                         dataset + '_subject_split.mat')
                mat = sio.loadmat(split_mat)
                split_ind_dict[phe][rng - 1, ik, :] = np.squeeze(
                    mat['sub_fold']['fold_index'][0][0])
        print(phe)

    for _, phe in enumerate(phes_tes):
        split_ind_dict[phe] = split_ind_dict[phe].astype(bool)

    npy = os.path.join(output_dir,
                       dataset + '_split_ind_' + str(rngs) + '.npz')
    np.savez(npy, split_ind_dict=split_ind_dict)

    # get split file from KRR base
    if dataset != 'HCP':
        rngs = 1
        for rng in range(1, rngs + 1):
            folder = os.path.join(mm_dir, dataset + '_rng_num_' + str(rng))
            split_mat = os.path.join(folder, dataset + '_subject_split.mat')
            target_mat = os.path.join(output_dir,
                                      'split_rng' + str(rng) + '.mat')
            copyfile(split_mat, target_mat)

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument("-d", "--dataset", type=str, default="ukbb")
    args = parser.parse_args()
    get_split(args.dataset)


text:
import numpy as np
import tensorflow as tf
import random

from rl_dnc.model.dnc_wrapper import DNCWrapper

sequence_length = 3
# used so that only the last value in the target is of importance
sequence_mask = np.array([[0] for i in range(sequence_length-1)] + [[1]])

# simple count(1s)%2 == ?
def generate_data(count):
    inputs = np.zeros((count, sequence_length, 1, 1))
    outputs = np.zeros_like(inputs)
    for i in range(count):
        sequence = [random.randint(0, 1) for i in range(sequence_length)]
        observation = np.zeros((sequence_length,1,1))
        observation[:, 0, 0] = sequence
        expected = np.zeros_like(observation)
        expected[-1, 0, 0] = sum(sequence) % 2
        inputs[i] = observation
        outputs[i] = expected
    return inputs, outputs

def cost_fun(prediction, target, mask):
    diff = tf.squared_difference(prediction, target)
    loss_time_batch = tf.reduce_sum(diff, axis=2)
    loss_batch = tf.reduce_sum(loss_time_batch * mask, axis=0)
    batch_size = tf.cast(tf.shape(prediction)[1], dtype=loss_time_batch.dtype)
    loss = tf.reduce_sum(loss_batch) / batch_size
    return loss

def xent_cost_fun(prediction, target, mask):
    xent = tf.nn.sigmoid_cross_entropy_with_logits(labels=target, logits=prediction)
    loss_time_batch = tf.reduce_sum(xent, axis=2)
    loss_batch = tf.reduce_sum(loss_time_batch * mask, axis=0)
    batch_size = tf.cast(tf.shape(prediction)[1], dtype=loss_time_batch.dtype)
    loss = tf.reduce_sum(loss_batch) / batch_size
    return loss

def main(unused_argv):
    input_size = 1
    output_size = 1
    batch_size = 1
    iterations = 5
    report_interval = 100

    inputs, outputs = generate_data(iterations)
    costfun_prediction = tf.placeholder(shape=inputs[0].shape, dtype=tf.float32)
    costfun_target = tf.placeholder(shape=outputs[0].shape, dtype=tf.float32)
    costfun_mask = tf.placeholder(shape=sequence_mask.shape, dtype=tf.float32)

    costfun_node = cost_fun(costfun_prediction, costfun_target, costfun_mask)
    xent_costfun_node = xent_cost_fun(costfun_prediction, costfun_target, costfun_mask)

    with tf.Session() as sess:

        for i in range(iterations):
            obs = inputs[i]
            target = outputs[i]
            out, xent_out = sess.run([costfun_node, xent_costfun_node], feed_dict={
                costfun_prediction: obs,
                costfun_target: target,
                costfun_mask: sequence_mask,
            })
            print('observation: ', obs[:,0,0])
            print('target: ', target[:,0,0])
            print('mask: ', sequence_mask[:,0])
            print('costfun: ', out)
            print('xentcostfun: ', xent_out)
            print('\n')

if __name__ == "__main__":
    tf.app.run()


text:
#!/usr/bin/env python
# -*- coding:utf-8 -*-
"""
@author: peter.s
@project: EmptyModel
@time: 2019/11/5 21:01
@desc:
"""
from abc import abstractmethod
import numpy as np

from aehn.lib import yield2batch_data, get_metrics_callback_from_names
from aehn.lib import DictScaler, VoidScaler, ZeroMaxScaler, SingletonStandScaler, SingletonZeroMaxScaler

class AbstractDataProvider(object):

    @abstractmethod
    def iterate_batch_data(self):
        """ Get batch model input of one epoch.
        Remark: batch -> partition -> epoch
        :return: yield a list containing batch inputs until the end of the epoch.
        """
        pass

    @abstractmethod
    def get_metrics(self, preds, labels):
        """ Calculate the metrics of preds and labels.
        :param preds:
        :param labels:
        :return: a dictionary of the metrics.
        """
        pass

    @abstractmethod
    def epoch_inverse_scaling(self, scaled_records):
        """ Inverse the scaled_records to real scale.
        :param scaled_records:
        :return: real scale records.
        """
        pass

class DataProvider(AbstractDataProvider):
    """
    Data provider for processing model inputs.
    """
    def __init__(self, data_source, data_config, scaler=None):
        self._data_source = data_source
        self._batch_size = data_config['batch_size']
        self._metrics_function = get_metrics_callback_from_names(data_config['metrics'])
        self._scaler = scaler if scaler else DictScaler(dtimes=VoidScaler, marks=SingletonStandScaler)
        self._is_first_iterate = True

        self._type_padding = data_config['process_dim']

    def get_scaler(self):
        return self._scaler

    def epoch_inverse_scaling(self, scaled_records):
        return self._scaler.inverse_scaling(scaled_records)

    def get_metrics(self, preds, labels):
        seq_mask = labels['types'] < self._type_padding
        return self._metrics_function(preds, labels, seq_mask=seq_mask)

    def iterate_batch_data(self):
        # record_data of a partition whose shape is [n_records, ...]
        for data in self._data_source.load_partition_data():
            if self._is_first_iterate:
                data_stats = self._dataset_statistics(data)
                print(f'Load dataset {self._data_source.data_name}: {data_stats}')

            inputs = self._process_model_input(data)
            if self._scaler.is_fit():
                scaled_inputs = self._scaler.scaling(inputs)
            else:
                scaled_inputs = self._scaler.fit_scaling(inputs)

            # yield records to batch data separately
            for batch_data in yield2batch_data(scaled_inputs, self._batch_size, keep_remainder=True):
                yield batch_data

        if self._is_first_iterate:
            self._is_first_iterate = False

    def _process_model_input(self, records):
        """ Process each item as model input.
        :param records:
        :return:
        """
        type_padding = self._type_padding
        dt_padding = 0.0

        type_seqs = records['types']
        time_seqs = records['timestamps']
        # dt_i = t_i - t_{i-1}
        dt_seqs = [[t_seq[i] - t_seq[max(i - 1, 0)] for i in range(len(t_seq))]
                   for t_seq in time_seqs]

        n_records = len(type_seqs)
        max_len = max([len(seq) for seq in type_seqs])

        # padding
        type_seqs_padded = np.ones([n_records, max_len]) * type_padding
        dt_seqs_padded = np.ones([n_records, max_len]) * dt_padding

        for i in range(n_records):
            len_seq = len(type_seqs[i])
            type_seqs_padded[i, :len_seq] = type_seqs[i]
            dt_seqs_padded[i, : len_seq] = dt_seqs[i]

        ret = dict()
        ret['types'] = type_seqs_padded
        ret['dtimes'] = dt_seqs_padded

        # add to model input if original input has marks
        if 'marks' in records.keys():
            marks_seqs = records['marks']
            marks_padding = 0.0
            marks_padded = np.ones([n_records, max_len, 1]) * marks_padding
            for i in range(n_records):
                len_seq = len(marks_seqs[i])
                marks_padded[i, : len_seq, 0] = marks_seqs[i]
            ret['marks'] = marks_padded

        return ret

    def _dataset_statistics(self, data):
        statistics = {}
        # get target seqs
        type_seqs = data['types']
        dt_seqs = [[t_seq[i] - t_seq[i - 1] for i in range(1, len(t_seq))]
                   for t_seq in data['timestamps']]

        event_num = self._type_padding

        if isinstance(type_seqs[0], list):
            types = sum(type_seqs, [])
        else:
            types = np.concatenate(type_seqs)
        dts = sum(dt_seqs, [])

        # get statistics
        statistics['n_records'] = len(type_seqs)
        statistics['max_len_of_record'] = max([len(seq) for seq in type_seqs])
        statistics['min_len_of_record'] = min([len(seq) for seq in type_seqs])
        statistics['max_dt'] = np.max(dts)
        statistics['mean_dt'] = np.mean(dts)

        type_count = [0] * event_num
        for t in types:
            type_count[int(t)] += 1

        type_ratio = np.divide(type_count, np.sum(type_count))
        statistics['max_type_ratio'] = np.max(type_ratio)

        # plot marks
        mark_seqs = data.get('marks')
        if mark_seqs is not None:

            import matplotlib.pyplot as plt
            if isinstance(mark_seqs, list):
                marks = sum(mark_seqs, [])
            else:
                marks = np.concatenate(mark_seqs)
            statistics['mean_marks'] = np.mean(marks)
            marks_rmse = np.sqrt(np.mean(np.subtract(marks, np.mean(marks)) ** 2))
            print(marks_rmse)
            # plt.figure()
            # plt.hist(marks, bins=100, range=(0, 2000))
            # plt.show()

        return statistics


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
'''
Given a similarity matrix on sparse or numpy format, it creates a
Laplacian Embedding, for various type of graph Laplacian as well as
normalization.
So far the similarity is assumed to represent a fully connected graph.
'''
# try:
#     np
# except:
import numpy as np
# graphical utilities
import matplotlib.pyplot as plt
import time
# import sys
from scipy import sparse as sp
# from sinkhorn_knopp import sinkhorn_knopp as skp
"""
if 'sinkhorn' not in sys.modules:
    '''
    TODO: may not be true for every-one
    '''
    sys.path.append("/usr/local/lib/python3.5/dist-packages")
"""

def check_matrix(matrix):
    '''
    check that the matrix is square and symmetric and return its
    type.
    '''
    # check for squareness
    (n, m) = matrix.shape
    if n != m:
        raise ValueError('the matrix is not square')
    # check for symmetry
    if sp.issparse(matrix):
        if not (sp.find(matrix)[2] >= 0).all():
            raise ValueError('the matrix is not nonnegative')
        if not (abs(matrix-matrix.T) > 1e-10).nnz == 0:
            raise ValueError('specified similarity matrix is not\
                symmetric')
        return('sparse')
    else:
        if not np.all(matrix >= 0):
            raise ValueError('matrix is not nonnegative')
        if not np.allclose(matrix, matrix.T, atol=1e-6):
            raise ValueError('specified similarity matrix is not\
                symmetric.')
        return('dense')

def normalize(matrix, type_normalization):
    '''
    normalize a similarity matrix with coifman
    or sinkhorn normalization
    '''
    type_matrix = check_matrix(matrix)
    n = matrix.shape[0]
    if type_normalization is None:
        if type_matrix == 'sparse':
            return(sp.csr_matrix(matrix, copy=True, dtype='float64'))
        else:
            return(matrix)

    if type_normalization not in ['coifman', 'sinkhorn']:
        raise ValueError('the normalisation is not admissible')

    if type_normalization == 'coifman' and type_matrix == 'dense':
        W = matrix
        D = np.diag(1./matrix.dot(np.ones(n)))
        W = D.dot(matrix.dot(D))
        return(W)
    if type_normalization == 'coifman' and type_matrix == 'sparse':
        W = sp.csr_matrix(matrix, copy=True, dtype='float64')
        d = W.sum(axis=1).getA1()
        # d = np.array(sim_mat.sum(axis=0, dtype=sim_mat.dtype))[0]
        Dinv = sp.diags(1./d)
        W = Dinv.dot(W.dot(Dinv))
        return(W)
    if type_normalization == 'sinkhorn' and type_matrix == 'dense':
        sk = skp.SinkhornKnopp()
        normalized_matrix = sk.fit(matrix)
        return(normalized_matrix)
    else:
        raise ValueError('There is a problem in the choice of \
        type_normalization')

def make_laplacian_emb_dense(matrix,
                             dim,
                             type_laplacian='unnormalized',
                             type_normalization=None,
                             scaled=False, verbose=0):
        '''
        Create a vector of dimension dim for every node whose pair-wise
        distances are stored in similarity_matrix.
        compute first k eigenvectors of matrix of the similarity_matrix
        L = diag(W 1) - W

        INPUT:
            matrix should be a numpy matrix here.
        '''
        admissible_type_laplacian = ['symmetric',
                                     'random_walk',
                                     'unnormalized']
        if type_laplacian not in admissible_type_laplacian:
            raise ValueError('the parameter type_laplacian is\
             not well specified!')

        n = matrix.shape[0]
        n_vec = dim+1
        # normalize similarity_matrix
        W = normalize(matrix, type_normalization)

        # compute a laplacian
        if type_laplacian == 'random_walk':
            Laplacian = np.eye(n) - W.dot(np.diag(1/W.dot(np.ones(n))))
        elif type_laplacian == 'symmetric':
            Laplacian = np.eye(n) - np.diag(1/W.dot(np.ones(n))**(0.5)).\
                dot(W.dot(np.diag(1/W.dot(np.ones(n))**(0.5))))
        elif type_laplacian == 'unnormalized':
            Laplacian = np.diag(W.dot(np.ones(n))) - W

        # compute embedding
        [d, Vec] = np.linalg.eig(Laplacian)
        idx = d.argsort()
        d = d[idx]
        Vec = Vec[:, idx]
        if scaled:
            if scaled == 'CTD':
                dsqrt = np.sqrt(d[1:n_vec])
            else:
                dsqrt = np.arange(1, n_vec)
            V = Vec[:, 1:n_vec].dot(np.diag(1./dsqrt))
        else:
            V = Vec[:, 1:n_vec]
        '''
        TODO: check if it not supposed to be Vec[:,2:k].dot(np.diag(d[2:k]))
        '''
        return(V)

def make_laplacian_emb_sparse(matrix,
                              dim,
                              type_laplacian='unnormalized',
                              type_normalization=None,
                              scaled=False, verbose=0):
    '''
    Create a vector of dimension dim for every node whose pair-wise
    distances are stored in similarity_matrix.
    compute first k eigenvectors of matrix of the similarity_matrix
    L = diag(W 1) - W
    '''
    admissible_type_laplacian = ['symmetric', 'random_walk', 'unnormalized']
    if type_laplacian not in admissible_type_laplacian:
        raise ValueError('the parameter type_laplacian is not well specified!'
                         '(choose between symmetric, random_walk,'
                         ' or unnormalized)')

    n = matrix.shape[0]
    n_vec = dim + 1
    W = normalize(matrix, type_normalization)
    W.dtype = 'float64'
    # compute a laplacian
    lap = sp.csgraph.laplacian(W, normed=False, return_diag=False)

    if type_laplacian == 'random_walk':
        d = W.sum(axis=1).getA1()
        Dinv = sp.diags(1./d)
        lap = Dinv.dot(lap)
    elif type_laplacian == 'symmetric':
        lap = sp.csgraph.laplacian(W, normed=True, return_diag=False)

    # lap = sp.csr_matrix(lap, copy=True, dtype='float64')
    lap.dtype = 'float64'

    # Compute embedding
    t0 = time.time()
    # Largest eigenvalue first
    (evals_max, _) = sp.linalg.eigsh(lap, 1, which='LA', tol=1e-15, ncv=20)
    t1 = time.time()
    if verbose > 1:
        print('Computed largest eigenvalue of A in {}s.\n'.format(t1-t0))
    maxval = float(evals_max)
    # Then largest eigenvalue of minus laplacian
    # (it is faster to compute largest eigenvalues)
    m_lap = lap
    m_lap *= -1
    m_lap += sp.diags(np.tile(maxval, (n)))

    # m_lap = maxval * sp.identity(n) - lap
    # evec0 = np.ones(n)  # , v0=evec0)
    evals_small, evecs_small = sp.linalg.eigsh(m_lap, n_vec, which='LA', tol=1e-15)
    # eval_s, evec_s = eigsh(lap, n_vec, which='SM', v0=evec0)
    t2 = time.time()
    if verbose > 0:
        print('Computed Laplacian embedding of dim. {} in {}s.\n'.format(
            dim, t2-t1))
    evals_small = maxval - evals_small
    idx = np.array(evals_small).argsort()
    d = evals_small[idx]
    Vec = evecs_small[:, idx]
    V = Vec[:, 1:n_vec]
    if scaled:
        if scaled == 'CTD':
            dsqrt = np.sqrt(d[1:n_vec])
        else:
            dsqrt = np.arange(1, n_vec)
        V = V.dot(np.diag(1./dsqrt))
    return(V)

def make_laplacian_emb(matrix,
                       dim,
                       type_laplacian='unnormalized',
                       type_normalization=None,
                       scaled=False,
                       verb=0):
    type_matrix = check_matrix(matrix)
    if type_matrix == 'dense':
        V = make_laplacian_emb_dense(matrix,
                                     dim,
                                     type_laplacian,
                                     type_normalization,
                                     scaled, verbose=verb)
    else:
        V = make_laplacian_emb_sparse(matrix,
                                      dim,
                                      type_laplacian,
                                      type_normalization,
                                      scaled, verbose=verb)
    return(V)

def vizualize_embedding(embedding, title=None, perm=None):
        embedding = np.real(embedding)
        n = embedding.shape[0]
        plt.figure(int(time.clock()*100))
        if perm is not None:
            plt.scatter(embedding[:, 0], embedding[:, 1], c=perm)
        else:
            plt.scatter(embedding[:, 0], embedding[:, 1], c=np.arange(n))
        if title is not None:
            plt.title(title)

if __name__ == '__main__':
    # Embedding_graph()
    D = np.abs(np.random.normal(0, 1, (5, 5)))
    D = np.transpose(D) + D
    S = sp.csr_matrix(D)

    embeddi = make_laplacian_emb(D,
                                 3,
                                 type_laplacian='unnormalized',
                                 type_normalization='coifman',
                                 scaled=False)

    print(embeddi.shape)
    vizualize_embedding(embeddi)

    embeddi = make_laplacian_emb(S,
                                 3,
                                 type_laplacian='unnormalized',
                                 type_normalization='coifman',
                                 scaled=False)
    print(embeddi.shape)
    vizualize_embedding(embeddi)

reason: import numpy,from scipy


text:
import museval
from tqdm import tqdm

import numpy as np
import torch

import data.utils
import model.utils as model_utils
import utils

def compute_model_output(model, inputs):
    '''
    Computes outputs of model with given inputs. Does NOT allow propagating gradients! See compute_loss for training.
    Procedure depends on whether we have one model for each source or not
    :param model: Model to train with
    :param compute_grad: Whether to compute gradients
    :return: Model outputs, Average loss over batch
    '''
    all_outputs = {}

    if model.separate:
        for inst in model.instruments:
            output = model(inputs, inst)
            all_outputs[inst] = output[inst].detach().clone()
    else:
        all_outputs = model(inputs)

    return all_outputs

def predict(audio, model):
    '''
    Predict sources for a given audio input signal, with a given model. Audio is split into chunks to make predictions on each chunk before they are concatenated.
    :param audio: Audio input tensor, either Pytorch tensor or numpy array
    :param model: Pytorch model
    :return: Source predictions, dictionary with source names as keys
    '''
    if isinstance(audio, torch.Tensor):
        is_cuda = audio.is_cuda()
        audio = audio.detach().cpu().numpy()
        return_mode = "pytorch"
    else:
        return_mode = "numpy"

    expected_outputs = audio.shape[1]

    # Pad input if it is not divisible in length by the frame shift number
    output_shift = model.shapes["output_frames"]
    pad_back = audio.shape[1] % output_shift
    pad_back = 0 if pad_back == 0 else output_shift - pad_back
    if pad_back > 0:
        audio = np.pad(audio, [(0,0), (0, pad_back)], mode="constant", constant_values=0.0)

    target_outputs = audio.shape[1]
    outputs = {key: np.zeros(audio.shape, np.float32) for key in model.instruments}

    # Pad mixture across time at beginning and end so that neural network can make prediction at the beginning and end of signal
    pad_front_context = model.shapes["output_start_frame"]
    pad_back_context = model.shapes["input_frames"] - model.shapes["output_end_frame"]
    audio = np.pad(audio, [(0,0), (pad_front_context, pad_back_context)], mode="constant", constant_values=0.0)

    # Iterate over mixture magnitudes, fetch network prediction
    with torch.no_grad():
        for target_start_pos in range(0, target_outputs, model.shapes["output_frames"]):
            # Prepare mixture excerpt by selecting time interval
            curr_input = audio[:, target_start_pos:target_start_pos + model.shapes["input_frames"]] # Since audio was front-padded input of [targetpos:targetpos+inputframes] actually predicts [targetpos:targetpos+outputframes] target range

            # Convert to Pytorch tensor for model prediction
            curr_input = torch.from_numpy(curr_input).unsqueeze(0)

            # Predict
            for key, curr_targets in compute_model_output(model, curr_input).items():
                outputs[key][:,target_start_pos:target_start_pos+model.shapes["output_frames"]] = curr_targets.squeeze(0).cpu().numpy()

    # Crop to expected length (since we padded to handle the frame shift)
    outputs = {key : outputs[key][:,:expected_outputs] for key in outputs.keys()}

    if return_mode == "pytorch":
        outputs = torch.from_numpy(outputs)
        if is_cuda:
            outputs = outputs.cuda()
    return outputs

def predict_song(args, audio_path, model):
    '''
    Predicts sources for an audio file for which the file path is given, using a given model.
    Takes care of resampling the input audio to the models sampling rate and resampling predictions back to input sampling rate.
    :param args: Options dictionary
    :param audio_path: Path to mixture audio file
    :param model: Pytorch model
    :return: Source estimates given as dictionary with keys as source names
    '''
    model.eval()

    # Load mixture in original sampling rate
    mix_audio, mix_sr = data.utils.load(audio_path, sr=None, mono=False)
    mix_channels = mix_audio.shape[0]
    mix_len = mix_audio.shape[1]

    # Adapt mixture channels to required input channels
    if args.channels == 1:
        mix_audio = np.mean(mix_audio, axis=0, keepdims=True)
    else:
        if mix_channels == 1: # Duplicate channels if input is mono but model is stereo
            mix_audio = np.tile(mix_audio, [args.channels, 1])
        else:
            assert(mix_channels == args.channels)

    # resample to model sampling rate
    mix_audio = data.utils.resample(mix_audio, mix_sr, args.sr)

    sources = predict(mix_audio, model)

    # Resample back to mixture sampling rate in case we had model on different sampling rate
    sources = {key : data.utils.resample(sources[key], args.sr, mix_sr) for key in sources.keys()}

    # In case we had to pad the mixture at the end, or we have a few samples too many due to inconsistent down- and upsamṕling, remove those samples from source prediction now
    for key in sources.keys():
        diff = sources[key].shape[1] - mix_len
        if diff > 0:
            print("WARNING: Cropping " + str(diff) + " samples")
            sources[key] = sources[key][:, :-diff]
        elif diff < 0:
            print("WARNING: Padding output by " + str(diff) + " samples")
            sources[key] = np.pad(sources[key], [(0,0), (0, -diff)], "constant", 0.0)

        # Adapt channels
        if mix_channels > args.channels:
            assert(args.channels == 1)
            # Duplicate mono predictions
            sources[key] = np.tile(sources[key], [mix_channels, 1])
        elif mix_channels < args.channels:
            assert(mix_channels == 1)
            # Reduce model output to mono
            sources[key] = np.mean(sources[key], axis=0, keepdims=True)

        sources[key] = np.asfortranarray(sources[key]) # So librosa does not complain if we want to save it

    return sources

def evaluate(args, dataset, model, instruments):
    '''
    Evaluates a given model on a given dataset
    :param args: Options dict
    :param dataset: Dataset object
    :param model: Pytorch model
    :param instruments: List of source names
    :return: Performance metric dictionary, list with each element describing one dataset sample's results
    '''
    perfs = list()
    model.eval()
    with torch.no_grad():
        for example in dataset:
            print("Evaluating " + example["mix"])

            # Load source references in their original sr and channel number
            target_sources = np.stack([data.utils.load(example[instrument], sr=None, mono=False)[0].T for instrument in instruments])

            # Predict using mixture
            pred_sources = predict_song(args, example["mix"], model)
            pred_sources = np.stack([pred_sources[key].T for key in instruments])

            # Evaluate
            SDR, ISR, SIR, SAR, _ = museval.metrics.bss_eval(target_sources, pred_sources)
            song = {}
            for idx, name in enumerate(instruments):
                song[name] = {"SDR" : SDR[idx], "ISR" : ISR[idx], "SIR" : SIR[idx], "SAR" : SAR[idx]}
            perfs.append(song)

    return perfs

def validate(args, model, criterion, test_data):
    '''
    Iterate with a given model over a given test dataset and compute the desired loss
    :param args: Options dictionary
    :param model: Pytorch model
    :param criterion: Loss function to use (similar to Pytorch criterions)
    :param test_data: Test dataset (Pytorch dataset)
    :return:
    '''
    # PREPARE DATA
    dataloader = torch.utils.data.DataLoader(test_data,
                                             batch_size=args.batch_size,
                                             shuffle=False,
                                             num_workers=args.num_workers)

    # VALIDATE
    model.eval()
    total_loss = 0.
    with tqdm(total=len(test_data) // args.batch_size) as pbar, torch.no_grad():
        for example_num, (x, targets) in enumerate(dataloader):
            if args.cuda:
                x = x.cuda()
                for k in list(targets.keys()):
                    targets[k] = targets[k].cuda()

            _, avg_loss = model_utils.compute_loss(model, x, targets, criterion)

            total_loss += (1. / float(example_num + 1)) * (avg_loss - total_loss)

            pbar.set_description("Current loss: {:.4f}".format(total_loss))
            pbar.update(1)

    return total_loss


text:
#! /usr/bin/python3
# -*- coding: utf-8 -*-
# @Author: Peng Xiang

import numpy as np
import torch.nn.functional as F
from torch import nn
import torch
from torch import nn, einsum
from SiaTrans.utils import MLP_Res, grouping_operation, query_knn

class SkipTransformer(nn.Module):
    def __init__(self, in_channel, dim=256, n_knn=16, pos_hidden_dim=64, attn_hidden_multiplier=4):
        super(SkipTransformer, self).__init__()
        self.mlp_v = MLP_Res(in_dim=in_channel*2,
                             hidden_dim=in_channel, out_dim=in_channel)
        self.n_knn = n_knn
        self.conv_key = nn.Conv1d(in_channel, dim, 1)
        self.conv_query = nn.Conv1d(in_channel, dim, 1)
        self.conv_value = nn.Conv1d(in_channel, dim, 1)

        self.pos_mlp = nn.Sequential(
            nn.Conv2d(3, pos_hidden_dim, 1),
            nn.BatchNorm2d(pos_hidden_dim),
            nn.ReLU(),
            nn.Conv2d(pos_hidden_dim, dim, 1)
        )

        self.attn_mlp = nn.Sequential(
            nn.Conv2d(dim, dim * attn_hidden_multiplier, 1),
            nn.BatchNorm2d(dim * attn_hidden_multiplier),
            nn.ReLU(),
            nn.Conv2d(dim * attn_hidden_multiplier, dim, 1)
        )

        self.conv_end = nn.Conv1d(dim, in_channel, 1)

    def forward(self, pos, key, query, include_self=True):
        """
        Args:
            pos: (B, 3, N)
            key: (B, in_channel, N)
            query: (B, in_channel, N)
            include_self: boolean

        Returns:
            Tensor: (B, in_channel, N), shape context feature
        """
        value = self.mlp_v(torch.cat([key, query], 1))
        identity = value
        key = self.conv_key(key)
        query = self.conv_query(query)
        value = self.conv_value(value)
        b, dim, n = value.shape

        pos_flipped = pos.permute(0, 2, 1).contiguous()
        idx_knn = query_knn(self.n_knn, pos_flipped,
                            pos_flipped, include_self=include_self)

        key = grouping_operation(key, idx_knn)  # b, dim, n, n_knn
        qk_rel = query.reshape((b, -1, n, 1)) - key

        pos_rel = pos.reshape((b, -1, n, 1)) - \
            grouping_operation(pos, idx_knn)  # b, 3, n, n_knn
        pos_embedding = self.pos_mlp(pos_rel)

        attention = self.attn_mlp(qk_rel + pos_embedding)  # b, dim, n, n_knn
        attention = torch.softmax(attention, -1)

        value = value.reshape((b, -1, n, 1)) + pos_embedding  #

        agg = einsum('b c i j, b c i j -> b c i',
                     attention, value)  # b, dim, n
        y = self.conv_end(agg)

        return y + identity

class Self_Attn(nn.Module):
    """ Self attention Layer"""

    def __init__(self, in_channel):
        super(Self_Attn, self).__init__()

        self.query_conv = nn.Conv1d(
            in_channels=in_channel, out_channels=in_channel//8, kernel_size=1)
        self.key_conv = nn.Conv1d(
            in_channels=in_channel, out_channels=in_channel//8, kernel_size=1)
        self.value_conv = nn.Conv1d(
            in_channels=in_channel, out_channels=in_channel, kernel_size=1)
        self.gamma = nn.Parameter(torch.zeros(1))

        self.softmax = nn.Softmax(dim=-1)

    def forward(self, x):
        """
            inputs :
                x : input feature maps( B X C X n_points)
            returns :
                out : self attention value + input feature
                attention: B X N X N (N is Width*Height)
        """
        m_batchsize, C, npoints = x.size()
        proj_query = self.query_conv(x).view(
            m_batchsize, -1, npoints).permute(0, 2, 1)  # B X CX(N)
        # print(proj_query.size())
        proj_key = self.key_conv(x).view(
            m_batchsize, -1, npoints)  # B X C x (*W*H)
        # print(proj_key.size())
        energy = torch.bmm(proj_query, proj_key)  # transpose check
        attention = self.softmax(energy)  # BX (N) X (N)
        proj_value = self.value_conv(x).view(
            m_batchsize, -1, npoints)  # B X C X N
        # print(proj_value.size())

        out = torch.bmm(proj_value, attention.permute(0, 2, 1))
        out = out.view(m_batchsize, C, npoints)

        out = self.gamma*out + x
        return out

# bs = 8
# input = torch.randn((bs, 3, 512))
# feat_global = torch.randn((bs, 32, 512))
# K_rev = torch.randn((bs, 128, 512))

# ps = nn.ConvTranspose1d(32, 128, 4, 4, bias=False)
# res = ps(feat_global)

# up_sampler = nn.Upsample(scale_factor=4)
# H = torch.randn((bs, 128, 512))
# H_up = up_sampler(H)
# print(res.size())
# print(H_up.size())

class PointNonLocalCell(nn.Module):
    """Input
        feature: (batch_size, ndataset, channel) Torch tensor
        new_point: (batch_size, npoint, nsample, channel)
        Output
        (batch_size, npoint, nsample, channel)
    """

    def __init__(self, feature_channel, in_channel, mlps):
        super(PointNonLocalCell, self).__init__()
        self.bottleneck_channel = mlps[0]
        self.transformed_feature = nn.Conv2d(
            feature_channel, self.bottleneck_channel * 2, 1)
        self.tf_bn = nn.BatchNorm2d(self.bottleneck_channel*2)
        self.transformed_new_point = nn.Conv2d(
            in_channel, self.bottleneck_channel, 1)
        self.tnp_bn = nn.BatchNorm2d(self.bottleneck_channel)
        self.new_nonlocal_point = nn.Conv2d(
            self.bottleneck_channel, mlps[-1], 1)
        self.nnp_bn = nn.BatchNorm2d(mlps[-1])

    def forward(self, feature, new_point):
        B, P, C, S = new_point.shape
        FB, FC, FD = feature.shape
        feature = feature.view(FB, FC, 1, FD).permute(0, 1, 3, 2)
        new_point = new_point.permute(0, 2, 1, 3)

        # B, P, S, C = new_point.shape
        # FB, FD, FC = feature.shape
        # feature = feature.view(FB, FD, 1, FC).permute(0,3,1,2) FB, FC, FD, 1
        # new_point = new_point.permute(0,3,1,2) B, C, P, S

        transformed_feature = self.tf_bn(
            self.transformed_feature(feature)).permute(0, 2, 3, 1)
        transformed_new_point = self.tnp_bn(
            self.transformed_new_point(new_point)).permute(0, 2, 3, 1)
        transformed_new_point = torch.reshape(
            transformed_new_point, (B, P*S, self.bottleneck_channel))
        transformed_feature1 = torch.squeeze(
            transformed_feature[:, :, :, :self.bottleneck_channel], 2)
        transformed_feature2 = torch.squeeze(
            transformed_feature[:, :, :, self.bottleneck_channel:], 2)

        attention_map = torch.matmul(
            transformed_new_point, transformed_feature1.transpose(1, 2))  # mode = 'dot'
        attention_map = attention_map / np.sqrt(self.bottleneck_channel)
        attention_map = F.softmax(attention_map, -1)

        new_nonlocal_point = torch.matmul(attention_map, transformed_feature2)
        new_nonlocal_point = torch.reshape(
            new_nonlocal_point, (B, P, S, self.bottleneck_channel)).permute(0, 3, 1, 2)
        new_nonlocal_point = self.nnp_bn(self.new_nonlocal_point(
            new_nonlocal_point)).permute(0, 2, 1, 3)
        # print(new_nonlocal_point.size())
        new_nonlocal_point = torch.squeeze(new_nonlocal_point, 1)

        return new_nonlocal_point

# bs = 16
# new_feature = torch.randn((bs, 128, 512))
# new_points = torch.randn((bs, 1, 128, 512))
# mlp = [128, 128, 256]
# PNL = PointNonLocalCell(128, 128, [max(32, 128+3//2), 128])
# new_nonlocal_point = PNL(new_feature, new_points)
# print(new_nonlocal_point.size())


text:
#%pylab inline
import pydotplus
from sklearn.utils.random import check_random_state
from sklearn.datasets import load_diabetes
from gplearn.genetic import SymbolicRegressor
import numpy as np
from sklearn.model_selection import train_test_split
from gplearn.complexity import complexity

rng = check_random_state(2)

(X, y) = load_diabetes(return_X_y=True)
print("diabtes dataset")
print("size:", X.shape)
ground_compl = complexity(X, y)
print("Ground complexity:", sum(ground_compl))

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.70, random_state=42)

ground_compl = complexity(X_train, y_train)
print("Ground complexity on training data:", sum(ground_compl))

est_gp = SymbolicRegressor(population_size=600,
                           generations=100, stopping_criteria=0.01,
                           p_crossover=0.8, p_subtree_mutation=0.1,
                           p_hoist_mutation=0, p_point_mutation=0,
                           parsimony_coefficient=0,
                            verbose=1,
                            random_state=42,
                           n_jobs=2,
                           safe_best_program_to_file=False,
                           tournament_size=10,
                           first_tournament="fitness",
                           #second_tournament="complexity",
                           second_tournament_size=1.4)
print("Run GP with parameters: ", est_gp.get_params())
est_gp.fit(X_train, y_train)

program = est_gp._program
graph = pydotplus.graphviz.graph_from_dot_data(program.export_graphviz())
filename = "final_best_program.pdf"
graph.write_pdf(filename)

y_pred = est_gp.predict(X_test)

print("test fitness:", np.average(np.abs(y_pred - y_test)))

score_gp = est_gp.score(X_test, y_test)
print(score_gp)

#print("final program complexity:", program.complexity_)

#print(est_gp._program.parents)
#
#idx = est_gp._program.parents['donor_idx']
#fade_nodes = est_gp._program.parents['donor_nodes']
#print(est_gp._programs[-2][idx])
#print('Fitness:', est_gp._programs[-2][idx].fitness_)


text:
import torch
import random
import numpy as np
import torch.nn as nn
import torch.nn.functional as F
import torchvision
from auto_LiRPA import BoundedModule, BoundedTensor
from auto_LiRPA.perturbations import *
import sys
sys.path.append('../examples/vision')
import models
from testcase import TestCase

class TestResnetPatches(TestCase):
    def __init__(self, methodName='runTest', generate=False):
        super().__init__(methodName,
            seed=1234, ref_path='data/resnet_patches_test_data',
            generate=generate)

    def test(self):
        model_oris = [
            models.model_resnet(width=1, mult=2),
            models.ResNet18(in_planes=2)
        ]
        self.result = []

        for model_ori in model_oris:
            conv_mode = 'patches' # conv_mode can be set as 'matrix' or 'patches'

            normalize = torchvision.transforms.Normalize(mean = [0.4914, 0.4822, 0.4465], std = [0.2023, 0.1994, 0.2010])
            test_data = torchvision.datasets.CIFAR10("./data", train=False, download=True,
                            transform=torchvision.transforms.Compose([torchvision.transforms.ToTensor(), normalize]))
            N = 1
            n_classes = 10

            image = torch.Tensor(test_data.data[:N]).reshape(N,3,32,32)
            image = image.to(torch.float32) / 255.0

            model = BoundedModule(model_ori, image, bound_opts={"conv_mode": conv_mode})

            ptb = PerturbationLpNorm(norm = np.inf, eps = 0.03)
            image = BoundedTensor(image, ptb)
            pred = model(image)
            lb, ub = model.compute_bounds(IBP=False, C=None, method='backward')
            self.result += [lb, ub]

        self.check()

if __name__ == '__main__':
    # Change to generate=True when genearting reference results
    testcase = TestResnetPatches(generate=True)
    testcase.test()


text:
#!/usr/bin/env python

#
# Test BLAS 1
#

from numpy.testing import assert_almost_equal
import numpy as np
import nipy.neurospin.bindings as fb

n = 15

def test_dnrm2():
    x = np.random.rand(n)
    assert_almost_equal(np.sqrt(np.sum(x**2)), fb.blas_dnrm2(x))

def test_dasum():
    x = np.random.rand(n)
    assert_almost_equal(np.sum(np.abs(x)), fb.blas_dasum(x))

def test_ddot():
    x = np.random.rand(n)
    y = np.random.rand(n)
    assert_almost_equal(np.dot(x,y), fb.blas_ddot(x, y))

def test_daxpy():
    x = np.random.rand(n)
    y = np.random.rand(n)
    alpha = np.random.rand()
    assert_almost_equal(alpha*x+y, fb.blas_daxpy(alpha, x, y))

def test_dscal():
    x = np.random.rand(n)
    alpha = np.random.rand()
    assert_almost_equal(alpha*x, fb.blas_dscal(alpha, x))

if __name__ == "__main__":
    import nose
    nose.run(argv=['', __file__])

max_stars_repo_name: fperez/nipy

max_issues_repo_name: fperez/nipy

max_forks_repo_name: fperez/nipy

reason: import numpy,from numpy


text:
#!/usr/bin/env python
"""
# Author: ******
# Created Time : Mon 23 Apr 2021 08:26:32
# File Name: integrated_feat.py
# Description:`

"""
"""
test:
import scanpy as sc
from integrated_feat import image_feature, image_crop
adata = sc.read(data_path)
image_crop(adata, save_path)
adata = image_feature(adata).extract_image_feat()

"""
import os
import math
import anndata
import numpy as np
import scanpy as sc
import pandas as pd
from PIL import Image
from pathlib import Path
from scipy.sparse import issparse, isspmatrix_csr, csr_matrix, spmatrix
from sklearn.metrics import pairwise_distances
from sklearn.linear_model import LinearRegression
from sklearn.decomposition import PCA
from tqdm import tqdm

import torch
import torch.nn
import torchvision.models as models
from torch.autograd import Variable
import torchvision.transforms as transforms

class image_feature:
    def __init__(self,
                 adata,
                 pca_components=50,
                 cnnType='ResNet50',
                 verbose=False,
                 seeds=88,
                ):
        self.device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
        self.adata = adata
        self.pca_components = pca_components
        self.verbose = verbose
        self.seeds = seeds
        self.cnnType = cnnType

    def load_cnn_model(self):

        if self.cnnType == 'ResNet50':
            cnn_pretrained_model = models.resnet50(pretrained=True)
            cnn_pretrained_model.to(self.device)
        elif self.cnnType == 'Resnet152':
            cnn_pretrained_model = models.resnet152(pretrained=True)
            cnn_pretrained_model.to(self.device)
        elif self.cnnType == 'Vgg19':
            cnn_pretrained_model = models.vgg19(pretrained=True)
            cnn_pretrained_model.to(self.device)
        elif self.cnnType == 'Vgg16':
            cnn_pretrained_model = models.vgg16(pretrained=True)
            cnn_pretrained_model.to(self.device)
        elif self.cnnType == 'DenseNet121':
            cnn_pretrained_model = models.densenet121(pretrained=True)
            cnn_pretrained_model.to(self.device)
        elif self.cnnType == 'Inception_v3':
            cnn_pretrained_model = models.inception_v3(pretrained=True)
            cnn_pretrained_model.to(self.device)
        else:
            raise ValueError(
                    f"""\
                        {self.cnnType} is not a valid type.
                        """)
        return cnn_pretrained_model

    def extract_image_feat(self,):

        transform_list = [transforms.ToTensor(),
                          transforms.Normalize(mean=[0.485, 0.456, 0.406],
                          std =[0.229, 0.224, 0.225])]
        # transform_list = [transforms.ToTensor(),
        #                   transforms.Normalize(mean=[0.54, 0.51, 0.68],
        #                   std =[0.25, 0.21, 0.16])]
        img_to_tensor = transforms.Compose(transform_list)

        feat_df = pd.DataFrame()
        model = self.load_cnn_model()
        #model.fc = torch.nn.LeakyReLU(0.1)
        model.eval()

        if "slices_path" not in self.adata.obs.keys():
             raise ValueError("Please run the function image_crop first")

        with tqdm(total=len(self.adata),
              desc="Extract image feature",
              bar_format="{l_bar}{bar} [ time left: {remaining} ]",) as pbar:
            for spot, slice_path in self.adata.obs['slices_path'].items():
                spot_slice = Image.open(slice_path)
                spot_slice = spot_slice.resize((224,224))
                spot_slice = np.asarray(spot_slice, dtype="int32")
                spot_slice = spot_slice.astype(np.float32)
                tensor = img_to_tensor(spot_slice)
                tensor = tensor.resize_(1,3,224,224)
                tensor = tensor.to(self.device)
                result = model(Variable(tensor))
                result_npy = result.data.cpu().numpy().ravel()
                feat_df[spot] = result_npy
                feat_df = feat_df.copy()
                pbar.update(1)
        self.adata.obsm["image_feat"] = feat_df.transpose().to_numpy()
        if self.verbose:
            print("The image feature is added to adata.obsm['image_feat'] !")
        pca = PCA(n_components=self.pca_components, random_state=self.seeds)
        pca.fit(feat_df.transpose().to_numpy())
        self.adata.obsm["image_feat_pca"] = pca.transform(feat_df.transpose().to_numpy())
        if self.verbose:
            print("The pca result of image feature is added to adata.obsm['image_feat_pca'] !")
        del model
        return self.adata

def image_crop(
                adata,
                save_image_path,
                library_id=None,
                crop_size=50,
                target_size=224,
                verbose=False,):
    if library_id is None:
       library_id = list(adata.uns["spatial"].keys())[0]

    image = adata.uns["spatial"][library_id]["images"][
            adata.uns["spatial"][library_id]["use_quality"]]
    if image.dtype == np.float32 or image.dtype == np.float64:
        image = (image * 255).astype(np.uint8)
    img_pillow = Image.fromarray(image)
    tile_names = []

    with tqdm(total=len(adata),
              desc="Tiling image",
              bar_format="{l_bar}{bar} [ time left: {remaining} ]") as pbar:
        for imagerow, imagecol in zip(adata.obs["imagerow"], adata.obs["imagecol"]):
            imagerow_down = imagerow - crop_size / 2
            imagerow_up = imagerow + crop_size / 2
            imagecol_left = imagecol - crop_size / 2
            imagecol_right = imagecol + crop_size / 2
            tile = img_pillow.crop(
                (imagecol_left, imagerow_down, imagecol_right, imagerow_up))
            tile.thumbnail((target_size, target_size), Image.ANTIALIAS) #####
            tile.resize((target_size, target_size)) ######
            tile_name = str(imagecol) + "-" + str(imagerow) + "-" + str(crop_size)
            out_tile = Path(save_image_path) / (tile_name + ".png")
            tile_names.append(str(out_tile))
            if verbose:
                print(
                    "generate tile at location ({}, {})".format(
                        str(imagecol), str(imagerow)))
            tile.save(out_tile, "PNG")
            pbar.update(1)

    adata.obs["slices_path"] = tile_names
    if verbose:
        print("The slice path of image feature is added to adata.obs['slices_path'] !")
    return adata

max_stars_repo_name: spatial-Transcriptomics/DeepST

max_issues_repo_name: spatial-Transcriptomics/DeepST

max_forks_repo_name: spatial-Transcriptomics/DeepST


text:
import matplotlib
matplotlib.use('Agg')
from Swing.util.BoxPlot import BoxPlot
from matplotlib.backends.backend_pdf import PdfPages
from scipy import stats

import pdb

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import sys
import os
import time
from Swing.util.mplstyle import style1
import seaborn as sns
from palettable.colorbrewer.qualitative import Set1_3

def get_df(df, fp, min_lag, max_lag, td_window, inftype = "RandomForest"):
    new_df = df[(df['file_path'] == fp) & (df['min_lag'] == min_lag) & (df['max_lag'] == max_lag) & (df['td_window'] == td_window) & (df['InfType'] == inftype)]
    return(new_df)

def load_data():
    input_folder_list = ["/projects/p20519/roller_output/gnw/RandomForest/"]
    agg_df_RF = read_tdr_results(input_folder_list, folder_str = "2017-09")
    agg_df_RF['InfType'] = 'RandomForest'

    input_folder_list = ["/projects/p20519/roller_output/gnw/Dionesus/"]
    agg_df_P = read_tdr_results(input_folder_list, folder_str = "2017-09")
    agg_df_P['InfType'] = 'PLSR'

    input_folder_list = ["/projects/p20519/roller_output/gnw/Lasso/"]
    agg_df_L = read_tdr_results(input_folder_list, folder_str = "2017-09")
    agg_df_L['InfType'] = 'Lasso'

    all_dfs = [agg_df_RF, agg_df_P, agg_df_L]
    merged_df = pd.concat(all_dfs)
    return(merged_df)

def read_tdr_results(folder_list, folder_str):
    agg_df = pd.DataFrame()
    for input_folder in folder_list:
        for file_path in os.listdir(input_folder):
            if folder_str in file_path:
              df = pd.read_csv(input_folder+file_path,sep='\t', engine='python')
              # check if the columns are misaligned.
              if type(df['permutation_n'].iloc[0]) is str:
                  new_col = df.columns.tolist()
                  new_col.pop(0)
                  new_df = df.iloc[:,0:len(df.iloc[0])-1]
                  new_df.columns = new_col
                  df=new_df

              agg_df = agg_df.append(df)
    return(agg_df)

def get_inf_df(network_1, inf_type):
    RFnet1 = network_1[network_1['InfType'] == inf_type]
    RFn1 = RFnet1.groupby('td_window').mean()
    return(RFn1)

def get_comparisons(merged_df, inftypes, window_sizes, network_list):
    overall_df = pd.DataFrame()
    network_1_df = pd.DataFrame()
    for inftype in inftypes:
        for td_window in window_sizes:
            for network in network_list:
                baseline = get_df(merged_df, network, 0, 0, 21, inftype = inftype)
                if len(baseline) == 0:
                    continue
                if 21-td_window > 2:
                    max_lag = 3
                else:
                    max_lag = 21-td_window
                if (td_window == 21):
                    min_lag = 0
                    max_lag = 0
                else:
                    min_lag = 1

                comparisons = get_df(merged_df, network, min_lag, max_lag, td_window, inftype = inftype)
                if len(comparisons) == 0:
                    continue

                # for each statistic, get the percent difference to the baseline comparison.
                stat = 'aupr'
                baseline_mean=baseline[stat].mean()
                comparisons['percent_{}'.format(stat)] = ((comparisons[stat]-baseline_mean)/baseline_mean)*100
                stat = 'auroc'
                baseline_mean=baseline[stat].mean()
                comparisons['percent_{}'.format(stat)] = ((comparisons[stat]-baseline_mean)/baseline_mean)*100
                overall_df = overall_df.append(comparisons.iloc[0:20,:], ignore_index = True)
                if network == network_list[6]:
                    network_1_df = network_1_df.append(comparisons.iloc[0:20,:], ignore_index = True)
                print(comparisons,len(comparisons))
    return(overall_df, network_1_df)

test_statistic = ['aupr', 'auroc']
save_tag = "window_scan"
n_trials = 100

#merged_df = load_data()
#merged_df.to_pickle("merged_window_scan.pkl")

#merged_df = pd.read_pickle("merged_window_scan.pkl")

#network_list = merged_df['file_path'].unique().tolist()
#window_sizes = range(2,22)
#inftypes = ['RandomForest', 'Lasso', 'PLSR']

#overall_df, network_1 = get_comparisons(merged_df, inftypes, window_sizes, network_list)

#overall_df.to_pickle("merged_window_scan_comparisons.pkl")
#network_1.to_pickle("merged_window_scan_comparisons_network1.pkl")
overall_df = pd.read_pickle("merged_window_scan_comparisons.pkl")
network_1 = pd.read_pickle("merged_window_scan_comparisons_network1.pkl")

stat = 'percent_aupr'

meanpointprops = dict(marker = 'o', markeredgecolor='black', markerfacecolor='white', markeredgewidth = 1, markersize = 10)
medianprops = dict(color='black', linewidth=1.5)
g = sns.FacetGrid(overall_df, col = "InfType", size = 12, aspect=0.6)
g = g.map(sns.boxplot, 'td_window', 'percent_aupr', showfliers = False, color = 'cornflowerblue', showmeans=True, meanprops = meanpointprops, medianprops=medianprops)
g.set(ylim=(-60, 200))
g.set(xlabel="Window Size")
g.fig.get_axes()[0].set_ylabel("% change AUPR")
g.fig.get_axes()[0].axhline(y = 0, c = "darkgrey")
g.fig.get_axes()[1].axhline(y = 0, c = "darkgrey")
g.fig.get_axes()[2].axhline(y = 0, c = "darkgrey")

## Plot the AUPR for 1 graph, dionesus
RFn1 = get_inf_df(network_1, inf_type = 'RandomForest')
net1color = 'firebrick'
g.fig.get_axes()[0].plot(RFn1[stat].values, marker='.', color = net1color)

Ln1 = get_inf_df(network_1, inf_type = 'Lasso')
g.fig.get_axes()[1].plot(Ln1[stat].values, marker='.', color = net1color)

Pn1 = get_inf_df(network_1, inf_type = 'PLSR')
g.fig.get_axes()[2].plot(Pn1[stat].values, marker='.', color = net1color)

for axes in g.axes.flat:
    axes.set_xticklabels(['{:d}'.format(x+2) for x in axes.get_xticks()])
g.savefig('combined_10_AUPR_window_scan.pdf')

## Plot AUROC

stat = 'percent_auroc'

meanpointprops = dict(marker = 'o', markeredgecolor='black', markerfacecolor='white', markeredgewidth = 1, markersize = 10)
medianprops = dict(color='black', linewidth=1.5)
g = sns.FacetGrid(overall_df, col = "InfType", size = 12, aspect=0.6)
g = g.map(sns.boxplot, 'td_window', 'percent_auroc', showfliers = False, color = 'darkorange', showmeans=True, meanprops = meanpointprops, medianprops=medianprops)
g.set(ylim=(-50, 60))
g.set(xlabel="Window Size")
g.fig.get_axes()[0].set_ylabel("% change AUROC")
g.fig.get_axes()[0].axhline(y = 0, c = "darkgrey")
g.fig.get_axes()[1].axhline(y = 0, c = "darkgrey")
g.fig.get_axes()[2].axhline(y = 0, c = "darkgrey")

## Plot the AUROC for 1 graph, all inf methods
RFn1 = get_inf_df(network_1, inf_type = 'RandomForest')
net1color = 'firebrick'
g.fig.get_axes()[0].plot(RFn1[stat].values, marker='.', color = net1color)

Ln1 = get_inf_df(network_1, inf_type = 'Lasso')
g.fig.get_axes()[1].plot(Ln1[stat].values, marker='.', color = net1color)

Pn1 = get_inf_df(network_1, inf_type = 'PLSR')
g.fig.get_axes()[2].plot(Pn1[stat].values, marker='.', color = net1color)

for axes in g.axes.flat:
    axes.set_xticklabels(['{:d}'.format(x+2) for x in axes.get_xticks()])

g.savefig('combined_10_AUROC_window_scan.pdf')

max_stars_repo_name: jiawu/Roller

max_issues_repo_name: jiawu/Roller

max_forks_repo_name: jiawu/Roller


text:
import os
import time
from PIL import Image
import numpy as np
import tensorflow as tf
import tensorflow_hub as hub
import matplotlib.pyplot as plt

os.environ['TFHUB_DOWNLOAD'] = 'True'

IMAGE_PATH = 'original.png'
SAVED_MODEL_PATH = "https://tfhub.dev/captain-pool/esrgan-tf2/1"

def preprocess_image(image_path):
    hr_image = tf.image.decode_image(tf.io.read_file(image_path))

    if hr_image.shape[-1] == 4:
        hr_image = hr_image[...,: -1]

    hr_size = (tf.convert_to_tensor(hr_image.shape[:-1]) // 4) * 4
    hr_image = tf.image.crop_to_bounding_box(hr_image, 0, 0, hr_size[0], hr_size[1])
    hr_image = tf.cast(hr_image, tf.float32)

    return tf.expand_dims(hr_image, 0)

def save_image(image, filename):
    if not isinstance(image, Image.Image):
        image = tf.clip_by_value(image, 0, 255)
        image = Image.fromarray(tf.cast(image, tf.uint8).numpy())

    image.save(f"{filename}.jpg")
    print(f"Saved as {filename}.jpg")

hr_image = preprocess_image(IMAGE_PATH)
model = hub.load(SAVED_MODEL_PATH)

sr_image = model(hr_image)
sr_image = tf.squeeze(sr_image)

save_image(tf.squeeze(sr_image), filename='sr_image')

max_stars_repo_name: amalrajan/super-resolution

max_issues_repo_name: amalrajan/super-resolution

max_forks_repo_name: amalrajan/super-resolution