text:
########################################################################
#
#  Copyright 2014 Johns Hopkins University
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Contact: turbulence@pha.jhu.edu
# Website: http://turbulence.pha.jhu.edu/
#
########################################################################

import os
import sys
import numpy as np
import ctypes
import inspect

import pyJHTDB
from pyJHTDB.dbinfo import interpolation_code

class ThresholdInfo(ctypes.Structure):
    _fields_ = [('x', ctypes.c_int),
                ('y', ctypes.c_int),
                ('z', ctypes.c_int),
                ('value', ctypes.c_float)]

class libJHTDB(object):
    def __init__(self,
            auth_token = pyJHTDB.auth_token):
        self.libname = 'libJHTDB'
        lib_location = os.path.dirname(inspect.getfile(pyJHTDB))
        self.lib = np.ctypeslib.load_library(
            self.libname,
            os.path.abspath(os.path.join(lib_location, os.path.pardir)))
        self.authToken = ctypes.c_char_p(auth_token.encode('ascii'))
        self.connection_on = False
        self.hdf5_file_list = []
        self.hdf5_file_desc = {}
        return None
    def initialize(self, exit_on_error = True):
        #initialize gSOAP
        self.lib.soapinit()
        if exit_on_error:
            #enable exit on error
            self.lib.turblibSetExitOnError(ctypes.c_int(1))
        self.connection_on = True
        return None
    def finalize(self):
        #free gSOAP resources
        self.lib.soapdestroy()
        self.connection_on = False
        return None
    def add_hdf5_file(self, filename):
        if pyJHTDB.found_h5py and (not filename in self.hdf5_file_list):
            self.hdf5_file_list.append(filename)
            data = pyJHTDB.h5py.File(filename + '.h5', mode = 'r')
            self.hdf5_file_desc[filename] = {}
            for key in ['_contents', '_dataset', '_size', '_start']:
                self.hdf5_file_desc[filename][key] = data[key][:]
            data.close()
            return self.lib.turblibAddLocalSource(ctypes.c_char_p((filename + '.h5').encode('ascii')))
        else:
            return 0
    def getData(self,
            time, point_coords,
            sinterp = 0, tinterp = 0,
            data_set = 'isotropic1024coarse',
            getFunction = 'getVelocity',
            make_modulo = False):
        if not self.connection_on:
            print('you didn\'t connect to the database')
            sys.exit()
        if not (point_coords.shape[-1] == 3):
            print ('wrong number of values for coordinates in getData')
            sys.exit()
            return None
        if not (point_coords.dtype == np.float32):
            print('point coordinates in getData must be floats. stopping.')
            sys.exit()
            return None
        if (type(sinterp) == str):
            sinterp = interpolation_code[sinterp]
        if (type(tinterp) == str):
            tinterp = interpolation_code[tinterp]
        npoints = point_coords.shape[0]
        for i in range(1, len(point_coords.shape)-1):
            npoints *= point_coords.shape[i]
        if make_modulo:
            pcoords = np.zeros(point_coords.shape, np.float64)
            pcoords[:] = point_coords
            np.mod(pcoords, 2*np.pi, point_coords)
        if not getFunction[0:3] == 'get':
            getFunction = 'get' + getFunction
        get_data = getattr(self.lib, getFunction)
        if getFunction in ['getVelocity',
                           'getForce',
                           'getMagneticField',
                           'getMagneticFieldDebug',
                           'getBunit',
                           'getVectorPotential',
                           'getPressureGradient',
                           'getVelocityLaplacian',
                           'getMagneticFieldLaplacian',
                           'getVectorPotentialLaplacian',
                           'getPressureGradient']:
            result_dim = 3
        elif getFunction in ['getVelocityAndPressure']:
            result_dim = 4
        elif getFunction in ['getPressureHessian']:
            result_dim = 6
        elif getFunction in ['getVelocityGradient',
                             'getMagneticFieldGradient',
                             'getVectorPotentialGradient']:
            result_dim = 9
        elif getFunction in ['getVelocityHessian',
                             'getMagneticFieldHessian',
                             'getVectorPotentialHessian']:
            result_dim = 18
        else:
            print(('wrong result type requested in getData\n'
                 + 'maybe it\'s just missing from the list?'))
            sys.exit()
            return None
        newshape = list(point_coords.shape[0:len(point_coords.shape)-1])
        newshape.append(result_dim)
        result_array = np.empty(newshape, dtype=np.float32)
        get_data(self.authToken,
                 ctypes.c_char_p(data_set.encode('ascii')),
                 ctypes.c_float(time),
                 ctypes.c_int(sinterp), ctypes.c_int(tinterp), ctypes.c_int(npoints),
                 point_coords.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))),
                 result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))))
        return result_array
    def getRawData(
            self,
            time,
            start = np.array([0, 0, 0], dtype = np.int),
            size  = np.array([8, 8, 8], dtype = np.int),
            data_set = 'channel',
            getFunction = 'Velocity'):
        if not self.connection_on:
            print('you didn\'t connect to the database')
            sys.exit()
        if getFunction in ['Velocity',
                           'MagneticField',
                           'VectorPotential']:
            result_dim = 3
        elif getFunction in ['Pressure']:
            result_dim = 1
        else:
            print(('wrong result type requested in getRawData\n'
                 + 'maybe it\'s just missing from the list?'))
            sys.exit()
            return None
        getFunction = 'getRaw' + getFunction
        get_data = getattr(self.lib, getFunction)
        result_array = np.empty(tuple(list(size[::-1]) + [result_dim]), dtype=np.float32)
        get_data(self.authToken,
                 ctypes.c_char_p(data_set.encode('ascii')),
                 ctypes.c_float(time),
                 ctypes.c_int(start[0]),
                 ctypes.c_int(start[1]),
                 ctypes.c_int(start[2]),
                 ctypes.c_int(size[0]),
                 ctypes.c_int(size[1]),
                 ctypes.c_int(size[2]),
                 result_array.ctypes.data_as(ctypes.POINTER(ctypes.c_char)))
        return result_array
    def getBoxFilter(self,
            time, point_coords,
            data_set = 'isotropic1024coarse',
            make_modulo = False,
            field = 'velocity',
            filter_width = 7*2*np.pi / 1024):
        if not self.connection_on:
            print('you didn\'t connect to the database')
            sys.exit()
        if not (point_coords.shape[-1] == 3):
            print ('wrong number of values for coordinates in getBoxFilter')
            sys.exit()
            return None
        if not (point_coords.dtype == np.float32):
            print('point coordinates in getBoxFilter must be floats. stopping.')
            sys.exit()
            return None
        npoints = point_coords.shape[0]
        for i in range(1, len(point_coords.shape)-1):
            npoints *= point_coords.shape[i]
        if make_modulo:
            pcoords = np.zeros(point_coords.shape, np.float64)
            pcoords[:] = point_coords
            np.mod(pcoords, 2*np.pi, point_coords)
        result_array = point_coords.copy()
        self.lib.getBoxFilter(self.authToken,
                 ctypes.c_char_p(data_set.encode('ascii')),
                 ctypes.c_char_p(field.encode('ascii')),
                 ctypes.c_float(time),
                 ctypes.c_float(filter_width),
                 ctypes.c_int(npoints),
                 point_coords.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))),
                 result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))))
        return result_array
    def getThreshold(
            self,
            data_set = 'isotropic1024coarse',
            field = 'vorticity',
            time = 0.1,
            threshold = 0.0,
            cx = 0, cy = 0, cz = 0,
            nx = 4, ny = 4, nz = 4,
            sinterp = 40,
            tinterp = 0):
        result = ctypes.POINTER(ThresholdInfo)()
        result_size = ctypes.c_int()
        call_result = self.lib.getThreshold(
                self.authToken,
                ctypes.c_char_p(data_set.encode('ascii')),
                ctypes.c_char_p(field.encode('ascii')),
                ctypes.c_float(time),
                ctypes.c_float(threshold),
                ctypes.c_int32(sinterp),
                ctypes.c_int32(cx),
                ctypes.c_int32(cy),
                ctypes.c_int32(cz),
                ctypes.c_int32(nx),
                ctypes.c_int32(ny),
                ctypes.c_int32(nz),
                ctypes.byref(result),
                ctypes.byref(result_size))
        data_type = np.dtype([('x', np.int32),
                              ('y', np.int32),
                              ('z', np.int32),
                              ('value', np.float32)])
        data = np.empty((result_size.value,), data_type)
        for i in range(result_size.value):
            data[i]['x'] = result[i].x
            data[i]['y'] = result[i].y
            data[i]['z'] = result[i].z
            data[i]['value'] = result[i].value
        return data
    def getPosition(self,
            starttime = 0.0,
            endtime = 0.1,
            dt = 0.0004,
            point_coords = None,
            sinterp = 4,
            steps_to_keep = 1,
            data_set = 'isotropic1024coarse'):
        if not self.connection_on:
            print('you didn\'t connect to the database')
            sys.exit()
        if not (point_coords.shape[-1] == 3):
            print ('wrong number of values for coordinates in getPosition')
            sys.exit()
            return None
        npoints = point_coords.shape[0]
        result_array = np.empty((npoints, 3), dtype=np.float32)
        traj_array = np.empty((steps_to_keep+1, npoints, 3), dtype = np.float32)
        traj_array[0] = point_coords
        time_array = np.linspace(starttime, endtime, steps_to_keep+1).astype(np.float32)
        time_array[0] = starttime
        evolver = self.lib.getPosition
        if data_set == 'custom':
            evolver = self.lib.getCustomPosition
        print('starting integration loop, dataset is ', data_set)
        for tstep in range(1, steps_to_keep + 1):
            print('at time step {0} out of {1}'.format(tstep, steps_to_keep))
            pcoords = traj_array[tstep - 1].copy()
            evolver(
                    self.authToken,
                    ctypes.c_char_p(data_set.encode('ascii')),
                    ctypes.c_float(time_array[tstep-1]),
                    ctypes.c_float(time_array[tstep  ]),
                    ctypes.c_float(dt),
                    ctypes.c_int(sinterp), ctypes.c_int(npoints),
                    pcoords.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))),
                    result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))))
            print('got next position for time step {0}'.format(tstep))
            traj_array[tstep] = result_array
        return traj_array, time_array
    def getFilteredPosition(self,
            starttime = 0.0,
            endtime = 0.1,
            dt = 0.0004,
            filterwidth = (2*np.pi / 1024) * 5,
            point_coords = None,
            sinterp = 4,
            steps_to_keep = 1,
            data_set = 'isotropic1024coarse'):
        if not self.connection_on:
            print('you didn\'t connect to the database')
            sys.exit()
        if not (point_coords.shape[-1] == 3):
            print ('wrong number of values for coordinates in getPosition')
            sys.exit()
            return None
        npoints = point_coords.shape[0]
        result_array = np.empty((npoints, 3), dtype=np.float32)
        traj_array = np.empty((steps_to_keep+1, npoints, 3), dtype = np.float32)
        traj_array[0] = point_coords
        time_array = np.empty((steps_to_keep+1), dtype = np.float32)
        integration_time = (endtime - starttime) / steps_to_keep
        time_array[0] = starttime
        evolver = self.lib.getFilteredPosition
        for tstep in range(1, steps_to_keep + 1):
            print('at time step {0} out of {1}'.format(tstep, steps_to_keep))
            pcoords = traj_array[tstep - 1].copy()
            evolver(
                    self.authToken,
                    ctypes.c_char_p(data_set.encode('ascii')),
                    ctypes.c_float(starttime + (tstep - 1)*integration_time),
                    ctypes.c_float(starttime +  tstep     *integration_time),
                    ctypes.c_float(dt),
                    ctypes.c_float(filterwidth),
                    ctypes.c_int(npoints),
                    pcoords.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))),
                    result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))))
            traj_array[tstep] = result_array
            time_array[tstep] = starttime + tstep * integration_time
        return traj_array, time_array
    def getBlines(self,
            time = 0.0,
            ds = 0.0004,
            S = 0.1,
            points = np.array([[0, 0, 0]]).astype(np.float32),
            sinterp = 4,
            tinterp = 0,
            data_set = 'mhd1024'):
        """ use Heun solver to integrate magnetic field lines.
            better solvers can be written in the future.
        """
        if not self.connection_on:
            print('you didn\'t connect to the database')
            sys.exit()
        if not (points.shape[1] == 3 and len(points.shape) == 2):
            print ('wrong shape of initial condition in getBlines')
            sys.exit()
            return None
        iterations = int(S / ds)
        npoints = points.shape[0]
        history = np.zeros(
                (iterations+1, npoints, 3),
                dtype = np.float32)
        history[0] = points
        result_array = points.copy()
        get_data = getattr(self.lib, 'getMagneticField')
        def getBunit(point_coords):
            get_data(
                    self.authToken,
                    ctypes.c_char_p(data_set.encode('ascii')),
                    ctypes.c_float(time),
                    ctypes.c_int(sinterp), ctypes.c_int(tinterp), ctypes.c_int(npoints),
                    point_coords.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))),
                    result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))))
            bsize = np.maximum(np.sqrt(np.sum(result_array**2, axis = 1)), 1e-18)
            return result_array / bsize[:, None]
        for s in range(1, iterations + 1):
            bhat0 = getBunit(history[s-1])
            y = history[s-1] + ds * bhat0
            history[s] = history[s-1] + .5*ds * (bhat0 + getBunit(y))
        return history
    def getBline(self,
            time = 0.0,
            ds = 0.0004,
            S = 0.1,
            point = np.array([0, 0, 0]).astype(np.float32),
            sinterp = 4,
            tinterp = 0,
            data_set = 'mhd1024',
            out_of_domain = None):
        """ use Heun solver to integrate magnetic 1 magnetic field line, within a given domain.
            better solvers can be written in the future.
        """
        if not self.connection_on:
            print('you didn\'t connect to the database')
            sys.exit()
        if not (point.shape[0] == 3 and len(point.shape) == 1):
            print ('wrong shape of initial condition in getBline')
            sys.exit()
            return None
        iterations = int(S / abs(ds))
        history = np.zeros(
                (iterations+1, 3),
                dtype = np.float32)
        history[0] = point
        result_array = point.copy()
        get_data = getattr(self.lib, 'getMagneticField')
        def getBunit(point_coords):
            #print point_coords
            get_data(
                    self.authToken,
                    ctypes.c_char_p(data_set.encode('ascii')),
                    ctypes.c_float(time),
                    ctypes.c_int(sinterp), ctypes.c_int(tinterp), ctypes.c_int(1),
                    point_coords.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))),
                    result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))))
            #print 'get_data finished'
            bsize = max(np.sqrt(np.sum(result_array**2)), 1e-18)
            return result_array / bsize
        for s in range(1, iterations + 1):
            if not out_of_domain == None:
                if out_of_domain(history[s-1]):
                    return history[:s]
            y = history[s-1].copy()
            bhat0 = getBunit(y)
            y = history[s-1] + ds * bhat0
            history[s] = history[s-1] + .5*ds * (bhat0 + getBunit(y))
        return history
    def get2wayBline(self,
            time = 0.0,
            ds = 0.0004,
            S = 0.1,
            point = np.array([0, 0, 0]).astype(np.float32),
            sinterp = 4,
            tinterp = 0,
            data_set = 'mhd1024',
            out_of_domain = None):
        l0 = self.getBline(
                time = time,
                ds = ds,
                S = S,
                point = point,
                sinterp = sinterp,
                tinterp = tinterp,
                data_set = data_set,
                out_of_domain = out_of_domain)
        l1 = self.getBline(
                time = time,
                ds = -ds,
                S = S,
                point = point,
                sinterp = sinterp,
                tinterp = tinterp,
                data_set = data_set,
                out_of_domain = out_of_domain)
        return np.concatenate((l1[::-1], l0[1:]), axis = 0)
    def getBlineAlt(self,
            time, nsteps, ds,
            x0,
            sinterp = 4,
            tinterp = 0,
            data_set = 'mhd1024'):
        if not self.connection_on:
            print('you didn\'t connect to the database')
            sys.exit()
        if not (x0.shape[1] == 3 and len(x0.shape) == 2):
            print(('wrong shape of initial condition in getBlineAlt, ', x0.shape))
            sys.exit()
            return None
        npoints = x0.shape[0]
        result_array = np.empty((nsteps+1, npoints, 3), dtype=np.float32)
        result_array[0] = x0
        self.lib.getBline(
                self.authToken,
                ctypes.c_char_p(data_set.encode('ascii')),
                ctypes.c_float(time),
                ctypes.c_int(nsteps),
                ctypes.c_float(ds),
                ctypes.c_int(sinterp), ctypes.c_int(tinterp), ctypes.c_int(npoints),
                result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))))
        return result_array
    def getBlineSphereBounded(self,
            time = 0,
            S = 0.05,
            ds = 0.005,
            x0 = np.array([[0, 0, 0]], dtype = np.float32),
            sinterp = 4,
            tinterp = 0,
            data_set = 'mhd1024',
            origin = [0, 0, 0],
            radius = 1,
            both_ways = False):
        if not self.connection_on:
            print('you didn\'t connect to the database')
            sys.exit()
        if not (x0.shape[1] == 3 and len(x0.shape) == 2):
            print(('wrong shape of initial condition in getBlineSphereBounded, ', x0.shape))
            sys.exit()
            return None
        nsteps = int(S / abs(ds))
        npoints = x0.shape[0]
        result_array = np.empty((npoints, nsteps+1, 3), dtype=np.float32)
        length_array = np.empty((npoints), dtype = np.int32)
        result_array[:, 0] = x0
        self.lib.getSphericalBoundedBline(
                self.authToken,
                ctypes.c_char_p(data_set.encode('ascii')),
                ctypes.c_float(time),
                ctypes.c_int(nsteps),
                ctypes.c_float(ds),
                ctypes.c_int(sinterp), ctypes.c_int(tinterp), ctypes.c_int(npoints),
                result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))),
                length_array.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
                ctypes.c_float(origin[0]),
                ctypes.c_float(origin[1]),
                ctypes.c_float(origin[2]),
                ctypes.c_float(radius))
        bline_list = [result_array[p, :length_array[p]].copy()
                      for p in range(x0.shape[0])]
        if both_ways:
            result_array[:, 0] = x0
            self.lib.getSphericalBoundedBline(
                    self.authToken,
                    ctypes.c_char_p(data_set.encode('ascii')),
                    ctypes.c_float(time),
                    ctypes.c_int(nsteps),
                    ctypes.c_float(-ds),
                    ctypes.c_int(sinterp), ctypes.c_int(tinterp), ctypes.c_int(npoints),
                    result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))),
                    length_array.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
                    ctypes.c_float(origin[0]),
                    ctypes.c_float(origin[1]),
                    ctypes.c_float(origin[2]),
                    ctypes.c_float(radius))
            for l in range(len(bline_list)):
                bline_list[l] = np.concatenate((result_array[l, 1:length_array[l]][::-1].copy(), bline_list[l]), axis = 0)
        return bline_list
    def getBlineSphereBoundedDebug(self,
            time = 0,
            S = 0.05,
            ds = 0.005,
            x0 = np.array([[0, 0, 0]], dtype = np.float32),
            sinterp = 4,
            tinterp = 0,
            data_set = 'mhd1024',
            origin = [0, 0, 0],
            radius = 1,
            both_ways = False):
        if not self.connection_on:
            print('you didn\'t connect to the database')
            sys.exit()
        if not (x0.shape[1] == 3 and len(x0.shape) == 2):
            print(('wrong shape of initial condition in getBlineSphereBoundedDebug, ', x0.shape))
            sys.exit()
            return None
        nsteps = int(S / abs(ds))
        npoints = x0.shape[0]
        result_array = np.empty((npoints, nsteps+1, 3), dtype=np.float32)
        length_array = np.empty((npoints), dtype = np.int32)
        result_array[:, 0] = x0
        self.lib.getSphericalBoundedBlineDebug(
                self.authToken,
                ctypes.c_char_p(data_set.encode('ascii')),
                ctypes.c_float(time),
                ctypes.c_int(nsteps),
                ctypes.c_float(ds),
                ctypes.c_int(sinterp), ctypes.c_int(tinterp), ctypes.c_int(npoints),
                result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))),
                length_array.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
                ctypes.c_float(origin[0]),
                ctypes.c_float(origin[1]),
                ctypes.c_float(origin[2]),
                ctypes.c_float(radius))
        bline_list = [result_array[p, :length_array[p]].copy()
                      for p in range(x0.shape[0])]
        if both_ways:
            result_array[:, 0] = x0
            self.lib.getSphericalBoundedBlineDebug(
                    self.authToken,
                    ctypes.c_char_p(data_set.encode('ascii')),
                    ctypes.c_float(time),
                    ctypes.c_int(nsteps),
                    ctypes.c_float(-ds),
                    ctypes.c_int(sinterp), ctypes.c_int(tinterp), ctypes.c_int(npoints),
                    result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))),
                    length_array.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
                    ctypes.c_float(origin[0]),
                    ctypes.c_float(origin[1]),
                    ctypes.c_float(origin[2]),
                    ctypes.c_float(radius))
            for l in range(len(bline_list)):
                bline_list[l] = np.concatenate((result_array[l, 1:length_array[l]][::-1].copy(), bline_list[l]), axis = 0)
        return bline_list
    def getBlineRectBounded(self,
            time = 0,
            S = 0.05,
            ds = 0.005,
            x0 = np.array([[0, 0, 0]], dtype = np.float32),
            sinterp = 4,
            tinterp = 0,
            data_set = 'mhd1024',
            xmin = 0, xmax = 1,
            ymin = 0, ymax = 1,
            zmin = 0, zmax = 1,
            both_ways = False):
        if not self.connection_on:
            print('you didn\'t connect to the database')
            sys.exit()
        if not (x0.shape[1] == 3 and len(x0.shape) == 2):
            print(('wrong shape of initial condition in getBlineRectBounded, ', x0.shape))
            sys.exit()
            return None
        nsteps = int(S / abs(ds))
        npoints = x0.shape[0]
        result_array = np.empty((npoints, nsteps+1, 3), dtype=np.float32)
        length_array = np.empty((npoints), dtype = np.int32)
        result_array[:, 0] = x0
        self.lib.getRectangularBoundedBline(
                self.authToken,
                ctypes.c_char_p(data_set.encode('ascii')),
                ctypes.c_float(time),
                ctypes.c_int(nsteps),
                ctypes.c_float(ds),
                ctypes.c_int(sinterp), ctypes.c_int(tinterp), ctypes.c_int(npoints),
                result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))),
                length_array.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
                ctypes.c_float(xmin), ctypes.c_float(xmax),
                ctypes.c_float(ymin), ctypes.c_float(ymax),
                ctypes.c_float(zmin), ctypes.c_float(zmax))
        bline_list = [result_array[p, :length_array[p]].copy()
                      for p in range(x0.shape[0])]
        if both_ways:
            result_array[:, 0] = x0
            self.lib.getRectangularBoundedBline(
                    self.authToken,
                    ctypes.c_char_p(data_set.encode('ascii')),
                    ctypes.c_float(time),
                    ctypes.c_int(nsteps),
                    ctypes.c_float(-ds),
                    ctypes.c_int(sinterp), ctypes.c_int(tinterp), ctypes.c_int(npoints),
                    result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))),
                    length_array.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
                    ctypes.c_float(xmin), ctypes.c_float(xmax),
                    ctypes.c_float(ymin), ctypes.c_float(ymax),
                    ctypes.c_float(zmin), ctypes.c_float(zmax))
            for l in range(len(bline_list)):
                bline_list[l] = np.concatenate((result_array[l, 1:length_array[l]][::-1].copy(), bline_list[l]), axis = 0)
        return bline_list

max_stars_repo_name: arthursoprano/pyJHTDB

max_issues_repo_name: arthursoprano/pyJHTDB

max_forks_repo_name: arthursoprano/pyJHTDB


text:
import argparse
from typing import Text
import os
import struct
import glob

import onnx
from onnx import numpy_helper

__author__ = "Alexander Jung (University of Tuebingen, Chair for Embedded Systems)"

def main():

    parser = argparse.ArgumentParser(
        description="Tool to parse the reference input and output files from the onnx model zoo.")
    parser.add_argument(
        "--input",
        type=Text,
        required=True,
        help="Path to reference files."
    )
    args = parser.parse_args()

    reference_data_dir = args.input

    print("Using reference data from:", reference_data_dir)

    # Load inputs
    inputs_num = len(glob.glob(os.path.join(reference_data_dir, 'input_*.pb')))

    if inputs_num == 0:
        print("ERROR: There seems to be no *.pb file in the specified directory.")
        exit(1)

    for i in range(inputs_num):
        input_file = os.path.join(reference_data_dir, 'input_{}.pb'.format(i))
        tensor = onnx.TensorProto()
        with open(input_file, 'rb') as f:
            tensor.ParseFromString(f.read())
            tensor_array = numpy_helper.to_array(tensor)

        if len(tensor_array.shape) != 4 or tensor_array.shape[0] != 1:
            print("ERROR: Unsupported input shape: {}".format(tensor_array.shape))
            return
        else:
            tensor_shape = tensor_array.shape

            packed_file = list(bytes())
            magic_number = bytes("FCI\n", "ascii")
            packed_file.append(struct.pack('{}s'.format(len(magic_number)), magic_number))
            # TODO: Maybe insert file name to identify the model for which the data is intended
            packed_file.append(struct.pack('i', tensor_shape[0]))
            packed_file.append(struct.pack('i', tensor_shape[1]))
            packed_file.append(struct.pack('i', tensor_shape[2]))
            packed_file.append(struct.pack('i', tensor_shape[3]))

            for channel in tensor_array[0]:
                for row in channel:
                    packed_file.append(struct.pack('f'*len(row), *row))

            tupac = bytes("end\n", "ascii")
            packed_file.append(struct.pack('{}s'.format(len(tupac)), tupac))

            with open(os.path.join(args.input, "input_{}.data".format(i)), "wb") as f:
                for packed_struct in packed_file:
                    f.write(packed_struct)

            print("Input data shape: {}".format(tensor_shape))

    # Load reference outputs
    ref_outputs_num = len(glob.glob(os.path.join(reference_data_dir, 'output_*.pb')))

    if ref_outputs_num == 0:
        print("ERROR: There seems to be no *.pb file in the specified directory.")
        exit(1)

    for i in range(ref_outputs_num):
        output_file = os.path.join(reference_data_dir, 'output_{}.pb'.format(i))
        tensor = onnx.TensorProto()
        with open(output_file, 'rb') as f:
            tensor.ParseFromString(f.read())
            tensor_array = numpy_helper.to_array(tensor)

        if len(tensor_array.shape) != 2 or tensor_array.shape[0] != 1:
            print("ERROR: Unsupported output shape: {}".format(tensor_array.shape))
            return
        else:
            tensor_shape = tensor_array.shape

            packed_file = list(bytes())
            magic_number = bytes("FCO\n", "ascii")
            packed_file.append(struct.pack('{}s'.format(len(magic_number)), magic_number))
            # TODO: Maybe insert file name to identify the model for which the data is intended
            packed_file.append(struct.pack('i', tensor_shape[0]))
            packed_file.append(struct.pack('i', tensor_shape[1]))

            packed_file.append(struct.pack('f'*len(tensor_array[0]), *tensor_array[0]))

            tupac = bytes("end\n", "ascii")
            packed_file.append(struct.pack('{}s'.format(len(tupac)), tupac))

            with open(os.path.join(args.input, "output_{}.data".format(i)), "wb") as f:
                for packed_struct in packed_file:
                    f.write(packed_struct)

            print("Output data shape: {}".format(tensor_shape))

if __name__ == '__main__':
    main()

ext: py

max_stars_repo_name: ekut-es/pico-cnn

max_issues_repo_name: ekut-es/pico-cnn

max_forks_repo_name: ekut-es/pico-cnn

reason: import numpy


text:
# This file is part of the mt5se package
#  mt5se home: https://github.com/paulo-al-castro/mt5se
# Author: Paulo Al Castro
# Date: 2020-11-17

import pandas as pd
import numpy as np
import mt5se.mt5se as se
from scipy import stats

def rsi(returns):
    if type(returns)==pd.core.frame.DataFrame:
        returns=se.get_returns(returns)
    u=0.0
    uc=0
    d=0.0
    dc=0
    for r in returns:
        if r>=0:
            u=u+r
            uc=uc+1
        else:
            d=d+abs(r)
            dc=dc+1
    if uc>0:
        u=u/uc
    if dc>0:
        d=d/dc
    if d==0:
        d=1.0
    ifr=100*( 1 - 1/(1+u/d))
    return ifr

#	Returns the angular coefficient of linear regression (slope)
#      for a serie of prices in regular intervals
def slope(serie):
    x=np.array(range(len(serie)))
    #y=np.array(serie)
    s=stats.linregress(x,serie)
    return s.slope

# equals to slope
def trend(serie):
    return slope(serie)

reason: import numpy,from scipy


text:
import sys
import os.path
import argparse
import datetime
from lib.trail import Trail
from lib.party import Party
from lib.world import World

from numpy import mean

def main(args):
    trials = 1000
    strategy = 'greatest need'
    start_datetime = datetime.datetime.strptime(args.start_date_time, '%Y-%m-%d %H:%M')

    data_path = os.path.abspath('data/')
    trail_file_name = os.path.join(data_path, 'belly_river_trail.json')
    terrain_file_name = os.path.join(data_path, 'belly_river_terrain.json')
    party_file_name = os.path.join(data_path, 'belly_river_party.json')

    stats = {'survivors': list(),
             'travel time': list(),
             'bear attacks': list()}

    for t in range(trials):

        trail = Trail(start_datetime, trail_file_name, terrain_file_name)
        party = Party(start_datetime, party_file_name, trail)
        world = World(start_datetime, party, trail)

        # simulate until destination or death
        continue_trail = True

        while continue_trail:
            world.update(strategy)

            if world.party.condition in ('arrived', 'dead'):
                # update Powell metrics
                continue_trail = False

            #print world.date_and_time, world.party.current_stop.mile_marker, world.party.current_stop.name
            #print world.party.condition, world.party.number_alive()
            #world.party.print_party_status()
            #print
        stats['travel time'].append((world.date_and_time - start_datetime).total_seconds())
        stats['survivors'].append(world.party.number_alive())
        stats['bear attacks'].append(world.event_counter['bear attack'].get(True, 0))

    #print 'Total trail time:', world.date_and_time - start_datetime

    print 'number of trips:', len(stats['survivors'])
    print 'number of bear attacks', sum(stats['bear attacks'])
    print 'mean survivors:', mean(stats['survivors'])
    print 'mean travel time (h):', mean(stats['travel time'])/60.0/60.0

    # present Powell metrics
    bears_per_mile = sum(stats['bear attacks']) / 25.67 / trials
    mortality = (4.0*trials - sum(stats['survivors'])) / trials
    powells = 1.25 * (25.67 + 4495 / 1000) * (mortality * 2) * 1000.0
    HARM = powells / 25.67

    print bears_per_mile, powells, HARM

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument("--start-date-time", help="party start date YYYY-MM-DD hh:mm")
    args = parser.parse_args()
    main(args)

max_stars_repo_name: rwsharp/TORT

max_issues_repo_name: rwsharp/TORT

max_forks_repo_name: rwsharp/TORT

reason: from numpy


text:
# -*- coding: utf-8 -*-
#/usr/bin/python2
'''
By kyubyong park. kbpark.linguist@gmail.com.
https://www.github.com/kyubyong/tacotron_asr
'''
from __future__ import print_function

import codecs
import copy
import re

import librosa

from hyperparams import Hyperparams as hp
import numpy as np
import tensorflow as tf

def get_spectrogram(sound_fpath):
    '''Extracts melspectrogram and magnitude from given `sound_file`.
    Args:
      sound_fpath: A string. Full path of a sound file.

    Returns:
      Transposed S: A 2d array. A transposed melspectrogram with shape of (T, n_mels)
      Transposed magnitude: A 2d array. A transposed magnitude spectrogram
        with shape of (T, 1+hp.n_fft//2)
    '''
    # Loading sound file
    y, sr = librosa.load(sound_fpath, sr=None) # or set sr to hp.sr.

    # stft. D: (1+n_fft//2, T)
    D = librosa.stft(y=y,
                     n_fft=hp.n_fft,
                     hop_length=hp.hop_length,
                     win_length=hp.win_length)

    # magnitude spectrogram
    magnitude = np.abs(D) #(1+n_fft/2, T)

    # power spectrogram
    power = magnitude**2

    # mel spectrogram
    S = librosa.feature.melspectrogram(S=power, n_mels=hp.n_mels) #(n_mels, T)

    return np.transpose(S.astype(np.float32))

def shift_by_one(inputs):
    '''Shifts the content of `inputs` to the right by one
      so that it becomes the decoder inputs.

    Args:
      inputs: A 3d tensor with shape of [N, T, C]

    Returns:
      A 3d tensor with the same shape and dtype as `inputs`.
    '''
    return tf.concat((tf.zeros_like(inputs[:, :1]), inputs[:, :-1]), 1)

def reduce_frames(arry, r):
    '''Reduces and adjust the shape and content of `arry` according to r.

    Args:
      arry: A 2d array with shape of [T, C]
      r: Reduction factor

    Returns:
      A 2d array with shape of [-1, C*r]
    '''
    T, C = arry.shape
    num_paddings = hp.r - (T % r) if T % r != 0 else 0

    padded = np.pad(arry, [[0, num_paddings], [0, 0]], 'constant')
    output = np.reshape(padded, (-1, C*r))
    return output

def restore_shape(arry, r):
    '''Restore and adjust the shape and content of `inputs` according to r.
    Args:
      arry: A 3d array with shape of [N, T, C]
      r: Reduction factor

    Returns:
      A 3d tensor with shape of [-1, C*r]
    '''
    N, T, C = arry.shape
    return arry.reshape((N, -1, C//r))


text:
import itertools
import pathlib

import numpy

from .learning import smape
from . import gaiting
from ..utils import Namespace, first_sign_change, make_progressbar

def terminated(rec, lookback, smooth, threshold, rolling='mean',
               average_window=None, gait_func=gaiting.smr1e):
    if average_window is None:
        average_window = lookback
    result = Namespace(**locals())

    gp = rec.flatten_gen_params()
    lookback_step = int(lookback / rec.rc.gen_learning_rate)
    smooth_step = int(smooth / rec.rc.gen_learning_rate)
    window_step = int(average_window / rec.rc.gen_learning_rate)

    diffs = gait_func(gp, lookback_step, smooth_step, rolling)
    diffs /= lookback
    i_cross = first_sign_change(diffs - threshold)
    if i_cross is None:
        result.first_cross = None
        result.param = numpy.empty_like(gp[0])
        result.param[:] = numpy.nan
        result.cost_norm2 = numpy.nan
        result.cost_smape = numpy.nan
    else:
        result.first_cross = first_cross = i_cross + lookback_step
        result.param = gp[first_cross - window_step: first_cross].mean(axis=0)
        p_gen = result.param
        p_true = rec.flatten_true_params()
        result.cost_norm2 = numpy.linalg.norm(p_gen - p_true)
        result.cost_smape = smape(p_gen, p_true)
    return result

def sweep_terminators(records, gait_funcs=[gaiting.smr1e, gaiting.pw1e],
                      lookbacks=[2, 1, 0.5], smooths=[2, 1, 0.5],
                      thresholds=[0.05, 0.01, 0.005, 0.001],
                      rollings=['mean']):
    names = ['gait_func', 'lookback', 'smooth', 'threshold', 'rolling']
    axes = [gait_funcs, lookbacks, smooths, thresholds, rollings]
    for rec in records:
        rec_id = str(rec.datastore.directory)
        for values in itertools.product(*axes):
            row = kwargs = dict(zip(names, values))
            result = terminated(rec, **kwargs)
            row['rec_id'] = rec_id
            row.update(vars(result))
            row['gait'] = row.pop('gait_func').__name__
            del row['rec']  # save memory, in case `records` is lazy
            yield row

def simulate_terminations(records=None, glob=None, progress=False, **kwargs):
    import pandas
    if records is None:
        assert glob is not None
        from ..loaders import load_records
        paths = make_progressbar(not progress)(list(pathlib.Path().glob(glob)))
        records = map(load_records, paths)
    else:
        records = make_progressbar(not progress)(records)
    return pandas.DataFrame(sweep_terminators(records, **kwargs))

def simulate_terminations_snr(*args,
                              gait_funcs=[gaiting.pw1e], rollings=['snr'],
                              thresholds=[1.0, 0.75, 0.5, 0.25],
                              **kwargs):
    """
    Same as `simulate_terminations` but with rolling SNR.

    Options are default to the values that work well with rolling SNR.
    """
    return simulate_terminations(*args,
                                 gait_funcs=gait_funcs,
                                 rollings=rollings,
                                 thresholds=thresholds,
                                 **kwargs)

max_stars_repo_name: ahmadianlab/tc-gan

max_issues_repo_name: ahmadianlab/tc-gan

max_forks_repo_name: ahmadianlab/tc-gan


text:
import numpy as np
import cv2
import sys
import math

def distance(x, y, i, j):
    return np.sqrt((x-i)**2 + (y-j)**2)

def gaussian(x, sigma):
    # print(sum(x**2))
    if type(x) == np.ndarray:
        sq = sum(x ** 2)
    else:
        sq = x**2

    return (1.0 / (2 * math.pi * (sigma ** 2))) * math.exp(- sq / (2 * sigma ** 2))

def apply_bilateral_filter(source, filtered_image, x, y, diameter, sigma_i, sigma_s):
    hl = int(diameter/2)
    i_filtered = 0
    Wp = 0
    i = 0
    while i < diameter:
        j = 0
        while j < diameter:
            neighbour_x = x - (hl - i)
            neighbour_y = y - (hl - j)
            if neighbour_x >= len(source):
                neighbour_x -= len(source)
            if neighbour_y >= len(source[0]):
                neighbour_y -= len(source[0])
            gi = gaussian(source[neighbour_x][neighbour_y] - source[x][y], sigma_i)
            gs = gaussian(distance(neighbour_x, neighbour_y, x, y), sigma_s)
            w = gi * gs
            i_filtered += source[neighbour_x][neighbour_y] * w
            Wp += w
            j += 1
        i += 1
    i_filtered = i_filtered / Wp
    filtered_image[x][y] = i_filtered

def bilateral_filter_own(source, filter_diameter, sigma_i, sigma_s):
    filtered_image = np.zeros(source.shape)

    i = 0
    while i < len(source):
        print(i)
        j = 0
        while j < len(source[0]):
            apply_bilateral_filter(source, filtered_image, i, j, filter_diameter, sigma_i, sigma_s)
            j += 1
        i += 1
    return filtered_image

if __name__ == "__main__":
    src = cv2.imread(str(sys.argv[1]), 0)
    filtered_image_OpenCV = cv2.bilateralFilter(src, 5, 12.0, 16.0)
    cv2.imwrite("original_image_grayscale.png", src)
    cv2.imwrite("filtered_image_OpenCV.png", filtered_image_OpenCV)
    filtered_image_own = bilateral_filter_own(src, 5, 12.0, 16.0)
    cv2.imwrite("filtered_image_own.png", filtered_image_own)


text:
from __future__ import print_function
import os, sys
import argparse
import json
import numpy as np

def convert_to_dict(data, is_test=False):
    ret_dict = {}
    if not is_test:
        ret_dict["Y"] = data[:, 0].tolist()
        ret_dict["X"] = data[:, 1:].tolist()
    else:
        ret_dict["X"] = data.tolist()
    return ret_dict

def extract_features(args):
    train_file = os.path.join(args.data_dir, args.training_file)
    test_file = os.path.join(args.data_dir, args.test_file)
    output_dir = args.output_dir
    train_output_file = os.path.join(output_dir, "train.json")
    test_output_file = os.path.join(output_dir, "test.json")

    if not os.path.exists(train_file):
        print("Training file does not exist")
        sys.exit(-1)
    if not os.path.exists(test_file):
        print("Test file does not exist")
        sys.exit(-1)

    training_data = np.genfromtxt(train_file, delimiter=',', skip_header=1)
    test_data = np.genfromtxt(test_file, delimiter=',', skip_header=1)

    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    with open(train_output_file, 'w') as fh:
        json.dump(convert_to_dict(training_data), fh)

    with open(test_output_file, 'w') as fh:
        json.dump(convert_to_dict(test_data), fh)

    if os.path.exists(train_output_file) and os.path.exists(test_output_file):
        sys.exit(0)
    else:
        print("Extract Features Fail")
        sys.exit(-1)

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument("--training_file", default="train.csv")
    parser.add_argument("--test_file", default="test.csv")
    parser.add_argument("--output_dir", default="output/features")
    parser.add_argument('data_dir')

    args = parser.parse_args()
    extract_features(args)

max_stars_repo_name: cjacoby/mnist-eval

max_issues_repo_name: cjacoby/mnist-eval

max_forks_repo_name: cjacoby/mnist-eval


text:
# -*- coding: utf-8 -*-

import numpy as np
import astropy.units as u
import pkg_resources

from astropy.io import ascii
from astropy.modeling.tabular import tabular_model

from .baseclasses import BaseAtttauVModel
from .helpers import _test_valid_x_range

__all__ = ["WG00"]

x_range_WG00 = [0.1, 3.0001]

class WG00(BaseAtttauVModel):
    r"""
    Attenuation curve of Witt & Gordon (2000)

    Parameters
    ----------
    tau_v: float
        optical depth in V band

    Raises
    ------
    InputParameterError
       Input Av values outside of defined range

    Notes
    -----
    From Witt & Gordon (2000, ApJ, Volume 528, pp. 799-816)

    Example:

    .. plot::
        :include-source:

        import numpy as np
        import matplotlib.pyplot as plt
        import astropy.units as u

        from dust_attenuation.radiative_transfer import WG00

        fig, ax = plt.subplots(1,2, figsize=(10,6))

        # generate the curves and plot them
        # Use 1/microns for a better sampling
        x = np.arange(0.35,10.0,0.1)/u.micron

        x_Vband = 0.55 # microns

        tau_Vs = [0.25,0.4,1.1,17.0,46.0]
        for tau_V in tau_Vs[::-1]:
           att_model = WG00(tau_V = tau_V, geometry = 'cloudy',
                            dust_type = 'mw', dust_distribution = 'clumpy')
           ax[0].plot(x,att_model(1/x),label=r'$\tau_V$ = %.2f mag' % (tau_V))
           ax[1].plot(x,att_model(1/x)/att_model(x_Vband),
                      label=r'$\tau_V$ = %.2f mag' % (tau_V))

        ax[0].set_xlabel(r'$x$ [$\mu m^{-1}$]')
        ax[0].set_ylabel(r'$Att(x)$ [mag]')
        ax[1].set_xlabel(r'$x$ [$\mu m^{-1}$]')
        ax[1].set_ylabel(r'$Att(x)/Att_V$')

        ax[0].legend(loc='best')
        ax[1].legend(loc='best')
        fig.suptitle(r'CLOUDY / MW / clumpy model',size=15)
        plt.tight_layout()
        fig.subplots_adjust(top=0.88)

        plt.show()

    """

    tau_V_range = [0.25, 50.0]
    x_range = x_range_WG00

    def __init__(
        self, tau_V, geometry="dusty", dust_type="mw", dust_distribution="clumpy"
    ):
        """
        Load the attenuation curves for a given geometry, dust type and
        dust distribution.

        Parameters
        ----------
        tau_V: float
           optical depth in V band

        geometry: string
           'shell', 'cloudy' or 'dusty'

        dust_type: string
           'mw' or 'smc'

        dust_distribution: string
           'homogeneous' or 'clumpy'

        Returns
        -------
        Attx: np array (float)
            Att(x) attenuation curve [mag]

        """
        # Ensure strings are lower cases
        self.geometry = geometry.lower()
        self.dust_type = dust_type.lower()
        self.dust_distribution = dust_distribution.lower()

        data_path = pkg_resources.resource_filename("dust_attenuation", "data/WG00/")

        data = ascii.read(data_path + self.geometry + ".txt", header_start=0)

        if self.dust_type == "mw":
            start = 0
        elif self.dust_type == "smc":
            start = 25

        # Column names
        tau_colname = "tau"
        tau_att_colname = "tau_att"
        fsca_colname = "f(sca)"
        fdir_colname = "f(dir)"
        fesc_colname = "f(esc)"

        if self.dust_distribution == "clumpy":
            tau_att_colname += "_c"
            fsca_colname += "_c"
            fdir_colname += "_c"
            fesc_colname += "_c"

        elif self.dust_distribution == "homogeneous":
            tau_att_colname += "_h"
            fsca_colname += "_h"
            fdir_colname += "_h"
            fesc_colname += "_h"

        tau_att_list = []
        tau_list = []
        fsca_list = []
        fdir_list = []
        fesc_list = []

        len_data = len(data["lambda"])
        # number of lines between 2 models
        steps = 25

        counter = start
        while counter < len_data:
            tau_att_list.append(
                np.array(data[tau_att_colname][counter : counter + steps])
            )
            tau_list.append(np.array(data[tau_colname][counter : counter + steps]))
            fsca_list.append(np.array(data[fsca_colname][counter : counter + steps]))
            fdir_list.append(np.array(data[fdir_colname][counter : counter + steps]))
            fesc_list.append(np.array(data[fesc_colname][counter : counter + steps]))

            counter += int(2 * steps)

        # Convert to np.array and take transpose to have (wvl, tau_V)
        tau_att_table = np.array(tau_att_list).T
        tau_table = np.array(tau_list).T
        fsca_table = np.array(fsca_list).T
        fdir_table = np.array(fdir_list).T
        fesc_table = np.array(fesc_list).T

        # wavelength grid. It is the same for all the models
        wvl = np.array(data["lambda"][0:25])
        self.wvl_grid = wvl

        # Grid for the optical depth
        tau_V_grid = np.array(
            [
                0.25,
                0.5,
                0.75,
                1.0,
                1.5,
                2.0,
                2.5,
                3.0,
                3.5,
                4.0,
                4.5,
                5.0,
                5.5,
                6.0,
                7.0,
                8.0,
                9.0,
                10.0,
                15.0,
                20.0,
                25.0,
                30.0,
                35.0,
                40.0,
                45.0,
                50.0,
            ]
        )

        # Create a 2D tabular model for tau_att and all flux fraction
        tab = tabular_model(2, name="2D_table")

        # Values corresponding to the x and y grid points
        gridpoints = (wvl, tau_V_grid)

        self.model = tab(
            gridpoints,
            lookup_table=tau_att_table,
            name="tau_att_WG00",
            bounds_error=False,
            fill_value=None,
            method="linear",
        )

        self.tau = tab(
            gridpoints,
            lookup_table=tau_table,
            name="tau_WG00",
            bounds_error=False,
            fill_value=None,
            method="linear",
        )

        self.fsca = tab(
            gridpoints,
            lookup_table=fsca_table,
            name="fsca_WG00",
            bounds_error=False,
            fill_value=None,
            method="linear",
        )

        self.fdir = tab(
            gridpoints,
            lookup_table=fdir_table,
            name="fdir_WG00",
            bounds_error=False,
            fill_value=None,
            method="linear",
        )

        self.fesc = tab(
            gridpoints,
            lookup_table=fesc_table,
            name="fesc_WG00",
            bounds_error=False,
            fill_value=None,
            method="linear",
        )

        # In Python 2: super(WG00, self)
        # In Python 3: super() but super(WG00, self) still works
        super(WG00, self).__init__(tau_V=tau_V)

    def evaluate(self, x, tau_V):
        """
        WG00 function

        Parameters
        ----------
        x: float
           expects either x in units of wavelengths or frequency
           or assumes wavelengths in [micron]

           internally microns are used

        tau_V: float
           optical depth in V band

        Returns
        -------
        Attx: np array (float)
            Att(x) attenuation curve [mag]

        Raises
        ------
        ValueError
           Input x values outside of defined range
        """
        # convert to wavenumbers (1/micron) if x input in units
        # otherwise, assume x in appropriate wavenumber units
        with u.add_enabled_equivalencies(u.spectral()):
            x_quant = u.Quantity(x, u.micron, dtype=np.float64)

        # strip the quantity to avoid needing to add units to all the
        #    polynomical coefficients
        x = x_quant.value

        # check that the wavenumbers are within the defined range
        _test_valid_x_range(x, self.x_range, "WG00")

        # setup the ax vectors
        n_x = len(x)

        xinterp = 1e4 * x
        yinterp = tau_V * np.ones(n_x)

        taux = self.model(xinterp, yinterp)

        # Convert optical depth to attenuation
        Attx = 1.086 * taux

        return Attx

    def get_extinction(self, x, tau_V):
        """
        Return the extinction at a given wavelength and
        V-band optical depth.

        Parameters
        ----------
        x: float
           expects either x in units of wavelengths or frequency
           or assumes wavelengths in [micron]

           internally microns are used

        tau_V: float
           optical depth in V band

        Returns
        -------
        ext: np array (float)
            ext(x) extinction curve [mag]

        Raises
        ------
        ValueError
           Input x values outside of defined range
        """
        # convert to wavenumbers (1/micron) if x input in units
        # otherwise, assume x in appropriate wavenumber units
        with u.add_enabled_equivalencies(u.spectral()):
            x_quant = u.Quantity(x, u.micron, dtype=np.float64)

        # strip the quantity to avoid needing to add units to all the
        #    polynomical coefficients
        x = x_quant.value

        # check that the wavenumbers are within the defined range
        _test_valid_x_range(x, self.x_range, "WG00")

        # setup the ax vectors
        x = np.atleast_1d(x)
        n_x = len(x)

        xinterp = 1e4 * x
        yinterp = tau_V * np.ones(n_x)

        return self.tau(xinterp, yinterp) * 1.086

    def get_fsca(self, x, tau_V):
        """
        Return the scattered flux fraction  at a given wavelength and
        V-band optical depth.

        Parameters
        ----------
        x: float
           expects either x in units of wavelengths or frequency
           or assumes wavelengths in [micron]

           internally microns are used

        tau_V: float
           optical depth in V band

        Returns
        -------
        fsca: np array (float)
            fsca(x) scattered flux fraction

        Raises
        ------
        ValueError
           Input x values outside of defined range
        """
        # convert to wavenumbers (1/micron) if x input in units
        # otherwise, assume x in appropriate wavenumber units
        with u.add_enabled_equivalencies(u.spectral()):
            x_quant = u.Quantity(x, u.micron, dtype=np.float64)

        # strip the quantity to avoid needing to add units to all the
        #    polynomical coefficients
        x = x_quant.value

        # check that the wavenumbers are within the defined range
        _test_valid_x_range(x, self.x_range, "WG00")

        # setup the ax vectors
        x = np.atleast_1d(x)
        n_x = len(x)

        xinterp = 1e4 * x
        yinterp = tau_V * np.ones(n_x)

        return self.fsca(xinterp, yinterp)

    def get_fdir(self, x, tau_V):
        """
        Return the direct attenuated stellar flux fraction  at a given
        wavelength and V-band optical depth.

        Parameters
        ----------
        x: float
           expects either x in units of wavelengths or frequency
           or assumes wavelengths in [micron]

           internally microns are used

        tau_V: float
           optical depth in V band

        Returns
        -------
        fsca: np array (float)
            fsca(x) scattered flux fraction

        Raises
        ------
        ValueError
           Input x values outside of defined range
        """
        # convert to wavenumbers (1/micron) if x input in units
        # otherwise, assume x in appropriate wavenumber units
        with u.add_enabled_equivalencies(u.spectral()):
            x_quant = u.Quantity(x, u.micron, dtype=np.float64)

        # strip the quantity to avoid needing to add units to all the
        #    polynomical coefficients
        x = x_quant.value

        # check that the wavenumbers are within the defined range
        _test_valid_x_range(x, self.x_range, "WG00")

        # setup the ax vectors
        x = np.atleast_1d(x)
        n_x = len(x)

        xinterp = 1e4 * x
        yinterp = tau_V * np.ones(n_x)

        return self.fdir(xinterp, yinterp)

    def get_fesc(self, x, tau_V):
        """
        Return the total escaping flux fraction  at a given wavelength and
        V-band optical depth.

        Parameters
        ----------
        x: float
           expects either x in units of wavelengths or frequency
           or assumes wavelengths in [micron]

           internally microns are used

        tau_V: float
           optical depth in V band

        Returns
        -------
        fsca: np array (float)
            fsca(x) scattered flux fraction

        Raises
        ------
        ValueError
           Input x values outside of defined range
        """
        # convert to wavenumbers (1/micron) if x input in units
        # otherwise, assume x in appropriate wavenumber units
        with u.add_enabled_equivalencies(u.spectral()):
            x_quant = u.Quantity(x, u.micron, dtype=np.float64)

        # strip the quantity to avoid needing to add units to all the
        #    polynomical coefficients
        x = x_quant.value

        # check that the wavenumbers are within the defined range
        _test_valid_x_range(x, self.x_range, "WG00")

        # setup the ax vectors
        x = np.atleast_1d(x)
        n_x = len(x)

        xinterp = 1e4 * x
        yinterp = tau_V * np.ones(n_x)

        return self.fesc(xinterp, yinterp)

    def get_albedo(self, x):
        """
        Return the albedo in function of wavelength for the corresponding
        dust type (SMC or MW). The albedo gives the probability a photon
        is scattered from a dust grain.

        Parameters
        ----------
        x: float
           expects either x in units of wavelengths or frequency
           or assumes wavelengths in [micron]

           internally microns are used

        Returns
        -------
        albedo: np array (float)
            alb(x) albedo

        Raises
        ------
        ValueError
           Input x values outside of defined range
        """
        # convert to wavenumbers (1/micron) if x input in units
        # otherwise, assume x in appropriate wavenumber units
        with u.add_enabled_equivalencies(u.spectral()):
            x_quant = u.Quantity(x, u.micron, dtype=np.float64)

        # strip the quantity to avoid needing to add units to all the
        #    polynomical coefficients
        x = x_quant.value

        # check that the wavenumbers are within the defined range
        _test_valid_x_range(x, self.x_range, "WG00")

        # setup the ax vectors
        x = np.atleast_1d(x)

        alb_MW = np.array(
            [
                0.320,
                0.409,
                0.481,
                0.526,
                0.542,
                0.536,
                0.503,
                0.432,
                0.371,
                0.389,
                0.437,
                0.470,
                0.486,
                0.499,
                0.506,
                0.498,
                0.502,
                0.491,
                0.481,
                0.500,
                0.473,
                0.457,
                0.448,
                0.424,
                0.400,
            ]
        )

        alb_SMC = np.array(
            [
                0.400,
                0.449,
                0.473,
                0.494,
                0.508,
                0.524,
                0.529,
                0.528,
                0.523,
                0.520,
                0.516,
                0.511,
                0.505,
                0.513,
                0.515,
                0.498,
                0.494,
                0.489,
                0.484,
                0.493,
                0.475,
                0.465,
                0.439,
                0.417,
                0.400,
            ]
        )

        if self.dust_type == "smc":
            albedo = alb_SMC
        elif self.dust_type == "mw":
            albedo = alb_MW

        tab = tabular_model(1, name="Tabular1D")
        alb_fit = tab(
            self.wvl_grid,
            lookup_table=albedo,
            name="albedo",
            bounds_error=False,
            fill_value=None,
            method="linear",
        )

        xinterp = 1e4 * x

        return alb_fit(xinterp)

    def get_scattering_phase_function(self, x):
        """
        Return the scattering phase function in function of wavelength for the
        corresponding dust type (SMC or MW). The scattering phase
        function gives the angle at which the photon scatters.

        Parameters
        ----------
        x: float
           expects either x in units of wavelengths or frequency
           or assumes wavelengths in [micron]

           internally microns are used

        Returns
        -------
        g: np array (float)
            g(x) scattering phase function

        Raises
        ------
        ValueError
           Input x values outside of defined range
        """
        # convert to wavenumbers (1/micron) if x input in units
        # otherwise, assume x in appropriate wavenumber units
        with u.add_enabled_equivalencies(u.spectral()):
            x_quant = u.Quantity(x, u.micron, dtype=np.float64)

        # strip the quantity to avoid needing to add units to all the
        #    polynomical coefficients
        x = x_quant.value

        # check that the wavenumbers are within the defined range
        _test_valid_x_range(x, self.x_range, "WG00")

        # setup the ax vectors
        x = np.atleast_1d(x)

        g_MW = np.array(
            [
                0.800,
                0.783,
                0.767,
                0.756,
                0.745,
                0.736,
                0.727,
                0.720,
                0.712,
                0.707,
                0.702,
                0.697,
                0.691,
                0.685,
                0.678,
                0.646,
                0.624,
                0.597,
                0.563,
                0.545,
                0.533,
                0.511,
                0.480,
                0.445,
                0.420,
            ]
        )

        g_SMC = np.array(
            [
                0.800,
                0.783,
                0.767,
                0.756,
                0.745,
                0.736,
                0.727,
                0.720,
                0.712,
                0.707,
                0.702,
                0.697,
                0.691,
                0.685,
                0.678,
                0.646,
                0.624,
                0.597,
                0.563,
                0.545,
                0.533,
                0.511,
                0.480,
                0.445,
                0.420,
            ]
        )

        if self.dust_type == "smc":
            g = g_SMC
        elif self.dust_type == "mw":
            g = g_MW

        tab = tabular_model(1, name="Tabular1D")
        g_fit = tab(
            self.wvl_grid,
            lookup_table=g,
            name="albedo",
            bounds_error=False,
            fill_value=None,
            method="linear",
        )

        xinterp = 1e4 * x

        return g_fit(xinterp)

reason: import numpy,import astropy,from astropy


text:
from importData import importMyData
import numpy as np
from abc import ABC

from mainFile import *

class KNN(Classifier):
	knn = 0
	def __init__(self, data, labels, knn):
		self.trainData = data
		self.trainLabels = labels
		self.knn = knn

	def validate(self, sample, label):
		neighbors = self.closestKNeighbors(sample)
		neighbors = np.sort(neighbors)
	    # find the max frequency using
	    # linear traversal
		max_count = 1
		res = neighbors[0]
		curr_count = 1
		n = neighbors.shape[0]
		for i in range(1, neighbors.shape[0]):
		    if (neighbors[i] == neighbors[i - 1]):
		        curr_count += 1

		    else:
		        if (curr_count > max_count):
		            max_count = curr_count
		            res = neighbors[i - 1]

		        curr_count = 1

		# If last element is most frequent
		if (curr_count > max_count):

		    max_count = curr_count
		    res = neighbors[n - 1]
		print(label, neighbors)
		if (res == label):
			return True
		else:
			return False

	def closestKNeighbors(self, sample):
		# print((self.trainData - sample).shape)
		errors_feature_vec = np.linalg.norm(self.trainData - sample, axis=1)
		# print(errors_feature_vec)
		idx = np.argsort(errors_feature_vec)
		# idx = idx[::-1]
		return self.trainLabels[idx][0:self.knn]

if __name__=="__main__":
	np.random.seed(0)
	dsplit = 0.8
	knn = 1
	threshold = 0.9
	task = 1 # 1 - person idenfitication, 2 - expression vs nuetral
	data_choice = 'dD' # dD, iD, pD
	preprocess = 'pca'
	data_train, labels_train, data_test, labels_test = importMyData(data_choice, task, preprocess, threshold, dsplit)
	# print(data_train.shape, labels_train.shape, data_test.shape, labels_test.shape)

	model = KNN(data_train, labels_train, knn)

	correct_preds = 0
	for i in range(data_test.shape[0]):
		# print(data_test[i].shape , labels_test[i].shape)
		# input()
		if(model.validate(data_test[i],labels_test[i])):
			correct_preds += 1

	print(" KNN accuracy: ", correct_preds/data_test.shape[0])


text:
"""
Some helper functions for unit testing
"""
import os
import subprocess
import time
import socket

from numpy import ndarray, float, float32, float64, complex, complex64, \
    complex128

RSERVE_PATH = subprocess.check_output(
    ['R', '--vanilla', '--slave', '-e',
     'cat(system.file(package="Rserve", "libs", '
     '.Platform$r_arch, "Rserve.dbg"))'])
HERE_PATH = os.path.dirname(os.path.realpath(__file__))

# Use different port from default to avoid clashes with regular Rserve
# running on same machine:
RPORT = 6355

def start_pyRserve():
    """Setup connection to remote Rserve for unittesting"""
    # Start Rserve
    rProc = subprocess.Popen(
        ['R', 'CMD', RSERVE_PATH, '--no-save', '--RS-conf',
         os.path.join(HERE_PATH, 'test.conf'),
         '--RS-port', str(RPORT)],
        stdout=open('/dev/null'), stderr=subprocess.PIPE)
    # wait a moment until Rserve starts listening on RPORT
    time.sleep(0.6)
    if rProc.poll():
        # process has already terminated, so provide its output on stderr
        # to the user:
        raise RuntimeError('Rserve has terminated prematurely with the '
                           'following message:  %s' % rProc.stderr.read())

    # store original socket timeout and set timeout to new value during startup
    # of Rserve:
    defaultTimeout = socket.getdefaulttimeout()
    socket.setdefaulttimeout(1)

    rserv = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
    cnt = 0
    # give it a maximum of 10 tries with some sleep in between to wait for
    # Rserve to come into action!
    while cnt < 10:
        try:
            # open a socket connection to Rserve
            rserv.connect(('', RPORT))
        except socket.error:
            time.sleep(0.3)
            cnt += 1
        else:
            # got a connection! Jump out of the loop
            break
    else:
        # after trying 10 times we still got no connection to Rserv - something
        # must be wrong.
        raise RuntimeError('Could not connect to Rserv over the network')

    # set back original default timeout value:
    socket.setdefaulttimeout(defaultTimeout)

    # make a simple test that Rserve really answers correctly by looking at the
    # first few bytes:
    hdr = rserv.recv(1024)
    rserv.close()
    if not hdr.startswith(b'Rsrv01'):
        rProc.terminate()
        raise RuntimeError('received wrong header information from '
                           'socket (was: "%s")' % hdr[:10])
    return rProc

def compareArrays(arr1, arr2):
    """Compare two (possibly nested) arrays"""
    def _compareArrays(xarr1, xarr2):
        assert xarr1.shape == xarr2.shape
        for idx in range(len(xarr1)):
            if isinstance(xarr1[idx], ndarray):
                _compareArrays(xarr1[idx], xarr2[idx])
            else:
                if type(xarr1[idx]) in [float, float32, float64, complex,
                                        complex64, complex128]:
                    # make a comparison which works for floats and complex
                    # numbers
                    assert abs(xarr1[idx] - xarr2[idx]) < 0.000001
                else:
                    assert xarr1[idx] == xarr2[idx]
    try:
        _compareArrays(arr1, arr2)
    except TypeError:
        return False
    return True

max_stars_repo_name: flying-sheep/pyRserve

max_issues_repo_name: flying-sheep/pyRserve

max_forks_repo_name: flying-sheep/pyRserve


text:
# Copyright 2020 Pulser Development Team
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import json
from unittest.mock import patch

import numpy as np
import pytest

import pulser
from pulser import Sequence, Pulse, Register
from pulser.devices import Chadoq2, MockDevice
from pulser.devices._device_datacls import Device
from pulser.sequence import _TimeSlot
from pulser.waveforms import (
    BlackmanWaveform,
    CompositeWaveform,
    RampWaveform,
    InterpolatedWaveform,
)

reg = Register.triangular_lattice(4, 7, spacing=5, prefix="q")
device = Chadoq2

def test_init():
    with pytest.raises(TypeError, match="must be of type 'Device'"):
        Sequence(reg, Device)

    fake_device = Device("fake", 2, 100, 100, 1, Chadoq2._channels)
    with pytest.warns(UserWarning, match="imported from 'pulser.devices'"):
        Sequence(reg, fake_device)

    seq = Sequence(reg, device)
    assert seq.qubit_info == reg.qubits
    assert seq.declared_channels == {}
    assert seq.available_channels.keys() == device.channels.keys()

def test_channel_declaration():
    seq = Sequence(reg, device)
    available_channels = set(seq.available_channels)
    seq.declare_channel("ch0", "rydberg_global")
    seq.declare_channel("ch1", "raman_local")
    with pytest.raises(ValueError, match="No channel"):
        seq.declare_channel("ch2", "raman")
    with pytest.raises(ValueError, match="not available"):
        seq.declare_channel("ch2", "rydberg_global")
    with pytest.raises(ValueError, match="name is already in use"):
        seq.declare_channel("ch0", "raman_local")

    chs = {"rydberg_global", "raman_local"}
    assert seq._schedule["ch0"][-1] == _TimeSlot(
        "target", -1, 0, set(seq.qubit_info.keys())
    )
    assert set(seq.available_channels) == available_channels - chs

    seq2 = Sequence(reg, MockDevice)
    available_channels = set(seq2.available_channels)
    seq2.declare_channel("ch0", "raman_local", initial_target="q1")
    seq2.declare_channel("ch1", "rydberg_global")
    seq2.declare_channel("ch2", "rydberg_global")
    assert set(seq2.available_channels) == available_channels - {"mw_global"}
    assert seq2._taken_channels == {
        "ch0": "raman_local",
        "ch1": "rydberg_global",
        "ch2": "rydberg_global",
    }
    assert seq2._taken_channels.keys() == seq2._channels.keys()
    with pytest.raises(ValueError, match="type 'Microwave' cannot work "):
        seq2.declare_channel("ch3", "mw_global")

    seq2 = Sequence(reg, MockDevice)
    seq2.declare_channel("ch0", "mw_global")
    assert set(seq2.available_channels) == {"mw_global"}
    with pytest.raises(
        ValueError,
        match="cannot work simultaneously with the declared 'Microwave'",
    ):
        seq2.declare_channel("ch3", "rydberg_global")

def test_target():
    seq = Sequence(reg, device)
    seq.declare_channel("ch0", "raman_local", initial_target="q1")
    seq.declare_channel("ch1", "rydberg_global")

    with pytest.raises(ValueError, match="name of a declared channel"):
        seq.target("q0", "ch2")
    with pytest.raises(ValueError, match="qubits must belong"):
        seq.target(0, "ch0")
        seq.target("0", "ch1")
    with pytest.raises(ValueError, match="Can only choose target of 'Local'"):
        seq.target("q3", "ch1")
    with pytest.raises(ValueError, match="can target at most 1 qubits"):
        seq.target(["q1", "q5"], "ch0")

    assert seq._schedule["ch0"][-1] == _TimeSlot("target", -1, 0, {"q1"})
    seq.target("q4", "ch0")
    retarget_t = seq.declared_channels["ch0"].retarget_time
    assert seq._schedule["ch0"][-1] == _TimeSlot(
        "target", 0, retarget_t, {"q4"}
    )
    seq.target("q4", "ch0")  # targets the same qubit
    seq.target("q20", "ch0")
    assert seq._schedule["ch0"][-1] == _TimeSlot(
        "target", retarget_t, 2 * retarget_t, {"q20"}
    )
    seq.delay(216, "ch0")
    seq.target("q2", "ch0")
    ti = 2 * retarget_t + 216
    tf = ti + 16
    assert seq._schedule["ch0"][-1] == _TimeSlot("target", ti, tf, {"q2"})

    seq.delay(220, "ch0")
    seq.target("q1", "ch0")
    ti = tf + 220
    tf = ti
    assert seq._schedule["ch0"][-1] == _TimeSlot("target", ti, tf, {"q1"})

    seq.delay(100, "ch0")
    seq.target("q10", "ch0")
    ti = tf + 100
    tf = ti + 120
    assert seq._schedule["ch0"][-1] == _TimeSlot("target", ti, tf, {"q10"})

    seq2 = Sequence(reg, MockDevice)
    seq2.declare_channel("ch0", "raman_local", initial_target={"q1", "q10"})
    seq2.phase_shift(1, "q2")
    with pytest.raises(ValueError, match="qubits with different phase"):
        seq2.target({"q3", "q1", "q2"}, "ch0")

def test_delay():
    seq = Sequence(reg, device)
    seq.declare_channel("ch0", "raman_local")
    with pytest.raises(ValueError, match="Use the name of a declared channel"):
        seq.delay(1e3, "ch01")
    with pytest.raises(ValueError, match="channel has no target"):
        seq.delay(100, "ch0")
    seq.target("q19", "ch0")
    seq.delay(388, "ch0")
    assert seq._last("ch0") == _TimeSlot("delay", 0, 388, {"q19"})

def test_delay_min_duration():
    # Check that a delay shorter than a channel's minimal duration
    # is automatically extended to that minimal duration
    seq = Sequence(reg, device)
    seq.declare_channel("ch0", "rydberg_global")
    seq.declare_channel("ch1", "rydberg_local")
    seq.target("q0", "ch1")
    pulse0 = Pulse.ConstantPulse(52, 1, 1, 0)
    pulse1 = Pulse.ConstantPulse(180, 1, 1, 0)
    seq.add(pulse1, "ch1")
    seq.add(pulse0, "ch0")
    seq.target("q1", "ch1")
    seq.add(pulse1, "ch1")
    min_duration = seq._channels["ch1"].min_duration
    assert seq._schedule["ch1"][3] == _TimeSlot(
        "delay", 220, 220 + min_duration, {"q1"}
    )

def test_min_pulse_duration():
    seq = Sequence(reg, device)
    seq.declare_channel("ch0", "rydberg_global")
    seq.declare_channel("ch1", "rydberg_local")
    seq.target("q0", "ch1")
    pulse0 = Pulse.ConstantPulse(60, 1, 1, 0)
    pulse1 = Pulse.ConstantPulse(80, 1, 1, 0)
    seq.add(pulse1, "ch1")
    assert seq._min_pulse_duration() == 80
    seq.add(pulse0, "ch0")
    seq.delay(52, "ch0")
    seq.target("q1", "ch1")
    seq.add(pulse1, "ch1")
    assert seq._min_pulse_duration() == 60

def test_phase():
    seq = Sequence(reg, device)
    seq.declare_channel("ch0", "raman_local", initial_target="q0")
    seq.phase_shift(-1, "q0", "q1")
    with pytest.raises(ValueError, match="id of a qubit declared"):
        seq.current_phase_ref(0, "digital")
    with pytest.raises(ValueError, match="targets the given 'basis'"):
        seq.current_phase_ref("q1", "ground-rydberg")
    with pytest.raises(ValueError, match="No declared channel targets"):
        seq.phase_shift(1, "q3", basis="hyperfine")
    assert seq.current_phase_ref("q0", "digital") == 2 * np.pi - 1

    # Phase shifts of 0
    seq.phase_shift(0, "q0")
    seq.phase_shift(-8 * np.pi, "q1")
    assert seq.current_phase_ref("q0", "digital") == 2 * np.pi - 1
    assert seq.current_phase_ref("q1", "digital") == 2 * np.pi - 1

    with pytest.raises(ValueError, match="targets have to be qubit ids"):
        seq.phase_shift(np.pi, "q1", "q4", "q100")

    seq.declare_channel("ch1", "rydberg_global")
    seq.phase_shift(1, *seq._qids, basis="ground-rydberg")
    for q in seq.qubit_info:
        assert seq.current_phase_ref(q, "ground-rydberg") == 1
    seq.phase_shift(1, *seq._qids)
    assert seq.current_phase_ref("q1", "digital") == 0
    assert seq.current_phase_ref("q10", "digital") == 1

def test_align():
    seq = Sequence(reg, device)
    seq.declare_channel("ch0", "raman_local", initial_target="q0")
    seq.declare_channel("ch1", "rydberg_global")
    with pytest.raises(ValueError, match="names must correspond to declared"):
        seq.align("ch0", "ch1", "ch2")
    with pytest.raises(ValueError, match="more than once"):
        seq.align("ch0", "ch1", "ch0")
    with pytest.raises(ValueError, match="at least two channels"):
        seq.align()
        seq.align("ch1")

def test_measure():
    pulse = Pulse.ConstantPulse(500, 2, -10, 0, post_phase_shift=np.pi)
    seq = Sequence(reg, MockDevice)
    seq.declare_channel("ch0", "rydberg_global")
    assert "XY" in MockDevice.supported_bases
    with pytest.raises(ValueError, match="not supported"):
        seq.measure(basis="XY")
    seq.measure()
    with pytest.raises(RuntimeError, match="already been measured"):
        seq.measure(basis="digital")
    with pytest.raises(RuntimeError, match="Nothing more can be added."):
        seq.add(pulse, "ch0")

    seq = Sequence(reg, MockDevice)
    seq.declare_channel("ch0", "mw_global")
    assert "digital" in MockDevice.supported_bases
    with pytest.raises(ValueError, match="not supported"):
        seq.measure(basis="digital")
    seq.measure(basis="XY")

def test_str():
    seq = Sequence(reg, device)
    seq.declare_channel("ch0", "raman_local", initial_target="q0")
    pulse = Pulse.ConstantPulse(500, 2, -10, 0, post_phase_shift=np.pi)
    seq.add(pulse, "ch0")
    seq.delay(200, "ch0")
    seq.target("q7", "ch0")
    seq.measure("digital")
    msg = (
        "Channel: ch0\nt: 0 | Initial targets: q0 | Phase Reference: 0.0 "
        + "\nt: 0->500 | Pulse(Amp=2 rad/µs, Detuning=-10 rad/µs, Phase=0) "
        + "| Targets: q0\nt: 500->700 | Delay \nt: 700->700 | Target: q7 | "
        + "Phase Reference: 0.0\n\nMeasured in basis: digital"
    )
    assert seq.__str__() == msg

def test_sequence():
    seq = Sequence(reg, device)
    assert seq.get_duration() == 0
    with pytest.raises(RuntimeError, match="empty sequence"):
        seq.draw()
    seq.declare_channel("ch0", "raman_local", initial_target="q0")
    seq.declare_channel("ch1", "rydberg_local", initial_target="q0")
    seq.declare_channel("ch2", "rydberg_global")
    assert seq.get_duration("ch0") == 0
    assert seq.get_duration("ch2") == 0
    seq.phase_shift(np.pi, "q0", basis="ground-rydberg")

    with patch("matplotlib.pyplot.show"):
        seq.draw()

    pulse1 = Pulse(
        InterpolatedWaveform(500, [0, 1, 0]),
        InterpolatedWaveform(500, [-1, 1, 0]),
        phase=0,
        post_phase_shift=np.pi,
    )
    pulse2 = Pulse.ConstantDetuning(
        BlackmanWaveform(1e3, np.pi / 4), 25, np.pi, post_phase_shift=1
    )
    with pytest.raises(TypeError):
        seq.add([1, 5, 3], "ch0")
    with pytest.raises(ValueError, match="amplitude goes over the maximum"):
        seq.add(
            Pulse.ConstantPulse(20, 2 * np.pi * 10, -2 * np.pi * 100, 0), "ch2"
        )
    with pytest.raises(
        ValueError, match="detuning values go out of the range"
    ):
        seq.add(
            Pulse.ConstantPulse(500, 2 * np.pi, -2 * np.pi * 100, 0), "ch0"
        )
    with pytest.raises(ValueError, match="qubits with different phase ref"):
        seq.add(pulse2, "ch2")
    with pytest.raises(ValueError, match="Invalid protocol"):
        seq.add(pulse1, "ch0", protocol="now")

    wf_ = CompositeWaveform(BlackmanWaveform(30, 1), RampWaveform(15, 0, 2))
    with pytest.raises(TypeError, match="Failed to automatically adjust"):
        with pytest.warns(UserWarning, match="rounded up to 48 ns"):
            seq.add(Pulse.ConstantAmplitude(1, wf_, 0), "ch0")

    pulse1_ = Pulse.ConstantPulse(499, 2, -10, 0, post_phase_shift=np.pi)
    with pytest.warns(UserWarning, match="rounded up to 500 ns"):
        seq.add(pulse1_, "ch0")
    seq.add(pulse1, "ch1")
    seq.add(pulse2, "ch2")

    assert seq._last("ch0").ti == 0
    assert seq._last("ch0").tf == seq._last("ch1").ti
    assert seq._last("ch2").tf == seq._last("ch2").ti + 1000
    assert seq.current_phase_ref("q0", "digital") == np.pi

    seq.add(pulse1, "ch2")
    assert seq.get_duration("ch2") == 2500
    seq.add(pulse2, "ch1", protocol="no-delay")
    assert seq.get_duration("ch1") == 3500
    seq.add(pulse1, "ch0", protocol="no-delay")
    assert seq._last("ch0").ti == 500
    assert seq.get_duration("ch0") == 1000
    assert seq.current_phase_ref("q0", "digital") == 0
    seq.phase_shift(np.pi / 2, "q1")
    seq.target("q1", "ch0")
    assert seq._last_used["digital"]["q1"] == 0
    assert seq._last_target["ch0"] == 1000
    assert seq._last("ch0").ti == 1000
    assert seq.get_duration("ch0") == 1000
    seq.add(pulse1, "ch0")
    assert seq._last("ch0").ti == 2500
    assert seq.get_duration("ch0") == 3000
    seq.add(pulse1, "ch0", protocol="wait-for-all")
    assert seq._last("ch0").ti == 3500
    assert seq.get_duration("ch2") != seq.get_duration("ch0")
    seq.align("ch0", "ch2")
    assert seq.get_duration("ch2") == seq.get_duration("ch0")

    with patch("matplotlib.pyplot.show"):
        seq.draw(draw_phase_shifts=True)

    assert seq.get_duration() == 4000

    seq.measure(basis="digital")

    with patch("matplotlib.pyplot.show"):
        seq.draw(draw_phase_area=True)

    s = seq.serialize()
    assert json.loads(s)["__version__"] == pulser.__version__
    seq_ = Sequence.deserialize(s)
    assert str(seq) == str(seq_)

max_stars_repo_name: rbstsai/Pulser

max_issues_repo_name: rbstsai/Pulser

max_forks_repo_name: rbstsai/Pulser


text:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import scipy.stats

def mean_confidence_interval(data, confidence=0.95):
    a = 1.0 * np.array(data)
    n = len(a)
    m, se = np.mean(a), scipy.stats.sem(a)
    h = se * scipy.stats.t.ppf((1 + confidence) / 2., n-1)
    return m, m-h, m+h

def PlotInfectionSubsystemBounds(tmin, CumIActual, CumDActual, INsim, IHsim, Dsim, file, SimLabel):
	obstime = np.linspace(tmin, len(CumIActual), len(CumIActual))
	simtime = np.linspace(tmin, len(INsim[:,0]), len(INsim[:,0]))
	CumIupper = np.zeros(len(INsim[:,0]))
	CumIlower = np.zeros(len(INsim[:,0]))
	CumImean = np.zeros(len(INsim[:,0]))
	Dupper = np.zeros(len(Dsim[:,0]))
	Dlower = np.zeros(len(Dsim[:,0]))
	Dmean = np.zeros(len(Dsim[:,0]))

	for i in range(len(INsim[:,0])):
		Dmean[i], Dlower[i], Dupper[i] = mean_confidence_interval(Dsim[i,:])
		CumImean[i], CumIlower[i], CumIupper[i] = mean_confidence_interval(IHsim[i,:] + INsim[i,:])

	plt.clf()
	plt.subplot(1, 1, 1)
	# plt.title('Infection R0 = ' + str(R0mean))
	plt.title(SimLabel)
	plt.plot(obstime, CumIActual, 'k', linewidth=2, label='Actual Cumulative Infections',color ='g')
	plt.plot(obstime, CumDActual, 'k', linewidth=2, label='Actual Deaths',color ='r')
	plt.plot(simtime, CumImean, 'k', linewidth=2, label='Model Cumulative Infections',color ='k')
	plt.fill_between(simtime, CumIupper, CumIlower, where= CumIupper > CumIlower, color='k', alpha=.25)
	plt.plot(simtime, Dmean, 'k', linewidth=2, label='Model Deaths',color ='b')
	plt.fill_between(simtime, Dupper, Dlower, where= Dupper > Dlower, color='b', alpha=.25)
	plt.xlabel('Days since 2020-01-22')
	plt.ylabel('Number of Individuals')
	plt.legend()
	plt.tight_layout()
	plt.savefig(''+str(file)+'/'+str(SimLabel)+'.png',dpi=300, bbox_inches='tight')
	# plt.show()

	# plt.clf()
	# plt.subplot(1, 1, 1)
	# # plt.title('Infection R0 = ' + str(R0mean))
	# plt.title(SimLabel)
	# plt.plot(obstime, CumIActual, 'k', linewidth=2, label='Actual Cumulative Infections',color ='g')
	# plt.plot(obstime, CumDActual, 'k', linewidth=2, label='Actual Deaths',color ='r')
	# plt.plot(simtime, IHsim[:,-1] + INsim[:,-1], 'k', linewidth=2, label='Model Cumulative Infections',color ='k')
	# plt.plot(simtime, Dsim[:,-1], 'k', linewidth=2, label='Model Deaths',color ='b')
	# plt.xlabel('Days since 2020-01-22')
	# plt.ylabel('Number of Individuals')
	# plt.legend()
	# plt.tight_layout()
	# plt.savefig(''+str(file)+'/BestFit'+str(SimLabel)+'.png',dpi=300, bbox_inches='tight')
	# plt.show()

max_stars_repo_name: yupengyanghuhu/USA-CDC-Bayesian-ODE

max_issues_repo_name: yupengyanghuhu/USA-CDC-Bayesian-ODE

max_forks_repo_name: yupengyanghuhu/USA-CDC-Bayesian-ODE

reason: import numpy,import scipy


text:
# This code snippet will load the referenced package and return a DataFrame.
# If the code is run in a PySpark environemnt, then the code will return a
# Spark DataFrame. If not, the code will return a Pandas DataFrame. You can
# copy this code snippet to another code file as needed.    C:\Users\zhouf\.amlenvrc.cmd

# Import python module

import pickle
import sys
import os

import pandas as pd
import numpy as np
import csv

from sklearn.preprocessing import LabelEncoder

from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
from sklearn.tree import DecisionTreeClassifier

from sklearn import grid_search
from sklearn import metrics

from sklearn.metrics import confusion_matrix
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score
from sklearn.cross_validation import train_test_split

from azureml.logging import get_azureml_logger

# initialize the logger

run_logger = get_azureml_logger()

# Step 1 - Data Preparation

# Create the outputs folder

os.makedirs('./outputs', exist_ok=True)

print('Python version: {}'.format(sys.version))
print()

# Load the bank dataset

df = pd.read_csv('./data/BankMarketCampaignTrainingSample.csv')

# Step 2 - Feature Engineering

# Encode columns

columns_to_encode = list(df.select_dtypes(include=['category','object']))
for column_to_encode in columns_to_encode:
    g = df.groupby(column_to_encode)
for column_to_encode in columns_to_encode:
    dummies = pd.get_dummies(df[column_to_encode])
    one_hot_col_names = []
    for col_name in list(dummies.columns):
        one_hot_col_names.append(column_to_encode + '_' + col_name)
    dummies.columns = one_hot_col_names
    df = df.drop(column_to_encode, axis=1)
    df = df.join(dummies)

# Keep only one response variable

df = df.drop('y_no', axis=1)

# Step 3 - Model Training and Evaluation

# Split Data for Training and Testing

train, test = train_test_split(df, test_size = 0.2, random_state=0)

# Specify the values of label and features in training and testing datasets.

train_y = train['y_yes'].values

train_x = train.drop('y_yes', axis=1)
train_x = train_x.values

test_y = test['y_yes'].values
test_x = test.drop('y_yes', axis=1)
test_x = test_x.values

# Initialize Sweep Parameter

# Logistic Regression

def make_sweep_parameter_lr_dict():
    '''
    :return: parameters of logistic regression model to tune
    '''
    return {"penalty": ("l1", "l2"), "C": [0.1, 1, 10]}

# Support Vector Machine

def make_sweep_parameter_svm_dict():
    '''
    :return: parameters of support vector machine model to tune
    '''
    return {"gamma": [0.0001, 0.01, 1, 100], "C": [1]}

# Decision Tree

def make_sweep_parameter_dt_dict():
    '''
    :return: parameters of decision tree model to tune
    '''
    return{'max_depth': [4, 8, 16, 32]}

# Sweep Parameters with Each Classifier using Cross Validation

def sweep_classifier_cross_val(train_feature, train_label, model_name, parameter_dict):
    '''
    :param train_feature: feature of training set
    :param train_label: label of training set
    :param model_name: model
    :param parameter_dict: parameters of the model to tune
    :return: trained model
    '''
    classifier = grid_search.GridSearchCV(model_name, parameter_dict)
    classifier.fit(train_feature, train_label)
    return classifier

# Validate Model using Testing Data

def prediction(classifier, test_feature):
    '''
    :param classifier: the trained classifier
    :param test_feature: feature of test set
    :return: scored label of test set
    '''
    predicted_label = classifier.predict(test_feature)
    return predicted_label

# Compute Metrics

def calc_metrics(test_label, predicted_label):
    '''
    :param test_label: label of test set
    :param predicted_label: scored label of test set
    :return: evaluation metrics
    '''
    cm = confusion_matrix(test_label, predicted_label)
    acc = accuracy_score(test_label, predicted_label)
    pre = precision_score(test_label, predicted_label)
    rec = recall_score(test_label, predicted_label)
    f1 = f1_score(test_label, predicted_label)
    return cm, acc, pre, rec, f1

# Put it all together

# Train and evaluate logistic regression

model_name = LogisticRegression()
parameter_dict = make_sweep_parameter_lr_dict()
lrf = sweep_classifier_cross_val(train_feature=train_x,
                                 train_label=train_y,
                                 model_name=model_name,
                                 parameter_dict=parameter_dict)

print('Logistic Regression Classifier:')
print(lrf)

predicted = prediction(lrf, test_x)

cm, acc, pre, rec, f1 = calc_metrics(test_y, predicted)

print('Confusion Matrix:')
print(cm)
print('Accuracy: {:.2f}'.format(acc))
print('Precision: {:.2f}'.format(pre))
print('Recall: {:.2f}'.format(rec))
print('F1: {:.2f}'.format(f1))

# Log lr model accuracy

run_logger.log("Logistic Regression Accuracy", acc)

# Train and evaluate support vector machine

model_name = SVC(kernel="rbf")
parameter_dict = make_sweep_parameter_svm_dict()
svf = sweep_classifier_cross_val(train_feature=train_x,
                                 train_label=train_y,
                                 model_name=model_name,
                                 parameter_dict=parameter_dict)

print('Support Vector Machine Classifier:')
print(svf)

predicted = prediction(svf, test_x)

cm, acc, pre, rec, f1 = calc_metrics(test_y, predicted)

print('Confusion Matrix:')
print(cm)
print('Accuracy: {:.2f}'.format(acc))
print('Precision: {:.2f}'.format(pre))
print('Recall: {:.2f}'.format(rec))
print('F1: {:.2f}'.format(f1))

# Log sv model accuracy

run_logger.log("SVM Accuracy", acc)

# Train and evaluate support vector machine

model_name = DecisionTreeClassifier(min_samples_split=20, random_state=0)
parameter_dict = make_sweep_parameter_dt_dict()
dtf = sweep_classifier_cross_val(train_feature=train_x,
                                 train_label=train_y,
                                 model_name=model_name,
                                 parameter_dict=parameter_dict)

print('Decision Tree Classifier:')
print(dtf)

predicted = prediction(dtf, test_x)

cm, acc, pre, rec, f1 = calc_metrics(test_y, predicted)

print('Confusion Matrix:')
print(cm)
print('Accuracy: {:.2f}'.format(acc))
print('Precision: {:.2f}'.format(pre))
print('Recall: {:.2f}'.format(rec))
print('F1: {:.2f}'.format(f1))

# Log dt model accuracy

run_logger.log("Decision Tree Accuracy", acc)

print("")
print("==========================================")
print("Serialize and deserialize using the outputs folder.")
print("")

# Serialize the model on disk in the special 'outputs' folder

print ("Export the lr model to lr.pkl")
f = open('./outputs/lr.pkl', 'wb')
pickle.dump(lrf, f)
f.close()

print ("Export the svf model to sv.pkl")
f = open('./outputs/sv.pkl', 'wb')
pickle.dump(svf, f)
f.close()

print ("Export the dtf model to dt.pkl")
f = open('./outputs/dt.pkl', 'wb')
pickle.dump(dtf, f)
f.close()

# Load the model back from the 'outputs' folder into memory

print("Import the model from lr.pkl")
f2 = open('./outputs/lr.pkl', 'rb')
lrf2 = pickle.load(f2)

print("Import the model from sv.pkl")
f2 = open('./outputs/sv.pkl', 'rb')
svf2 = pickle.load(f2)

print("Import the model from dt.pkl")
f2 = open('./outputs/dt.pkl', 'rb')
dtf2 = pickle.load(f2)


text:
import numpy as np

def add_null(data:np.ndarray, null_size=100):
    n_row, _ = data.shape
    n = (data > 0).sum(axis=0)
    pr = n / n_row
    pr_mu = pr.mean()
    pr_sigma = pr.std()

    null_pr = np.random.normal(size=null_size, loc=pr_mu, scale=pr_sigma)
    null_uni = np.random.random((n_row, null_size))
    null_exp = -np.log(np.random.random((n_row ,null_size)))
    null = null_exp
    null[null_uni >= null_pr] = 0
    return np.hstack((data, null))

max_stars_repo_name: IdlingBrainUT/NeuroScale

max_issues_repo_name: IdlingBrainUT/NeuroScale

max_forks_repo_name: IdlingBrainUT/NeuroScale