text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Author: Corentin Ravoux

Description : Utilities including coordinate conversion and various
interfaces with other codes.
"""

#############################################################################
#############################################################################
########################## MODULE IMPORTATION ###############################
#############################################################################
#############################################################################

import numpy as np
import logging, time, os
import fitsio
from scipy.ndimage.filters import gaussian_filter
from scipy import interpolate
from scipy.ndimage import map_coordinates
from scipy.optimize import leastsq
import multiprocessing as mp
from multiprocessing.reduction import ForkingPickler, AbstractReducer
import matplotlib.pyplot as plt
from scipy.stats import binned_statistic
from picca import constants

#############################################################################
#############################################################################
############################### FUNCTIONS ###################################
#############################################################################
#############################################################################

lambdaLy = 1215.673123130217

def mpc_per_pixel(size,shape):
    return(np.array(size)/(np.array(shape)-1))

def pixel_per_mpc(size,shape):
    return((np.array(shape)-1)/np.array(size))

def get_map_shape(size,mpc_per_pixel):
    return((np.array(size)/np.array(mpc_per_pixel)) + 1)

def get_map_size(shape,mpc_per_pixel):
    return((np.array(shape) - 1)*mpc_per_pixel)

def get_cosmo_function(Omega_m,Omega_k=0.):

    Cosmo = constants.Cosmo(Omega_m,Ok=Omega_k)
    rcomov = Cosmo.get_r_comov
    distang = Cosmo.get_dist_m
    redshift_array = np.linspace(0,5,10000)
    R_array = rcomov(redshift_array)
    Dm_array = rcomov(redshift_array)
    inv_rcomov = interpolate.interp1d(R_array,redshift_array)
    inv_distang = interpolate.interp1d(Dm_array,redshift_array)
    return(rcomov,distang,inv_rcomov,inv_distang)

def return_suplementary_parameters(mode,property=None,zmin=None,zmax=None):
    if(mode == "middle"):
        if(property is not None):
            zmin = property.boundary_sky_coord[0][2]
            zmax = property.boundary_sky_coord[1][2]
            suplementary_parameters = [(zmin + zmax)/2]
        elif((zmin is not None)&(zmin is not None)):
            suplementary_parameters = [(zmin + zmax)/2]
    else:
        suplementary_parameters = None
    return(suplementary_parameters)

def convert_cartesian_to_sky(X,Y,Z,method,inv_rcomov=None,inv_distang=None,distang=None,suplementary_parameters=None):
    """Mpc.h-1 to radians"""
    if(method == "full_angle"):
        (RA,DEC,z)= convert_cartesian_to_sky_full_angle(X,Y,Z,inv_rcomov)
    if(method == "full"):
        (RA,DEC,z)= convert_cartesian_to_sky_full(X,Y,Z,inv_rcomov)
    elif(method == "middle"):
        (RA,DEC,z)= convert_cartesian_to_sky_middle(X,Y,Z,inv_rcomov,distang,suplementary_parameters[0])
    return(RA,DEC,z)

def convert_sky_to_cartesian(RA,DEC,z,method,rcomov=None,distang=None,suplementary_parameters=None):
    """Radians to Mpc.h-1"""
    if(method == "full_angle"):
        (X,Y,Z)= convert_sky_to_cartesian_full_angle(RA,DEC,z,rcomov)
    elif(method == "full"):
        (X,Y,Z)= convert_sky_to_cartesian_full(RA,DEC,z,rcomov)
    elif(method == "middle"):
        (X,Y,Z)= convert_sky_to_cartesian_middle(RA,DEC,z,rcomov,distang,suplementary_parameters[0])
    return(X,Y,Z)

def convert_cartesian_to_sky_full_angle(X,Y,Z,inv_rcomov):
    RA = np.arctan2(X,Z)
    DEC = np.arcsin(Y/np.sqrt(X**2 + Y**2 + Z**2))
    z = inv_rcomov(np.sqrt(X**2 + Y**2 + Z**2))
    return(RA,DEC,z)

def convert_sky_to_cartesian_full_angle(RA,DEC,z,rcomov):
    X = rcomov(z)*np.sin(RA)*np.cos(DEC)
    Y = rcomov(z)*np.sin(DEC)
    Z = rcomov(z)*np.cos(RA)*np.cos(DEC)
    return(X,Y,Z)

def convert_cartesian_to_sky_full(X,Y,Z,inv_rcomov):
    z = inv_rcomov(np.sqrt(X**2 + Y**2 + Z**2))
    RA = X/Z
    DEC = Y/Z
    return(RA,DEC,z)

def convert_sky_to_cartesian_full(RA,DEC,z,rcomov):
    X = rcomov(z)*RA
    Y = rcomov(z)*DEC
    Z = rcomov(z)
    return(X,Y,Z)

def convert_cartesian_to_sky_middle(X,Y,Z,inv_rcomov,distang,middle_z):
    RA = X/distang(middle_z)
    DEC = Y/distang(middle_z)
    z = inv_rcomov(Z)
    return(RA,DEC,z)

def convert_sky_to_cartesian_middle(RA,DEC,z,rcomov,distang,middle_z):
    X = distang(middle_z)*RA
    Y = distang(middle_z)*DEC
    Z = rcomov(z)
    return(X,Y,Z)

def convert_z_cartesian_to_sky_middle(Z,inv_rcomov):
    z = inv_rcomov(Z)
    return(z)

def convert_z_sky_to_cartesian_middle(z,rcomov):
    Z = rcomov(z)
    return(Z)

def get_direction_indexes(direction,rotate):
    if (direction.lower() == "x")|(direction.lower() == "ra"):
        x_index, y_index, index_direction = 2, 1, 0
    elif (direction.lower() == "y")|(direction.lower() == "dec"):
        x_index, y_index, index_direction = 2, 0, 1
    elif (direction.lower() == "z")|(direction.lower() == "redshift"):
        x_index, y_index, index_direction = 1, 0, 2
    if(rotate): x_index,y_index = y_index,x_index
    index_dict = {0:"x",1:"y",2:"z","x_lab":"X","y_lab":"Y","z_lab":"Z"}
    return(x_index,y_index,index_direction,index_dict)

saclay_mock_lines_per_box = {"box" : 1,"vx" : 5,"vy" : 5,"vz" : 5,"eta_xx":1,"eta_xy":1,"eta_xz":1,"eta_yx":1,"eta_yy":1,"eta_yz":1,"eta_zx":1,"eta_zy":1,"eta_zz":1}

def saclay_mock_box_cosmo_parameters(box_shape,size_cell):
    import cosmolopy.distance as dist
    try :
        from SaclayMocks import constant as saclay_mock_constant
    except :
        from lelantos.saclaymocks import constant as saclay_mock_constant
        print("SaclayMocks might be updated, we suggest to install SaclayMocks independently")
    NZ = box_shape[2]
    DZ = size_cell
    LZ = NZ*DZ

    h = saclay_mock_constant.h
    Om = saclay_mock_constant.omega_M_0
    OL = saclay_mock_constant.omega_lambda_0
    Ok = saclay_mock_constant.omega_k_0
    z0 = saclay_mock_constant.z0

    cosmo_fid = {'omega_M_0':Om, 'omega_lambda_0':OL, 'omega_k_0':Ok, 'h':h}
    R_of_z, z_of_R = dist.quick_distance_function(dist.comoving_distance, return_inverse=True, **cosmo_fid)
    R0 = h * R_of_z(z0)
    Rmin = R0 - LZ/2
    Rmax = R0 + LZ/2
    return(R0,z0,R_of_z,z_of_R,Rmin,Rmax,h)

def saclay_mock_center_of_the_box(box_bound):
    ra0_box = (box_bound[0] + box_bound[1])/2
    dec0_box = (box_bound[2] + box_bound[3])/2
    return(ra0_box,dec0_box)

def saclay_mock_coord_dm_map(X,Y,Z,Rmin,size_cell,box_shape,interpolation_method):
    size_cell = size_cell
    center_x = (box_shape[0]-1) / 2
    center_y = (box_shape[1]-1) / 2
    if(interpolation_method.upper() == "NEAREST"):
        n_i = (np.round(center_x + X/size_cell,0)).astype(int)
        n_j = (np.round(center_y + Y/size_cell,0)).astype(int)
        n_k = np.round((Z - Rmin)/size_cell,0).astype(int)
    else:
        n_i = center_x + X/size_cell
        n_j = center_y + Y/size_cell
        n_k = (Z - Rmin)/size_cell
    return(n_i,n_j,n_k)

def saclay_mock_read_box(box_dir,n_x,name_box):
    name = "{}-{}.fits".format(name_box,str(n_x))
    box = fitsio.FITS(os.path.join(box_dir, name))[0][:,:,:][0]
    return(box)

def saclay_mock_get_box(box_dir,box_shape,name_box="box"):
    line_per_box = saclay_mock_lines_per_box[name_box]
    DM_mocks = np.zeros((box_shape[0],box_shape[1],box_shape[2]))
    for i in range(box_shape[0]//line_per_box):
        DM_mocks[i*line_per_box:(i+1)*line_per_box,:,:] = saclay_mock_read_box(box_dir,i,name_box)[:,:]
    return(DM_mocks)

def saclay_mock_sky_to_cartesian(ra,dec,R,ra0,dec0):
    '''
    XYZ of a point P (ra,dec,R) in a frame with
    observer at O, Z along OP, X along ra0, Y along dec0
    angles in radians
    tested that ra,dec, R = box.ComputeRaDecR(R0,ra0,dec0,X,Y,Z)
    x,y,z = box.ComputeXYZ(ra[0],dec[0],R,ra0,dec0)
    print x-X,y-Y,z-R0-Z        prints ~1E-13  for random inputs
    '''
    try :
        from SaclayMocks import box as saclay_mock_box
    except :
        from lelantos.saclaymocks import box as saclay_mock_box
        print("SaclayMocks might be updated, we suggest to install SaclayMocks independently")
    X,Y,Z = saclay_mock_box.ComputeXYZ2(ra*(np.pi/180),dec*(np.pi/180),R,ra0*(np.pi/180),dec0*(np.pi/180))
    return(X,Y,Z)

def cut_sky_catalog(ra,dec,z,ramin=None,ramax=None,decmin=None,decmax=None,zmin=None,zmax=None):
    mask = np.full(ra.shape,True)
    if(ramin is not None):
        mask &= ra > np.radians(ramin)
    if(ramax is not None):
        mask &= ra < np.radians(ramax)
    if(decmin is not None):
        mask &= dec > np.radians(decmin)
    if(decmax is not None):
        mask &= dec < np.radians(decmax)
    if(zmin is not None):
        mask &= z > zmin
    if(zmax is not None):
        mask &= z < zmax
    return(mask)

def init_shared_array(shape,full_value=np.inf):
    distance_array = np.full(shape,full_value)
    shared_arr = mp.Array('d', distance_array.flatten())
    del distance_array
    return(shared_arr)

def mp_array_to_numpyarray(mp_arr):
    return np.frombuffer(mp_arr.get_obj())

def bin_ndarray(ndarray, new_shape, operation='mean'):
    """
    From : https://stackoverflow.com/questions/8090229/resize-with-averaging-or-rebin-a-numpy-2d-array/29042041
    Bins an ndarray in all axes based on the target shape, by summing or
    averaging.
    Number of output dimensions must match number of input dimensions.
    Example
    -------
    >>> m = np.arange(0,100,1).reshape((10,10))
    >>> n = bin_ndarray(m, new_shape=(5,5), operation='sum')
    >>> print(n)
    [[ 22  30  38  46  54]
    [102 110 118 126 134]
    [182 190 198 206 214]
    [262 270 278 286 294]
    [342 350 358 366 374]]
    """
    if not operation.lower() in ['sum', 'mean', 'average', 'avg','gauss']:
        raise ValueError("Operation {} not supported.".format(operation))
    if ndarray.ndim != len(new_shape):
        raise ValueError("Shape mismatch: {} -> {}".format(ndarray.shape,new_shape))
    compression_pairs = [(d, c//d) for d, c in zip(new_shape,ndarray.shape)]
    flattened = [l for p in compression_pairs for l in p]
    ndarray = ndarray.reshape(flattened)
    for i in range(len(new_shape)):
        if operation.lower() == "sum":
            ndarray = ndarray.sum(-1*(i+1))
        elif operation.lower() in ["mean", "average", "avg"]:
            ndarray = ndarray.mean(-1*(i+1))
        elif operation.lower() in ["gauss"]:
            if i!=0 : raise KeyError("gaussian mean is not available for dim higher than 1")
            from scipy import signal
            newndarray = np.zeros(new_shape)
            gaussian_weights = signal.gaussian(int(ndarray.shape[1]),int(ndarray.shape[1])/4)
            for j in range(len(ndarray)):
                newndarray[j] = np.average(ndarray[j],axis=0,weights=gaussian_weights)
            ndarray =newndarray
    return ndarray

def interpolate_map(interpolation_method,map_array,coord):
    if(interpolation_method.upper() == "NEAREST"):
        coord = np.around(coord,decimals=0).astype(int)
        DM_map = map_array[coord[:,:,:,0],coord[:,:,:,1],coord[:,:,:,2]]
    elif(interpolation_method.upper() == "LINEAR"):
        points =  coord.reshape(coord.shape[0]*coord.shape[1]*coord.shape[2],3)
        DM_map = map_coordinates(map_array, np.transpose(points), order=1).reshape((coord.shape[0],coord.shape[1],coord.shape[2]))
    elif(interpolation_method.upper() == "SPLINE"):
        points =  coord.reshape(coord.shape[0]*coord.shape[1]*coord.shape[2],3)
        DM_map = map_coordinates(map_array, np.transpose(points), order=2).reshape((coord.shape[0],coord.shape[1],coord.shape[2]))
    else :
        raise ValueError("Please select NEAREST, LINEAR or SPLINE as interpolation_method")
    return(DM_map)

def interpolate_and_fill_map(interpolation_method,map_array,coord):
    if(interpolation_method.upper() == "NEAREST"):
        coord_nearest = np.around(coord,decimals=0).astype(int)
        map_to_fill = map_array[coord_nearest[:,0],coord_nearest[:,1],coord_nearest[:,2]]
        del coord_nearest
    elif(interpolation_method.upper() == "LINEAR"):
        map_to_fill = map_coordinates(map_array, np.transpose(coord), order=1)
    elif(interpolation_method.upper() == "SPLINE"):
        map_to_fill = map_coordinates(map_array, np.transpose(coord), order=2)
    else :
        raise ValueError("Please select NEAREST, LINEAR or SPLINE as interpolation_method")
    return(map_to_fill)

def gaussian_smoothing(mapdata,sigma):
    gaussian_map = gaussian_filter(mapdata,sigma)
    return(gaussian_map)

class gaussian_fitter_2d(object):

    def __init__(self,inpdata=None):

        self.inpdata = inpdata

    def moments2D(self):
        """ Returns the (amplitude, xcenter, ycenter, xsigma, ysigma, rot, bkg, e) estimated from moments in the 2d input array Data """

        bkg=np.median(np.hstack((self.inpdata[0,:],self.inpdata[-1,:],self.inpdata[:,0],self.inpdata[:,-1])))  #Taking median of the 4 edges points as background
        Data=np.ma.masked_less(self.inpdata-bkg,0)   #Removing the background for calculating moments of pure 2D gaussian
        #We also masked any negative values before measuring moments

        amplitude=Data.max()

        total= float(Data.sum())
        Xcoords,Ycoords= np.indices(Data.shape)

        xcenter= (Xcoords*Data).sum()/total
        ycenter= (Ycoords*Data).sum()/total

        RowCut= Data[int(xcenter),:]  # Cut along the row of data near center of gaussian
        ColumnCut= Data[:,int(ycenter)]  # Cut along the column of data near center of gaussian
        xsigma= np.sqrt(np.ma.sum(ColumnCut* (np.arange(len(ColumnCut))-xcenter)**2)/ColumnCut.sum())
        ysigma= np.sqrt(np.ma.sum(RowCut* (np.arange(len(RowCut))-ycenter)**2)/RowCut.sum())

        #Ellipcity and position angle calculation
        Mxx= np.ma.sum((Xcoords-xcenter)*(Xcoords-xcenter) * Data) /total
        Myy= np.ma.sum((Ycoords-ycenter)*(Ycoords-ycenter) * Data) /total
        Mxy= np.ma.sum((Xcoords-xcenter)*(Ycoords-ycenter) * Data) /total
        e= np.sqrt((Mxx - Myy)**2 + (2*Mxy)**2) / (Mxx + Myy)
        pa= 0.5 * np.arctan(2*Mxy / (Mxx - Myy))
        rot= np.rad2deg(pa)

        return amplitude,xcenter,ycenter,xsigma,ysigma, rot,bkg, e

    def Gaussian2D(self,amplitude, xcenter, ycenter, xsigma, ysigma, rot,bkg):
        """ Returns a 2D Gaussian function with input parameters. rotation input rot should be in degress """
        rot=np.deg2rad(rot)  #Converting to radians
        Xc=xcenter*np.cos(rot) - ycenter*np.sin(rot)  #Centers in rotated coordinates
        Yc=xcenter*np.sin(rot) + ycenter*np.cos(rot)
        #Now lets define the 2D gaussian function
        def Gauss2D(x,y) :
            """ Returns the values of the defined 2d gaussian at x,y """
            xr=x * np.cos(rot) - y * np.sin(rot)  #X position in rotated coordinates
            yr=x * np.sin(rot) + y * np.cos(rot)
            return amplitude*np.exp(-(((xr-Xc)/xsigma)**2 +((yr-Yc)/ysigma)**2)/2) +bkg

        return Gauss2D

    def FitGauss2D(self,ip=None):
        """ Fits 2D gaussian to Data with optional Initial conditions ip=(amplitude, xcenter, ycenter, xsigma, ysigma, rot, bkg)
        Example:
        >>> X,Y=np.indices((40,40),dtype=np.float)
        >>> Data=np.exp(-(((X-25)/5)**2 +((Y-15)/10)**2)/2) + 1
        >>> FitGauss2D(Data)
        (array([  1.00000000e+00,   2.50000000e+01,   1.50000000e+01, 5.00000000e+00,   1.00000000e+01,   2.09859373e-07, 1]), 2)
         """
        if ip is None:   #Estimate the initial parameters form moments and also set rot angle to be 0
            ip=self.moments2D()[:-1]   #Remove ellipticity from the end in parameter list

        Xcoords,Ycoords= np.indices(self.inpdata.shape)
        def errfun(ip):
            dXcoords= Xcoords-ip[1]
            dYcoords= Ycoords-ip[2]
            Weights=np.sqrt(np.square(dXcoords)+np.square(dYcoords)) # Taking radius as the weights for least square fitting
            return np.ravel((self.Gaussian2D(*ip)(*np.indices(self.inpdata.shape)) - self.inpdata)/np.sqrt(Weights))  #Taking a sqrt(weight) here so that while scipy takes square of this array it will become 1/r weight.

        p, success = leastsq(errfun, ip)

        return p,success

def create_log(log_level="info"):
    log = Logger(log_level=log_level)
    log.setup_logging()
    return(log)

def create_report_log(name="Python_Report",log_level="info"):
    log = Logger(name=name,log_level=log_level)
    log.setup_report_logging()
    return(log)

_logging_handler = None

class Logger(object):

    def __init__(self,name="Python_Report",log_level="info"):
        self.name = name
        self.log_level = log_level

    def setup_logging(self):
    	"""
        Taken from https://nbodykit.readthedocs.io/
    	Turn on logging, with the specified level.
    	Parameters
    	----------
    	log_level : 'info', 'debug', 'warning'
    		the logging level to set; logging below this level is ignored
    	"""

    	# This gives:
    	#
    	# [ 000000.43 ]   0: 06-28 14:49  measurestats	INFO	 Nproc = [2, 1, 1]
    	# [ 000000.43 ]   0: 06-28 14:49  measurestats	INFO	 Rmax = 120

    	levels = {"info" : logging.INFO,"debug" : logging.DEBUG,"warning" : logging.WARNING}

    	logger = logging.getLogger();
    	t0 = time.time()

    	class Formatter(logging.Formatter):
    		def format(self, record):
    			s1 = ('[ %09.2f ]: ' % (time.time() - t0))
    			return s1 + logging.Formatter.format(self, record)

    	fmt = Formatter(fmt='%(asctime)s %(name)-15s %(levelname)-8s %(message)s',
    					datefmt='%m-%d %H:%M ')

    	global _logging_handler
    	if _logging_handler is None:
    		_logging_handler = logging.StreamHandler()
    		logger.addHandler(_logging_handler)

    	_logging_handler.setFormatter(fmt)
    	logger.setLevel(levels[self.log_level])

    def setup_report_logging(self):
        levels = {"info" : logging.INFO,"debug" : logging.DEBUG,"warning" : logging.WARNING}
        logging.basicConfig(filename=self.name, filemode='w',level=levels[self.log_level],format='%(asctime)s :: %(levelname)s :: %(message)s')

    @staticmethod
    def add(line,level="info"):
        if(level=="info"):
            logging.info(line)
        if(level=="warning"):
            logging.warning(line)
        if(level=="debug"):
            logging.debug(line)

    @staticmethod
    def add_array_statistics(arr,char):
        if(arr is not None):
            Logger.add(f"Min of {char}: {arr.min()}")
            Logger.add(f"Max of {char}: {arr.max()}")
            Logger.add(f"Mean of {char}: {arr.mean()}")
            Logger.add(f"Standard deviation of {char}: {arr.std()}")

    @staticmethod
    def close():
        logging.shutdown()

def latex_float(float_input,decimals_input="{0:.2g}"):
    """
    example use:
    import matplotlib.pyplot as plt
    plt.figure(),plt.clf()
    plt.plot(np.array([1,2.]),'ko-',label="$P_0="+latex_float(7.63e-5)+'$'),
    plt.legend()
    """
    float_str = decimals_input.format(float_input)
    if "e" in float_str:
        base, exponent = float_str.split("e")
        return r"{0} \times 10^{{{1}}}".format(base, int(exponent))
    else:
        return float_str

def return_key(dictionary,string,default_value):
    return(dictionary[string] if string in dictionary.keys() else default_value)

### https://stackoverflow.com/questions/51562221/python-multiprocessing-overflowerrorcannot-serialize-a-bytes-object-larger-t/51568084
### https://stackoverflow.com/questions/47776486/python-struct-error-i-format-requires-2147483648-number-2147483647/47776649#47776649
### Solve the pickle issue OverflowError('cannot serialize a bytes object larger than 4GiB')
### When too large array are considered in multiprocessing
### Add the following lines to you code:
### ctx = mp.get_context()
### ctx.reducer = utils.Pickle4Reducer()
### utils.patch_mp_connection_bpo_17560(log = log)
### You might want to change multiprocessing context with: mp.get_context('fork')

class ForkingPickler4(ForkingPickler):
    def __init__(self, *args):
        if len(args) > 1:
            args[1] = 2
        else:
            args.append(2)
        super().__init__(*args)

    @classmethod
    def dumps(cls, obj, protocol=4):
        return ForkingPickler.dumps(obj, protocol)

def dump(obj, file, protocol=4):
    ForkingPickler4(file, protocol).dump(obj)

class Pickle4Reducer(AbstractReducer):
    ForkingPickler = ForkingPickler4
    register = ForkingPickler4.register
    dump = dump

import functools
import struct
import sys

def patch_mp_connection_bpo_17560(log = None):
    """Apply PR-10305 / bpo-17560 connection send/receive max size update

    See the original issue at https://bugs.python.org/issue17560 and
    https://github.com/python/cpython/pull/10305 for the pull request.

    This only supports Python versions 3.3 - 3.7, this function
    does nothing for Python versions outside of that range.

    """
    patchname = "Multiprocessing connection patch for bpo-17560"
    if not (3, 3) < sys.version_info < (3, 8):
        if log is not None:
            log.add(
            patchname + " not applied, not an applicable Python version: %s",
            sys.version
        )
        return

    from multiprocessing.connection import Connection

    orig_send_bytes = Connection._send_bytes
    orig_recv_bytes = Connection._recv_bytes
    if (
        orig_send_bytes.__code__.co_filename == __file__
        and orig_recv_bytes.__code__.co_filename == __file__
    ):
        if log is not None:
            log.add(patchname + " already applied, skipping")
        return

    @functools.wraps(orig_send_bytes)
    def send_bytes(self, buf):
        n = len(buf)
        if n > 0x7fffffff:
            pre_header = struct.pack("!i", -1)
            header = struct.pack("!Q", n)
            self._send(pre_header)
            self._send(header)
            self._send(buf)
        else:
            orig_send_bytes(self, buf)

    @functools.wraps(orig_recv_bytes)
    def recv_bytes(self, maxsize=None):
        buf = self._recv(4)
        size, = struct.unpack("!i", buf.getvalue())
        if size == -1:
            buf = self._recv(8)
            size, = struct.unpack("!Q", buf.getvalue())
        if maxsize is not None and size > maxsize:
            return None
        return self._recv(size)

    Connection._send_bytes = send_bytes
    Connection._recv_bytes = recv_bytes

    if log is not None:
        log.add(patchname + " applied")

def hist_profile(x, y, bins, range_x,range_y,outlier_insensitive=False):
    w = (y>range_y[0]) & (y<range_y[1])
    if(outlier_insensitive):
        means_result = binned_statistic(x[w], y[w], bins=bins, range=range_x, statistic='median')
        outlier_insensitive_std = lambda x : (np.nanpercentile(x,84.135,axis=0)-np.nanpercentile(x,15.865,axis=0))/2
        std_result = binned_statistic(x[w], y[w], bins=bins, range=range_x, statistic=outlier_insensitive_std)
        nb_entries_result = binned_statistic(x[w], y[w], bins=bins, range=range_x, statistic='count')

    else:
        means_result = binned_statistic(x[w], y[w], bins=bins, range=range_x, statistic='mean')
        std_result = binned_statistic(x[w], y[w], bins=bins, range=range_x, statistic='std')
        nb_entries_result = binned_statistic(x[w], y[w], bins=bins, range=range_x, statistic='count')

    means = means_result.statistic
    std = std_result.statistic
    nb_entries = nb_entries_result.statistic

    errors = std/np.sqrt(nb_entries)

    bin_edges = means_result.bin_edges
    bin_centers = (bin_edges[:-1] + bin_edges[1:])/2.
    return(bin_centers, means, errors)

#### PLOT ####

def plot_histo(value,value_name,name,dotted=False,**kwargs):
    nb_bins = return_key(kwargs,f"{value_name}_bins",50)
    value_min = return_key(kwargs,f"{value_name}_value_min",np.min(value))
    value_max = return_key(kwargs,f"{value_name}_value_max",np.max(value))
    alpha = return_key(kwargs,f"{value_name}_alpha",1.0)
    histtype=return_key(kwargs,f"{value_name}_histtype",'bar')
    linestyle=return_key(kwargs,f"{value_name}_linestyle",None)
    lw=return_key(kwargs,f"{value_name}_linewidth",0)
    ec=return_key(kwargs,f"{value_name}_edgecolor",None)
    color=return_key(kwargs,f"{value_name}_color",None)
    if(dotted):
        linestyle = "dashed"
        color = None
        facecolor="None"
        lw = 2
        ec = "k"

    norm = return_key(kwargs,f"{value_name}_norm",False)
    cumulative = return_key(kwargs,f"{value_name}_cumulative",False)
    log = return_key(kwargs,f"{value_name}_log",False)

    if(norm):
        name = name + "_normalized"
    if(cumulative):
        name = name + "_cumulative"

    if((value_min is None)|(value_min is None)):
        bins = nb_bins
    else:
        bins = np.linspace(value_min,value_max, nb_bins)
    if(dotted):
        (n, bins, patches) = plt.hist(value, bins, alpha=alpha,histtype=histtype,
                                      linestyle=linestyle,ec=ec,
                                      density=norm,cumulative=cumulative,facecolor=facecolor,
                                      log=log,color=color,lw=lw)
    else:
        (n, bins, patches) = plt.hist(value, bins, alpha=alpha,histtype=histtype,
                                      linestyle=linestyle,ec=ec,
                                      density=norm,cumulative=cumulative,
                                      log=log,color=color,lw=lw)
    return(name, n, bins, patches)

def save_histo(pwd,value,value_name,name,comparison=None,comparison_legend=None,**kwargs):
    xlabel = return_key(kwargs,f"{value_name}_xlabel",value_name)
    ylabel = return_key(kwargs,f"{value_name}_ylabel","#")
    min_lim = return_key(kwargs,f"{value_name}_min_lim",np.min(value))
    max_lim = return_key(kwargs,f"{value_name}_max_lim",np.max(value))
    fontsize = return_key(kwargs,f"{value_name}_fontsize",14)
    fontsize_scale = return_key(kwargs,f"{value_name}_fontscalesize",14)
    figsize = return_key(kwargs,f"{value_name}_figsize",(8,5))
    comparison_dotted = return_key(kwargs,f"{value_name}_comparison_dotted",None)
    plt.figure(figsize=figsize)
    name_out, n, bins, patches = plot_histo(value,value_name,name,**kwargs)
    if(comparison is not None):
        for i in range(len(comparison)):
            dotted=False
            if(comparison_dotted is not None):
                if(comparison_dotted == i):
                    dotted=True
            plot_histo(comparison[i],value_name,name,dotted=dotted,**kwargs)
        if(comparison_legend is not None):
            plt.legend(comparison_legend,fontsize=fontsize)
    plt.ylabel(ylabel,fontsize=fontsize)
    plt.xlabel(xlabel,fontsize=fontsize)
    ax = plt.gca()
    ax.tick_params(axis='x', labelsize=fontsize_scale)
    ax.tick_params(axis='y', labelsize=fontsize_scale)

    if((min_lim is not None)&(max_lim is not None)):
        plt.xlim([min_lim,max_lim])

    plt.tight_layout()
    plt.savefig(os.path.join(pwd,f"{name_out}_histo_{value_name}.pdf"), format ="pdf")

def plot_mean_redshift_dependence(value,redshift,value_name,name,**kwargs):
    ax = return_key(kwargs,f"{value_name}_ax",None)
    nb_bins = return_key(kwargs,f"{value_name}_z_bins",50)
    ls=return_key(kwargs,f"{value_name}_linestyle",None)
    color=return_key(kwargs,f"{value_name}_color",None)
    marker=return_key(kwargs,f"{value_name}_marker",'.')
    lambda_obs=return_key(kwargs,f"{value_name}_lambda_obs",False)
    lambda_rest=return_key(kwargs,f"{value_name}_lambda_rest",False)
    multiplicative_coef=return_key(kwargs,f"{value_name}_multiplicative_coef",None)

    outlier_insensitive = return_key(kwargs,f"{value_name}_outlier_insensitive",False)

    if(ax is None):
        ax = plt.gca()

    if(multiplicative_coef is not None):
        value = value * multiplicative_coef

    if(lambda_obs):
        redshift = (1 + redshift)*lambdaLy

    elif(lambda_rest):
        redshift_qso=return_key(kwargs,f"{value_name}_redshift_qso",None)
        if(redshift_qso is None):
            raise ValueError("""You have chosen to convert your mean redshift
                                dependence plot rest frame wavelength but no
                                QSO redshift was provided""")
        redshift = ((1 + redshift)/(1+ redshift_qso))*lambdaLy

    z_min = return_key(kwargs,f"{value_name}_zmin",np.min(redshift))
    z_max = return_key(kwargs,f"{value_name}_zmax",np.max(redshift))
    range_x = np.array([z_min,z_max])

    bin_centers, means, errors = hist_profile(redshift, value, nb_bins, range_x,
                                              [np.min(value),np.max(value)],
                                              outlier_insensitive=outlier_insensitive)

    ax.errorbar(bin_centers,means, errors,color=color,ls=ls, marker=marker)

    if(outlier_insensitive):
        name = name + "_outlier_insensitive"

    return(name)

def save_mean_redshift_dependence(pwd,value,redshift,value_name,name,
                                  comparison=None,
                                  comparison_redshift=None,
                                  comparison_legend=None,**kwargs):
    ylabel = return_key(kwargs,f"{value_name}_xlabel",value_name)
    lambda_obs=return_key(kwargs,f"{value_name}_lambda_obs",False)
    lambda_rest=return_key(kwargs,f"{value_name}_lambda_rest",False)
    if((lambda_obs)&(lambda_rest)):
        raise ValueError(f"""You have chosen to convert your mean redshift
                             dependence plot to both rest frame and observed
                             wavelength. Please choose {value_name}_lambda_obs
                             or {value_name}_lambda_rest""")
    if(lambda_obs):
        default_xlabel = "observed wavelength"
    elif(lambda_rest):
        default_xlabel = "rest frame wavelength"
    else:
        default_xlabel = "redshift"
    xlabel = return_key(kwargs,f"{value_name}_ylabel",default_xlabel)

    plt.figure()
    name_out = plot_mean_redshift_dependence(value,redshift,value_name,name,**kwargs)
    if(comparison is not None):
        for i in range(len(comparison)):
            plot_mean_redshift_dependence(comparison[i],comparison_redshift[i],value_name,name,**kwargs)
        if(comparison_legend is not None):
            plt.legend(comparison_legend)
    plt.ylabel(ylabel)
    plt.xlabel(xlabel)
    plt.savefig(os.path.join(pwd,f"{name_out}_mean_redshift_dependence_{value_name}.pdf"), format ="pdf")

def plot_redshift_dependence(value,redshift,value_name,name,**kwargs):
    z_min = return_key(kwargs,f"{value_name}_zmin",np.min(redshift))
    z_max = return_key(kwargs,f"{value_name}_zmax",np.max(redshift))

    mask = (redshift >= z_min)&(redshift < z_max)
    plt.scatter(redshift[mask],value[mask])
    return(name)

def save_redshift_dependence(pwd,value,redshift,value_name,name,
                             comparison=None,
                             comparison_redshift=None,
                             comparison_legend=None,**kwargs):
    ylabel = return_key(kwargs,f"{value_name}_xlabel",value_name)
    xlabel = return_key(kwargs,f"{value_name}_ylabel","redshift")

    plt.figure()
    name_out = plot_redshift_dependence(value,redshift,value_name,name,**kwargs)
    if(comparison is not None):
        for i in range(len(comparison)):
            plot_redshift_dependence(comparison[i],comparison_redshift[i],value_name,name,**kwargs)
        if(comparison_legend is not None):
            plt.legend(comparison_legend)
    plt.ylabel(ylabel)
    plt.xlabel(xlabel)
    plt.savefig(os.path.join(pwd,f"{name_out}_redshift_dependence_{value_name}.pdf"), format ="pdf")

def plot_ra_dec(ra,dec,name,**kwargs):

    nb_cut = return_key(kwargs,"nb_cut",None)
    figsize = return_key(kwargs,"ra_dec_figsize",(7,3.5))
    plt.figure(figsize=figsize)

    if(nb_cut is not None):
        ramax,ramin,decmax,decmin = np.max(ra),np.min(ra),np.max(dec),np.min(dec)
        interval_ra_array = []
        for cut in range(nb_cut):
            interval_ra_array.append([((cut)/(nb_cut))*(ramax-ramin) + ramin  , ((cut + 1)/(nb_cut))*(ramax-ramin) + ramin])
        for i in range(len(interval_ra_array)):
            plt.plot([interval_ra_array[i][0],interval_ra_array[i][1]],[decmin,decmin],color="orange", linewidth=2)
            plt.plot([interval_ra_array[i][0],interval_ra_array[i][1]],[decmax,decmax],color="orange", linewidth=2)
            plt.plot([interval_ra_array[i][0],interval_ra_array[i][0]],[decmin,decmax],color="orange", linewidth=2)
            plt.plot([interval_ra_array[i][1],interval_ra_array[i][1]],[decmin,decmax],color="orange", linewidth=2)

    plt.plot(ra,dec,'b.', markersize=1.5)
    return(name)

def save_ra_dec(pwd,ra,dec,name,
                comparison_ra=None,comparison_dec=None,
                comparison_legend=None,**kwargs):

    grid = return_key(kwargs,"ra_dec_grid",True)
    fontsize = return_key(kwargs,"ra_dec_fontsize",13)
    fontsize_scale = return_key(kwargs,"ra_dec_fontscalesize",13)
    ra_min_lim = return_key(kwargs,"ra_dec_ra_min_lim",np.min(ra))
    ra_max_lim = return_key(kwargs,"ra_dec_ra_max_lim",np.max(ra))
    dec_min_lim = return_key(kwargs,"ra_dec_dec_min_lim",np.min(dec))
    dec_max_lim = return_key(kwargs,"ra_dec_dec_max_lim",np.max(dec))
    deg = return_key(kwargs,"deg",True)
    if(deg):
        label_angle = "deg"
    else:
        label_angle = "rad"
    ylabel = return_key(kwargs,"ra_dec_ylabel",f"DEC [{label_angle}] (J2000)")
    xlabel = return_key(kwargs,"ra_dec_xlabel",f"RA [{label_angle}] (J2000)")

    name_out = plot_ra_dec(ra,dec,name,**kwargs)
    if(comparison_ra is not None):
        for i in range(len(comparison_ra)):
            plot_redshift_dependence(comparison_ra[i],comparison_dec[i],name,**kwargs)
        if(comparison_legend is not None):
            plt.legend(comparison_legend)
    plt.xlabel(xlabel, fontsize=fontsize)
    plt.ylabel(ylabel, fontsize=fontsize)
    ax = plt.gca()
    ax.tick_params(axis='x', labelsize=fontsize_scale)
    ax.tick_params(axis='y', labelsize=fontsize_scale)
    plt.xlim([ra_min_lim,ra_max_lim])
    plt.ylim([dec_min_lim,dec_max_lim])
    if(grid):
        plt.grid()
    plt.tight_layout()
    plt.savefig(os.path.join(pwd,f"{name_out}_RA-DEC_diagram.pdf"), format ="pdf")

max_stars_repo_name: corentinravoux/lelantos

max_issues_repo_name: corentinravoux/lelantos

max_forks_repo_name: corentinravoux/lelantos

reason: import numpy,from scipy


text:
# coding: utf-8

"""
    description: NLP utils
    author: Suraj Iyer
"""

__all__ = [
    'words',
    'count_words',
    'edits_1',
    'edits_n',
    'cluster_words_by_edit_distance',
    'join_bigrams_with_replacements',
    'correct_word_compounding',
    # 'tf_idf',
    'knn_name_matching',
    'bigram_context',
    'cluster_urls',
    'corpus_level_tfidf',
    'extract_phrases',
    'remove_emojis'
]

import numpy as np
import pandas as pd
import re
from collections import Counter, defaultdict
import nltk
from python_data_utils.utils import load_artifact
import difflib
from typing import Iterable, Tuple

def words(text: str) -> str:
    return re.findall(r'\w+', text.lower())

def count_words(
        corpus: str = None, file_path: str = None) -> Counter:
    if corpus:
        return Counter(words(corpus))
    elif file_path:
        cnt = Counter()
        with open(file_path, 'r', encoding='utf-8') as f:
            for line in f:
                for word in line.split():
                    cnt[word] += 1
        return cnt

def edits_1(word: str) -> set:
    """All edits that are one edit away from `word`."""
    assert isinstance(word, str), '{} must be a string'.format(word)
    letters = 'abcdefghijklmnopqrstuvwxyz'
    splits = [(word[:i], word[i:]) for i in range(len(word) + 1)]
    deletes = [L + R[1:] for L, R in splits if R]
    transposes = [L + R[1] + R[0] + R[2:] for L, R in splits if len(R) > 1]
    replaces = [L + c + R[1:] for L, R in splits if R for c in letters]
    inserts = [L + c + R for L, R in splits for c in letters]
    return set(deletes + transposes + replaces + inserts)

def edits_n(word: str, dist: int = 2, k: int = None) -> tuple:
    """All edits that are `n` edits away from `word`."""
    if k is None:
        k = dist
    assert dist > 0 and k > 0
    if k == 1:
        return edits_1(word)
    new_words = tuple(
        w2 for w1 in edits_1(word)
        for w2 in edits_n(w1, dist, k - 1))
    return tuple(set(new_words)) if k == dist else new_words

def cluster_words_by_edit_distance(
        words: Iterable[str], dist: int = 2) -> dict:
    assert isinstance(words, (list, tuple))
    assert any(isinstance(w, str) for w in words), 'words must be an iterable of strings'
    assert isinstance(dist, int)
    cluster = dict()
    import distance
    for i, w in enumerate(words):
        x = [words[j] for j in range(
            i + 1, len(words)) if distance.levenshtein(w, words[j]) == dist]
        if len(x) > 0:
            cluster[w] = x
    return cluster

def join_bigrams_with_replacements(
        s: Iterable[Iterable['str']],
        replacements=lambda w, default: default) -> str:
    """
    Function to join a list of bigrams back into a sentence taking into account
    unigram replacements for bigrams.

    working version: 1.0
    # does not work with generator input
    skip = False
    string = ""
    n = len(s)
    for i, w in enumerate(s):
        if skip:
            skip = False
            if i == n-1:
                string = string + " " + w[1]
            continue
        w = replacements(" ".join(w), w)
        if isinstance(w, str):
            string = string + " " + w
            skip = True
        else:
            string = string + " " + w[0]
    return string

    param:
        s: iterable of list of strings
            Example: [("I", "have"), ("have", "to"), ("to", "go")]
        replacements:
            Function that takes the following parameters. By default,
            it returns the default value which simply joins all the bigrams
            together.
                w: string
                    A bigram string. Example: "wi fi"
                default:
                    A default value to return in case no replacement available.
                returns: string
                    A replacement string, e.g., "wifi"
        returns: string
            A string sentence joining all bigrams correctly including the
            replacements (if any).
    """
    s = iter(s)
    first = next(s, None)
    if first is None:
        return ""
    first = replacements(" ".join(first), first)
    string = first if isinstance(first, str) else first[0] + " " + first[1]
    prev_string = "" if isinstance(first, str) else first[0]
    for w in s:
        w = replacements(" ".join(w), w)
        if isinstance(w, str):
            string = prev_string + " " + w
            prev_string = string
        else:
            prev_string = string
            string += " " + w[1]
    return string

def correct_word_compounding(
        dfs: pd.Series, language_dictionary,
        split_on_both_accurate: bool = True,
        threshold: int = 3) -> pd.DataFrame:
    """
    List all bigrams which are also written as unigrams
    by removing one space character
    and filter the ones which contains words not in any dictionary

    Assumptions: every word is unique is only written in
    one way and there are no spelling errors

    Heuristic:
        if bigram counts ≈ unibigram counts:
            # Assumption: both versions are accurate enough
            if not unibigram in language dictionary:
                # Assumption: only bigram version accurate
                split unibigrams to unigrams (make bigram)
            else:
                # Assumption: both versions are still accurate enough
                if both unigrams in language dictionary:
                    # Assumption: both versions are still accurate enough
                    # so depend on {:split_on_both_accurate:}
                    if split_on_both_accurate:
                        make bigram
                    else:
                        convert bigram to unigram by removing space
                        (make unibigram)
                else:
                    # Assumption: only unibigram version accurate
                    make unibigram
        else:
            # Assumption: Only one version is accurate
            if unibigram count >> bigram count:
                # Assumption: Unibigram version is more accurate
                if unibigram in language dictionary:
                    # Assumption: Unibigram version is still more accurate
                    make unibigram
                else:
                    # Assumption: Unibigram version is not accurate so
                    # split instead
                    if both unigrams in language dictionary:
                        make bigram
                    else:
                        make unibigram
            else:
                # Assumption: Bigram version is more accurate
                if both unigrams in language dictionary:
                        make bigram
                    else:
                        make unibigram
    """
    assert isinstance(threshold, int) and threshold > 0
    assert isinstance(split_on_both_accurate, bool)
    tokens = dfs.str.split(' ')
    unigrams = set(w for s in tokens for w in s)
    bigrams = set(" ".join(w) for s in tokens.apply(nltk.bigrams) for w in s)
    bigrams_also_in_unigrams = set(
        b for b in bigrams if b.replace(" ", "") in unigrams)

    bigrams = pd.DataFrame([(
        b,
        b.replace(" ", ""),
        b.replace(" ", "-"),
        dfs.str.contains(b).sum(),
        dfs.str.contains(b.replace(" ", "") + "|" + b.replace(" ", "-")).sum()
    ) for b in bigrams_also_in_unigrams], columns=[
        'bigram', 'unibigram0', 'unibigram1',
        'count_occurence_bigram', 'count_occurence_unibigram'])
    bigrams['near_equal'] = abs(
        bigrams["count_occurence_bigram"] -
        bigrams["count_occurence_unibigram"]) < threshold
    bigrams['unibigram_in_dict?'] = bigrams['bigram'].apply(
        lambda bi: bi.replace(" ", "") in language_dictionary)
    bigrams['unigram0_in_dict?'] = bigrams['bigram'].apply(
        lambda bi: bi.split(' ')[0] in language_dictionary)
    bigrams['unigram1_in_dict?'] = bigrams['bigram'].apply(
        lambda bi: bi.split(' ')[1] in language_dictionary)
    bigrams['both_unigrams_in_dict?'] = (
        bigrams['unigram0_in_dict?']) & (bigrams['unigram1_in_dict?'])
    make_bigram_mask = (
        (
            (bigrams["near_equal"]) &
            (
                (~bigrams["unibigram_in_dict?"]) |
                (bigrams["unibigram_in_dict?"] &
                    bigrams["both_unigrams_in_dict?"] &
                    split_on_both_accurate)
            )
        ) |
        (
            (~bigrams["near_equal"]) &
            (
                (
                    (bigrams["count_occurence_unibigram"] >
                        bigrams["count_occurence_bigram"]) &
                    (~bigrams["unibigram_in_dict?"]) &
                    (bigrams['both_unigrams_in_dict?'])
                ) |
                (
                    (bigrams["count_occurence_unibigram"] <=
                        bigrams["count_occurence_bigram"]) &
                    (bigrams['both_unigrams_in_dict?'])
                )
            )
        )
    )

    # List of word replacements to perform on the text
    replacements = {}
    replacements.update({
        w[1][f"unibigram{i}"]: w[1]["bigram"]
        for w in bigrams[make_bigram_mask].iterrows() for i in range(2)})
    replacements.update({
        w[1]["bigram"]: w[1]["unibigram0"]
        for w in bigrams[~make_bigram_mask].iterrows()})

    # Make replacements
    tokens = tokens.apply(lambda s: [replacements.get(w, w) for w in s])
    tokens.loc[tokens.apply(len) > 1] = tokens.loc[tokens.apply(len) > 1]\
        .apply(nltk.bigrams)\
        .apply(join_bigrams_with_replacements, replacements=replacements.get)
    tokens.loc[tokens.apply(len) == 1] = tokens.loc[tokens.apply(len) == 1]\
        .apply(lambda s: "".join(s))
    return tokens

# def tf_idf(documents: Iterable[str]) -> Tuple[dict, dict, dict]:
#     N = len(documents)
#     assert N > 0, "Count of documents must be at least 1."
#     corpus = " ".join(documents)
#     corpus_deduplicated = " ".join(
#         [" ".join(frozenset(d.split(" "))) for d in documents])

#     tf = count_words(corpus=corpus)
#     df = count_words(corpus=corpus_deduplicated)
#     tf_idf = {t: tf * np.log(1 + N / df[t]) for t, tf in tf.items()}

#     return tf_idf, tf, df

def knn_name_matching(
        A: Iterable[str], B: Iterable[str],
        vectorizer_kws: dict = {}, nn_kws: dict = {},
        max_distance: float = None, return_B=True) -> list:
    """
    Nearest neighbor name matching of sentences in B to A.
    """
    from sklearn.neighbors import NearestNeighbors
    from sklearn.feature_extraction.text import TfidfVectorizer

    # vectorize the B documents after fitting on A
    vectorizer = TfidfVectorizer(**vectorizer_kws)
    Xa = vectorizer.fit_transform(A)
    Xb = vectorizer.transform(B)

    # find nearest neighbor matching
    neigh = NearestNeighbors(n_neighbors=1, **nn_kws)
    neigh.fit(Xa)
    if max_distance is None:
        indices = neigh.kneighbors(Xb, return_distance=False).flatten()
    else:
        indices, distances = neigh.kneighbors(Xb)
        indices, distances = indices.flatten(), distances.flatten()
        indices = indices[distances <= max_distance]

    if return_B:
        result = [(B[i], A[idx]) for i, idx in enumerate(indices)]
    else:
        result = [A[idx] for idx in indices]

    return result

def bigram_context(
        documents: Iterable[str], window_size: int = 2,
        unique: bool = False, sort: bool = True) -> pd.DataFrame:
    """
    Get bigram context (word pairs) within given window size
    from given documents
    """
    words = []
    for d in documents:
        words.extend(d)
        words.extend(["|"] * (window_size - 1))
    finder = nltk.collocations.BigramCollocationFinder\
        .from_words(words, window_size=window_size)
    finder.apply_ngram_filter(lambda *w: "|" in w)
    df = pd.DataFrame(finder.ngram_fd.items(), columns=["pair", "count"])
    df[["left", "right"]] = df['pair'].apply(pd.Series)
    df.drop("pair", axis=1, inplace=True)
    if unique:
        df.loc[df.left >= df.right, ['left', 'right']] =\
            df.loc[df.left >= df.right, ['right', 'left']].values
        df = df.groupby(['left', 'right'])['count'].sum().reset_index()
    if sort:
        df.sort_values("count", ascending=False, inplace=True)
        df.reset_index(drop=True, inplace=True)
    return df

def cluster_urls(
        urls: Iterable[str],
        min_cluster_size: int = 10) -> pd.DataFrame:
    """
    Cluster URLs by regex rules defined in this package:
    https://pypi.org/project/urlclustering/

    urls: list
        List of urls.
    min_clustre_size: int
        Minimum cluster size
    """
    import urlclustering
    clusters = urlclustering.cluster(urls, min_cluster_size)
    tmp = {v0: [k[1], 0] for k, v in clusters['clusters'].items() for v0 in v}
    tmp.update({k: [k, 1] for k in clusters['unclustered']})
    clusters = pd.DataFrame.from_dict(
        tmp, orient='index', columns=['cluster', 'unclustered'])
    return clusters

def corpus_level_tfidf(
        texts: Iterable[str], **vectorizer_kwargs) -> Tuple[Iterable[float], Iterable[str]]:
    """
    See book: Ziegler, C.N., 2012. Mining for strategic competitive
    intelligence. Springer. pp. 128-130.
    """
    from sklearn.feature_extraction.text import CountVectorizer

    with np.errstate(divide='ignore'):
        # Compute 1 + log(tf+(t_i)) where tf+ computes the summed term
        # frequency for term t_i and t_i occurs in some text `d` \in texts.
        vectorizer_kwargs.pop('min_df', None)
        vectorizer_kwargs.pop('max_df', None)
        tf_vect = CountVectorizer(**vectorizer_kwargs)
        tf = np.array(tf_vect.fit_transform(texts).sum(0)).reshape(-1)
        tf = (tf > 0).astype(float) + np.log(tf, where=tf > 0)

        # Get the order of the terms in the tf_vect
        terms = tf_vect.get_feature_names()

        # Compute log(|texts| / df(t_i)) where df computes the document
        # frequency of t_i in texts.
        vectorizer_kwargs.pop('vocabulary', None)
        vectorizer_kwargs.pop('binary', None)
        idf_vect = CountVectorizer(
            binary=True, vocabulary=terms, **vectorizer_kwargs)
        idf = np.array(np.log(
            len(texts) / (1. + idf_vect.fit_transform(texts).sum(0))
        )).reshape(-1)
        return tf * idf, terms

def extract_phrases(
        list_of_strings: Iterable[str], min_phrase_length: int = 2,
        max_phrase_length: int = float("inf"), min_freq: int = 1,
        delimiter: str = "||") -> dict:
    N = len(list_of_strings)
    assert N > 0
    phrases = defaultdict(int)
    string = list_of_strings[0]
    if N == 1:
        phrases[string] = 1
    else:
        for i in range(1, len(list_of_strings)):
            matches = difflib\
                .SequenceMatcher(None, list_of_strings[i], string)\
                .get_matching_blocks()[:-1]
            string += delimiter
            end = 0
            for match in matches:
                string += " " + list_of_strings[i][end:match.a]
                end = match.a + match.size
                phrases[list_of_strings[i][match.a:end]] += 1
            string += " " + list_of_strings[i][end:]
    phrases = {
        k: v for k, v in map(lambda i: (i[0].strip(), i[1]), phrases.items())
        if (v >= min_freq) and
        (min_phrase_length <= len(k) <= max_phrase_length)}
    return phrases

def remove_emojis(inputString):
    # https://stackoverflow.com/questions/33404752/removing-emojis-from-a-string-in-python
    return inputString.encode('ascii', 'ignore').decode('ascii')

ext: py

max_stars_repo_name: surajiyer/python-data-utils

max_issues_repo_name: surajiyer/python-data-utils

max_forks_repo_name: surajiyer/python-data-utils

reason: import numpy


text:
# coding: utf-8

# # Nengo Example: Two Neurons
#
# This demo shows how to construct and manipulate a complementary pair of neurons.
#
# These are leaky integrate-and-fire (LIF) neurons. The neuron tuning properties have been selected so there is one ‘on’ and one ‘off’ neuron.
#
# One neuron will increase for positive input, and the other will decrease. This can be thought of as the simplest population that is able to give a reasonable representation of a scalar value.

# In[ ]:

import numpy as np
import matplotlib.pyplot as plt

import nengo

# ## Step 1: Create the neurons

# In[ ]:

from nengo.dists import Uniform

model = nengo.Network(label='Two Neurons')
with model:
    neurons = nengo.Ensemble(2, dimensions=1,  # Representing a scalar
                             intercepts=Uniform(-.5, -.5),  # Set the intercepts at .5
                             max_rates=Uniform(100,100),  # Set the max firing rate at 100hz
                             encoders=[[1],[-1]])  # One 'on' and one 'off' neuron

# ## Step 2: Create input for the model
#
# Create an input node generating a sine wave.

# In[ ]:

with model:
    sin = nengo.Node(lambda t: np.sin(8 * t))

# ## Step 3: Connect the network elements

# In[ ]:

with model:
    nengo.Connection(sin, neurons, synapse=0.01)

# ## Step 4: Probe outputs
#
# Anything that is probed will collect the data it produces over time, allowing us to analyze and visualize it later.

# In[ ]:

with model:
    sin_probe = nengo.Probe(sin)  # The original input
    spikes = nengo.Probe(neurons.neurons)  # Raw spikes from each neuron
    voltage = nengo.Probe(neurons.neurons, 'voltage')  # Subthreshold soma voltages of the neurons
    filtered = nengo.Probe(neurons, synapse=0.01)  # Spikes filtered by a 10ms post-synaptic filter

# ## Step 5: Run the model

# In[ ]:

with nengo.Simulator(model) as sim:  # Create a simulator
    sim.run(1)  # Run it for 1 second

# ## Step 6: Plot the results

# In[ ]:

# Plot the decoded output of the ensemble
plt.plot(sim.trange(), sim.data[filtered])
plt.plot(sim.trange(), sim.data[sin_probe])
plt.xlim(0, 1)

# Plot the spiking output of the ensemble
from nengo.utils.matplotlib import rasterplot
plt.figure(figsize=(10, 8))
plt.subplot(221)
rasterplot(sim.trange(), sim.data[spikes], colors=[(1, 0, 0), (0, 0, 0)])
plt.xlim(0, 1)
plt.yticks((0, 1), ("On neuron", "Off neuron"))

# Plot the soma voltages of the neurons
plt.subplot(222)
plt.plot(sim.trange(), sim.data[voltage][:, 0] + 1, 'r')
plt.plot(sim.trange(), sim.data[voltage][:, 1], 'k')
plt.yticks(())
plt.axis([0, 1, 0, 2])
plt.subplots_adjust(wspace=0.05);

# The top graph shows that the input signal in green and the filtered output spikes from the two neurons population in blue.  The spikes (that are filtered) from the 'on' and 'off' neurons are shown in the bottom graph on the left.  On the right are the subthreshold voltages for the neurons.

max_stars_repo_name: tbekolay/nengodocs-rtd

max_issues_repo_name: tbekolay/nengodocs-rtd

max_forks_repo_name: tbekolay/nengodocs-rtd


text:
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF

def generate_regression_task_dataset(n_tasks):
    kernel = RBF(length_scale=1., length_scale_bounds=(1e-5, 1e5))
    gpr = GaussianProcessRegressor(kernel=kernel)
    X = np.linspace(0, 100, 1000)
    X = np.atleast_2d(X).T
    f_samples = []
    for i in range(n_tasks):
        # sample a function from a zero mean GP
        f = gpr.sample_y(X, random_state=i)
        f_samples.append(f)
    return f_samples


text:
import xarray as xr
import numpy as np

def build_rot_mat(roll: xr.DataArray, pitch: xr.DataArray, yaw: xr.DataArray, order: str = 'rpy', degrees: bool = True):
    """
    Make the rotation matrix for a set of angles and return the matrix.
    All file angles are in degrees, so incoming angles are degrees.

    Intrinsic - each rotation performed on coordinate system as rotated by previous operation
    Intrinsic rotation, rot(rpy) = rot(y)*rot(p)*rot(r)

    Parameters
    ----------
    roll
        array of floating point numbers from roll sensor
    pitch
        array of floating point numbers from pitch sensor
    yaw
        array of floating point numbers from heading sensor
    order
        order of rotation, either 'rpy' or 'ypr'
    degrees
        True if incoming angles are in degrees, False if radians

    Returns
    -------
    xr.DataArray
        rotation matrix composed of rpy rotations
    """

    if type(roll) != xr.DataArray or type(pitch) != xr.DataArray or type(yaw) != xr.DataArray:
        raise TypeError('Expected xarray DataArray object')

    if order == 'ypr':
        r = yaw
        p = pitch
        y = roll
    elif order == 'rpy':
        r = roll
        p = pitch
        y = yaw
    else:
        raise ValueError('Order provided is not rpy or ypr.')

    if degrees:
        r = np.deg2rad(r)
        p = np.deg2rad(p)
        y = np.deg2rad(y)

    rcos = np.cos(r)
    pcos = np.cos(p)
    ycos = np.cos(y)
    rsin = np.sin(r)
    psin = np.sin(p)
    ysin = np.sin(y)

    r00 = (ycos * pcos).assign_coords({'x': 0, 'y': 0})
    r01 = (ycos * psin * rsin - ysin * rcos).assign_coords({'x': 0, 'y': 1})
    r02 = (ycos * psin * rcos + ysin * rsin).assign_coords({'x': 0, 'y': 2})
    r0 = xr.concat([r00, r01, r02], dim='y')
    r10 = (ysin * pcos).assign_coords({'x': 1, 'y': 0})
    r11 = (ysin * psin * rsin + ycos * rcos).assign_coords({'x': 1, 'y': 1})
    r12 = (ysin * psin * rcos - ycos * rsin).assign_coords({'x': 1, 'y': 2})
    r1 = xr.concat([r10, r11, r12], dim='y')
    r20 = (-psin).assign_coords({'x': 2, 'y': 0})
    r21 = (pcos * rsin).assign_coords({'x': 2, 'y': 1})
    r22 = (pcos * rcos).assign_coords({'x': 2, 'y': 2})
    r2 = xr.concat([r20, r21, r22], dim='y')

    rmat = xr.concat([r0, r1, r2], dim='x').transpose('time', 'x', 'y')

    return rmat

def build_mounting_angle_mat(roll: float, pitch: float, yaw: float, tstmp: str):
    """
    Feeds build_rot_mat, difference being this takes in single floating point numbers for rpy as you get from a
    surveyed mount angle data point.

    Assumes angles are in degrees and rpy rotation is desired.

    Parameters
    ----------
    roll
        roll angle for rotation matrix
    pitch
        pitch angle for rotation matrix
    yaw
        yaw angle for rotation matrix
    tstmp
        time relevant installation parameter showed up in the multibeam file

    Returns
    -------
    xr.DataArray
        rotation matrix composed of data provided
    """

    if type(roll) != float or type(pitch) != float or type(yaw) != float:
        raise TypeError('Expected floating point values for roll,pitch,yaw')

    time_coord = np.array([float(tstmp)])
    roll_xarr = xr.DataArray(np.array([roll], dtype=np.float32), dims=['time'], coords={'time': time_coord}).chunk()
    pitch_xarr = xr.DataArray(np.array([pitch], dtype=np.float32), dims=['time'], coords={'time': time_coord}).chunk()
    yaw_xarr = xr.DataArray(np.array([yaw], dtype=np.float32), dims=['time'], coords={'time': time_coord}).chunk()
    return build_rot_mat(roll_xarr, pitch_xarr, yaw_xarr, order='rpy', degrees=True)

def combine_rotation_matrix(mat_one: xr.DataArray, mat_two: xr.DataArray):
    """
    Composing two rotation matrices is performed by taking the product of the two matrices

    Assumes one of the input matrices is of size one (the mounting angle matrix, attitude changes over time)

    Order is important here.

    Parameters
    ----------
    mat_one
        3dim rotation matrix (time, x, y)
    mat_two
        3dim rotation matrix (time, x, y)

    Returns
    -------
    xr.DataArray
        3dim rotation matrix (time, x, y) for each time in input matrices
    """

    # This is apparently close, but not the right expression.  I can't figure this out right now
    # final_rot = einsum('ijk,jkl->ijk', mat_one, mat_two)

    # find the one element matrix
    if mat_one.shape[0] == 1:
        mat_one = mat_one.values
    elif mat_two.shape[0] == 1:
        mat_two = mat_two.values
    else:
        raise NotImplementedError('One of the input matrices must only have one value in the time dimension')

    # we'll just brute force it for now
    r00 = (mat_one[:, 0, 0] * mat_two[:, 0, 0]) + (mat_one[:, 1, 0] * mat_two[:, 0, 1]) + (mat_one[:, 2, 0] * mat_two[:, 0, 2])
    r00['y'] = 0
    r01 = (mat_one[:, 0, 1] * mat_two[:, 0, 0]) + (mat_one[:, 1, 1] * mat_two[:, 0, 1]) + (mat_one[:, 2, 1] * mat_two[:, 0, 2])
    r01['y'] = 1
    r02 = (mat_one[:, 0, 2] * mat_two[:, 0, 0]) + (mat_one[:, 1, 2] * mat_two[:, 0, 1]) + (mat_one[:, 2, 2] * mat_two[:, 0, 2])
    r02['y'] = 2
    r0 = xr.concat([r00, r01, r02], dim='y')
    r0['x'] = 0

    r10 = (mat_one[:, 0, 0] * mat_two[:, 1, 0]) + (mat_one[:, 1, 0] * mat_two[:, 1, 1]) + (mat_one[:, 2, 0] * mat_two[:, 1, 2])
    r10['y'] = 0
    r11 = (mat_one[:, 0, 1] * mat_two[:, 1, 0]) + (mat_one[:, 1, 1] * mat_two[:, 1, 1]) + (mat_one[:, 2, 1] * mat_two[:, 1, 2])
    r11['y'] = 1
    r12 = (mat_one[:, 0, 2] * mat_two[:, 1, 0]) + (mat_one[:, 1, 2] * mat_two[:, 1, 1]) + (mat_one[:, 2, 2] * mat_two[:, 1, 2])
    r12['y'] = 2
    r1 = xr.concat([r10, r11, r12], dim='y')
    r1['x'] = 1

    r20 = (mat_one[:, 0, 0] * mat_two[:, 2, 0]) + (mat_one[:, 1, 0] * mat_two[:, 2, 1]) + (mat_one[:, 2, 0] * mat_two[:, 2, 2])
    r20['y'] = 0
    r21 = (mat_one[:, 0, 1] * mat_two[:, 2, 0]) + (mat_one[:, 1, 1] * mat_two[:, 2, 1]) + (mat_one[:, 2, 1] * mat_two[:, 2, 2])
    r21['y'] = 1
    r22 = (mat_one[:, 0, 2] * mat_two[:, 2, 0]) + (mat_one[:, 1, 2] * mat_two[:, 2, 1]) + (mat_one[:, 2, 2] * mat_two[:, 2, 2])
    r22['y'] = 2
    r2 = xr.concat([r20, r21, r22], dim='y')
    r2['x'] = 2

    rmat = xr.concat([r0, r1, r2], dim='x').transpose('time', 'x', 'y')
    return rmat

def return_attitude_rotation_matrix(attitude: xr.Dataset, time_index: np.array = None):
    """
    We start off doing calculations in array relative reference frame.  To get to a geographic reference frame, we
    need rotation matrices for attitude and mounting angles.  Here we construct the attitude rotation matrix, at
    specific times (time of ping, time of receive, etc.).  We also allow for additional selection to accommodate
    for when we get duplicate times due to pingtime + twtt sometimes getting the same time for some beams.

    Parameters
    ----------
    attitude
        1 dimensional dataset representing the full attitude record
    time_index
        optional, if provided, is a 1 dimensional array of integers representing the index of values you want to return
        from the attitude Dataset

    Returns
    -------
    xr.DataArray
        1 dimensional array representing attitude times
    xr.DataArray
        3 dims (time, x, y) containing rot matrix at each time

    """
    # generate rotation matrices for the transmit array at time of ping
    if time_index is not None:
        attitude = attitude.isel(time=time_index)
    att_times = attitude.time
    attitude_rotation = build_rot_mat(attitude['roll'], attitude['pitch'], attitude['heading'], order='rpy', degrees=True)
    return att_times, attitude_rotation

def return_mounting_rotation_matrix(roll: float, pitch: float, yaw: float, tstmp: str):
    """
    Using the xyzrph record from xarray_conversion generated rangeangle DataSet, build the static mounting angle
    rotation matrix, pulling from the appropriate timestamp

    Parameters
    ----------
    roll
        roll sonar mounting angle
    pitch
        pitch sonar mounting angle
    yaw
        yaw sonar mounting angle
    tstmp
        timestamp associated with this installation parameters entry

    Returns
    -------
    xr.DataArray
        3dim (time, x, y) rotation matrix composed of rpy rotations at time of xyzrph record (length of time dim is
        always 1)

    """
    mount_rotation = build_mounting_angle_mat(float(roll), float(pitch), float(yaw), tstmp)
    return mount_rotation

max_stars_repo_name: giumas/kluster

max_issues_repo_name: giumas/kluster

max_forks_repo_name: giumas/kluster


text:
import torch
from torch.autograd import Variable
import utils
import mydataset
from PIL import Image
import numpy as np
import crnn as crnn
import cv2
import torch.nn.functional as F
import keys
import config
gpu = True
if not torch.cuda.is_available():
    gpu = False

model_path = './crnn_models/CRNN-0618-10w_21_990.pth'
alphabet = keys.alphabet
print(len(alphabet))
imgH = config.imgH
imgW = config.imgW
model = crnn.CRNN(imgH, 1, len(alphabet) + 1, 256)
if gpu:
    model = model.cuda()
print('loading pretrained model from %s' % model_path)
if gpu:
    model.load_state_dict( torch.load( model_path ) )
else:
    model.load_state_dict(torch.load(model_path,map_location=lambda storage,loc:storage))

converter = utils.strLabelConverter(alphabet)
transformer = mydataset.resizeNormalize((imgW, imgH),is_test=True)

def recognize_downline(img,crnn_model=model):
    img = cv2.cvtColor( img, cv2.COLOR_BGR2RGB )
    image = Image.fromarray(np.uint8(img)).convert('L')
    image = transformer( image )
    if gpu:
        image = image.cuda()
    image = image.view( 1, *image.size() )
    image = Variable( image )

    model.eval()
    preds = model( image )

    preds = F.log_softmax(preds,2)
    conf, preds = preds.max( 2 )
    preds = preds.transpose( 1, 0 ).contiguous().view( -1 )

    preds_size = Variable( torch.IntTensor( [preds.size( 0 )] ) )
    raw_pred = converter.decode( preds.data, preds_size.data, raw=True )
    sim_pred = converter.decode( preds.data, preds_size.data, raw=False )
    return sim_pred.upper()

if __name__ == '__main__':
    import shutil
    saved_path = 'test_imgs/'
    wrong_results = list()
    with open('data_set/infofile_test.txt') as f:
        content = f.readlines()
        num_all = 0
        num_correct = 0
        for line in content:
            fname, label = line.split('g:')
            fname += 'g'
            label = label.replace('\r', '').replace('\n', '')
            img = cv2.imread(fname)
            res = recognize_downline(img)
            if res==label:
                num_correct+=1
            else:
                # new_name = saved_path + fname.split('/')[-1]
                # shutil.copyfile(fname, new_name)
                wrong_results.append('res:{} / label:{}'.format(res,label))
            num_all+=1
            print(fname,res==label,res,label)

        print(num_correct/num_all)
        # print(wrong_results)

max_stars_repo_name: edgardeng/machine-learning-pytorch

max_issues_repo_name: edgardeng/machine-learning-pytorch

max_forks_repo_name: edgardeng/machine-learning-pytorch


text:
# Title: RADDet
# Authors: Ao Zhang, Erlik Nowruzi, Robert Laganiere
import json
import os
import cv2
import numpy as np
import pickle
import random
import colorsys
from glob import glob

def readConfig(config_file_name = "./config.json"):
    """ Read the configure file (json). """
    with open(config_file_name) as json_file:
        config = json.load(json_file)
    return config

def readAnchorBoxes(anchor_boxes_file = "./anchors.txt"):
    """ Read the anchor boxes found by k means """
    anchor_boxes = []
    with open(anchor_boxes_file) as txt_file:
        lines = txt_file.readlines()
        for line in lines:
            box = []
            characters = line.split(" ")
            for charact in characters:
                if charact != '\n':
                    box.append(int(charact))
            anchor_boxes.append(box)
    if len(anchor_boxes) == 0:
        anchor_boxes = None
    else:
        anchor_boxes = np.array(anchor_boxes)
        if anchor_boxes_file is not "./anchors.txt":
            anchor_boxes = sortAnchorBoxes(anchor_boxes, "2D")
        else:
            anchor_boxes = sortAnchorBoxes(anchor_boxes)
    return anchor_boxes

def sortAnchorBoxes(anchor_boxes, mode="3D"):
    """ Sort anchor boxes according to its area """
    if mode == "3D":
        anchor_box_areas = anchor_boxes[:, 0] * anchor_boxes[:, 1] * anchor_boxes[:, 2]
    else:
        anchor_box_areas = anchor_boxes[:, 0] * anchor_boxes[:, 1]
    anchor_boxes_order = np.argsort(anchor_box_areas)
    # print(anchor_box_areas[anchor_boxes_order]) # print out the boxes' areas
    return anchor_boxes[anchor_boxes_order]

def getSequenceNumbers(radar_dir, data_format):
    """ Get all the numbers from input file names """
    assert isinstance(data_format, list)
    assert all(isinstance(x, str) for x in data_format)
    assert (len(data_format) == 2 or len(data_format) == 3)
    sequence_numbers = []
    for this_file in glob(os.path.join(radar_dir, "*"+data_format[-1])):
        digits = int("".join((s for s in this_file.replace(radar_dir, "") \
                    if s.isdigit())))
        sequence_numbers.append(digits)
    if len(sequence_numbers) == 0:
        sequence_numbers = None
    else:
        sequence_numbers = np.sort(sequence_numbers)
    return sequence_numbers

def randomColors(classes, bright=True):
    """ Define colors for all categories. """
    assert isinstance(classes, list)
    assert all(isinstance(x, str) for x in classes)
    N = len(classes)
    brightness = 1.0 if bright else 0.7
    hsv = [(i / N, 1, brightness) for i in range(N)]
    colors = list(map(lambda c: colorsys.hsv_to_rgb(*c), hsv))
    random.seed(8888)
    random.shuffle(colors)
    return colors

def gtfileFromRADfile(RAD_file, prefix):
    """ Transfer RAD filename to gt filename """
    RAD_file_spec = RAD_file.split("RAD")[-1]
    gt_file = os.path.join(prefix, "gt") + RAD_file_spec.replace("npy", "pickle")
    return gt_file

def imgfileFromRADfile(RAD_file, prefix):
    """ Transfer RAD filename to gt filename """
    RAD_file_spec = RAD_file.split("RAD")[-1]
    gt_file = os.path.join(prefix, "stereo_image") + RAD_file_spec.replace("npy", "jpg")
    return gt_file

def readRAD(filename):
    """ read input RAD matrices """
    if os.path.exists(filename):
        return np.load(filename)
    else:
        return None

def readRadarInstances(pickle_file):
    """ read output radar instances. """
    if os.path.exists(pickle_file):
        with open(pickle_file, "rb") as f:
            radar_instances = pickle.load(f)
        if len(radar_instances['classes']) == 0:
            radar_instances = None
    else:
        radar_instances = None
    return radar_instances

def readStereoLeft(img_filename):
    """ read stereo left image for verification. """
    if os.path.exists(img_filename):
        stereo_image = cv2.imread(img_filename)
        left_image = stereo_image[:, :stereo_image.shape[1]//2, ...][..., ::-1]
        return left_image
    else:
        return None

def readSingleImage(img_filename):
    """ read stereo left image for verification. """
    if os.path.exists(img_filename):
        image = cv2.imread(img_filename)
        image = image[..., ::-1]
        return image
    else:
        return None

max_stars_repo_name: BaiLiping/RADDet

max_issues_repo_name: BaiLiping/RADDet

max_forks_repo_name: BaiLiping/RADDet


text:
###############################################################################
# WaterTAP Copyright (c) 2021, The Regents of the University of California,
# through Lawrence Berkeley National Laboratory, Oak Ridge National
# Laboratory, National Renewable Energy Laboratory, and National Energy
# Technology Laboratory (subject to receipt of any required approvals from
# the U.S. Dept. of Energy). All rights reserved.
#
# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
# information, respectively. These files are also available online at the URL
# "https://github.com/watertap-org/watertap/"
#
###############################################################################
"""
Tests for zero-order brine concentrator model
"""
import pytest
from io import StringIO

from pyomo.environ import (
    check_optimal_termination, ConcreteModel, Constraint, value, Var, Param)
from pyomo.util.check_units import assert_units_consistent

from idaes.core import FlowsheetBlock
from idaes.core.util import get_solver
from idaes.core.util.model_statistics import degrees_of_freedom
from idaes.core.util.testing import initialization_tester

from watertap.unit_models.zero_order import BrineConcentratorZO
from watertap.core.wt_database import Database
from watertap.core.zero_order_properties import WaterParameterBlock

solver = get_solver()

class TestBrineConcentratorZO_w_o_default_removal:
    @pytest.fixture(scope="class")
    def model(self):
        m = ConcreteModel()
        m.db = Database()

        m.fs = FlowsheetBlock(default={"dynamic": False})
        m.fs.params = WaterParameterBlock(
            default={"solute_list": ["tds",]})

        m.fs.unit = BrineConcentratorZO(default={
            "property_package": m.fs.params,
            "database": m.db})

        m.fs.unit.inlet.flow_mass_comp[0, "H2O"].fix(10000)
        m.fs.unit.inlet.flow_mass_comp[0, "tds"].fix(250)

        return m

    @pytest.mark.unit
    def test_build(self, model):
        assert model.fs.unit.config.database is model.db
        assert isinstance(model.fs.unit.power_consumption_constraint, Constraint)
        assert isinstance(model.fs.unit.power_consumption, Var)

    @pytest.mark.component
    def test_load_parameters(self, model):
        data = model.db.get_unit_operation_parameters("brine_concentrator")

        model.fs.unit.load_parameters_from_database()
        assert model.fs.unit.recovery_frac_mass_H2O[0].fixed
        assert model.fs.unit.recovery_frac_mass_H2O[0].value == \
            data["recovery_frac_mass_H2O"]["value"]

        for (t, j), v in model.fs.unit.removal_frac_mass_solute.items():
            assert v.fixed
            assert v.value == data["removal_frac_mass_solute"][j]["value"]

    @pytest.mark.component
    def test_degrees_of_freedom(self, model):
        assert degrees_of_freedom(model.fs.unit) == 0

    @pytest.mark.component
    def test_unit_consistency(self, model):
        assert_units_consistent(model.fs.unit)

    @pytest.mark.component
    def test_initialize(self, model):
        initialization_tester(model)

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_solve(self, model):
        results = solver.solve(model)

        # Check for optimal solution
        assert check_optimal_termination(results)

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_solution(self, model):
        assert (pytest.approx(9.005, rel=1e-5) ==
                value(model.fs.unit.properties_treated[0].flow_vol))
        assert (pytest.approx(0.555247, rel=1e-5) == value(
            model.fs.unit.properties_treated[0].conc_mass_comp["tds"]))
        assert (pytest.approx(196.78715, rel=1e-5) == value(
            model.fs.unit.properties_byproduct[0].conc_mass_comp["tds"]))
        assert (pytest.approx(855570.663, rel=1e-5) ==
                value(model.fs.unit.power_consumption[0]))
        assert (pytest.approx(23.186197, rel=1e-5) ==
                value(model.fs.unit.electricity_intensity[0]))

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_conservation(self, model):
        for j in model.fs.params.component_list:
            assert 1e-6 >= abs(value(
                model.fs.unit.inlet.flow_mass_comp[0, j] -
                model.fs.unit.treated.flow_mass_comp[0, j] -
                model.fs.unit.byproduct.flow_mass_comp[0, j]))

    @pytest.mark.component
    def test_report(self, model):
        stream = StringIO()

        model.fs.unit.report(ostream=stream)

        output = """
====================================================================================
Unit : fs.unit                                                             Time: 0.0
------------------------------------------------------------------------------------
    Unit Performance

    Variables:

    Key                                                : Value      : Fixed : Bounds
    Electricity intensity per Inlet Flowrate  (kWh/m3) :     23.186 : False : (None, None)
                                Power Consumption (kW) : 8.5557e+05 : False : (None, None)
                                  Solute Removal [tds] :    0.98000 :  True : (0, None)
                                        Water Recovery :    0.90000 :  True : (1e-08, 1.0000001)

------------------------------------------------------------------------------------
    Stream Table
                            Inlet  Treated  Byproduct
    Volumetric Flowrate    10.250  9.0050    1.2450
    Mass Concentration H2O 975.61  999.44    803.21
    Mass Concentration tds 24.390 0.55525    196.79
====================================================================================
"""

        assert output in stream.getvalue()

class Testbrine_concentratorZO_w_default_removal:
    @pytest.fixture(scope="class")
    def model(self):
        m = ConcreteModel()
        m.db = Database()

        m.fs = FlowsheetBlock(default={"dynamic": False})
        m.fs.params = WaterParameterBlock(
            default={"solute_list": ["tds", "foo"]})

        m.fs.unit = BrineConcentratorZO(default={
            "property_package": m.fs.params,
            "database": m.db})

        m.fs.unit.inlet.flow_mass_comp[0, "H2O"].fix(10000)
        m.fs.unit.inlet.flow_mass_comp[0, "tds"].fix(250)
        m.fs.unit.inlet.flow_mass_comp[0, "foo"].fix(1)

        return m

    @pytest.mark.unit
    def test_build(self, model):
        assert model.fs.unit.config.database is model.db
        assert isinstance(model.fs.unit.power_consumption_constraint, Constraint)
        assert isinstance(model.fs.unit.power_consumption, Var)

    @pytest.mark.component
    def test_load_parameters(self, model):
        data = model.db.get_unit_operation_parameters("brine_concentrator")
        model.fs.unit.load_parameters_from_database(use_default_removal=True)
        assert model.fs.unit.recovery_frac_mass_H2O[0].fixed
        assert model.fs.unit.recovery_frac_mass_H2O[0].value == \
            data["recovery_frac_mass_H2O"]["value"]
        for (t, j), v in model.fs.unit.removal_frac_mass_solute.items():
            assert v.fixed
            if j == "foo":
                assert v.value == data["default_removal_frac_mass_solute"]["value"]
            else:
                assert v.value == data["removal_frac_mass_solute"][j]["value"]

    @pytest.mark.component
    def test_degrees_of_freedom(self, model):
        assert degrees_of_freedom(model.fs.unit) == 0

    @pytest.mark.component
    def test_unit_consistency(self, model):
        assert_units_consistent(model.fs.unit)

    @pytest.mark.component
    def test_initialize(self, model):
        initialization_tester(model)

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_solve(self, model):
        results = solver.solve(model)

        # Check for optimal solution
        assert check_optimal_termination(results)

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_solution(self, model):
        assert (pytest.approx(9.006, rel=1e-5) ==
                value(model.fs.unit.properties_treated[0].flow_vol))
        assert (pytest.approx(0.555185, rel=1e-5) == value(
            model.fs.unit.properties_treated[0].conc_mass_comp["tds"]))
        assert (pytest.approx(0.1110371, rel=1e-5) == value(
            model.fs.unit.properties_treated[0].conc_mass_comp["foo"]))
        assert (pytest.approx(196.787149, rel=1e-5) == value(
            model.fs.unit.properties_byproduct[0].conc_mass_comp["tds"]))
        assert (pytest.approx(8.0321285e-09, rel=1e-5) == value(
            model.fs.unit.properties_byproduct[0].conc_mass_comp["foo"]))
        assert (pytest.approx(855649.563040, rel=1e-5) ==
                value(model.fs.unit.power_consumption[0]))
        assert (pytest.approx(23.186, rel=1e-5) ==
                value(model.fs.unit.electricity_intensity[0]))

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_conservation(self, model):
        for j in model.fs.params.component_list:
            assert 1e-6 >= abs(value(
                model.fs.unit.inlet.flow_mass_comp[0, j] -
                model.fs.unit.treated.flow_mass_comp[0, j] -
                model.fs.unit.byproduct.flow_mass_comp[0, j]))

    @pytest.mark.component
    def test_report(self, model):
        stream = StringIO()

        model.fs.unit.report(ostream=stream)

        output = """
====================================================================================
Unit : fs.unit                                                             Time: 0.0
------------------------------------------------------------------------------------
    Unit Performance

    Variables:

    Key                                                : Value      : Fixed : Bounds
    Electricity intensity per Inlet Flowrate  (kWh/m3) :     23.186 : False : (None, None)
                                Power Consumption (kW) : 8.5565e+05 : False : (None, None)
                                  Solute Removal [foo] :     0.0000 :  True : (0, None)
                                  Solute Removal [tds] :    0.98000 :  True : (0, None)
                                        Water Recovery :    0.90000 :  True : (1e-08, 1.0000001)

------------------------------------------------------------------------------------
    Stream Table
                             Inlet   Treated  Byproduct
    Volumetric Flowrate      10.251  9.0060      1.2450
    Mass Concentration H2O   975.51  999.33      803.21
    Mass Concentration tds   24.388 0.55519      196.79
    Mass Concentration foo 0.097551 0.11104  8.0321e-09
====================================================================================
"""
        assert output in stream.getvalue()

db = Database()
params = db._get_technology("brine_concentrator")

class Testbrine_concentratorZOsubtype:
    @pytest.fixture(scope="class")
    def model(self):
        m = ConcreteModel()

        m.fs = FlowsheetBlock(default={"dynamic": False})
        m.fs.params = WaterParameterBlock(
            default={"solute_list": ["tds"]})

        m.fs.unit = BrineConcentratorZO(default={
            "property_package": m.fs.params,
            "database": db})

        return m

    @pytest.mark.parametrize("subtype", [params.keys()])
    @pytest.mark.component
    def test_load_parameters(self, model, subtype):
        model.fs.unit.config.process_subtype = subtype
        data = db.get_unit_operation_parameters("brine_concentrator", subtype=subtype)

        model.fs.unit.load_parameters_from_database()

        for (t, j), v in model.fs.unit.removal_frac_mass_solute.items():
            assert v.fixed
            assert v.value == data["removal_frac_mass_solute"][j]["value"]

@pytest.mark.unit
def test_no_tds_in_solute_list_error():
    m = ConcreteModel()
    m.fs = FlowsheetBlock(default={"dynamic": False})
    m.fs.params = WaterParameterBlock(
        default={"solute_list": ["foo"]})

    with pytest.raises(KeyError,
                       match="TDS must be included in the solute list for determining "
                             "electricity intensity and power consumption of the brine concentrator unit."):
        m.fs.unit = BrineConcentratorZO(default={"property_package": m.fs.params,
                                                 "database": db})

max_stars_repo_name: jalving/watertap

max_issues_repo_name: jalving/watertap

max_forks_repo_name: jalving/watertap

reason: from pyomo


text:
import os
import sys
sys.path.append(os.path.abspath(os.path.dirname(__file__) + '/' + '../..'))

import numpy as np
from numpy.testing import assert_allclose
import pandas as pd

def tests_positive_infected():
    # Check that the time series for
    # infected individuals is always weakly positive
    df = pd.read_csv('../../../fixed_revisions/derived_data/dataset/dataset.csv')
    for days_infectious in range(5, 10 + 1):
        if days_infectious == 5:
            mask = df['infected_{}'.format(days_infectious)] < 0.0
        else:
            mask_temp = df['infected_{}'.format(days_infectious)] < 0.0
            mask = mask | mask_temp
    test = df.loc[mask, ].shape[0]
    print (df.loc[mask, ])
    assert_allclose(0, test)

def test_new_cases():
    # Check that new cases are always
    # weakly positive
    df = pd.read_csv('../../../fixed_revisions/derived_data/dataset/dataset.csv')
    mask = df['new_cases'] < 0
    test = df.loc[mask, ].shape[0]
    print (df.loc[mask, ])
    assert_allclose(0, test)

def tests_growth_rates():
    # Check that the growth rate of infected
    # individuals is equal to (-gamma) when
    # new cases are zero
    df = pd.read_csv('../../../fixed_revisions/derived_data/dataset/dataset.csv')
    for days_infectious in range(5, 10 + 1):
        mask = (df['new_cases'] == 0.0) & ~np.isnan((df['gr_infected_{}'.format(days_infectious)]))
        assert_allclose(df.loc[mask, 'gr_infected_{}'.format(days_infectious)], - 1 / float(days_infectious))

def tests_growth_rates_2():
    # Check that the growth rate of infected
    # individuals is bounded below by gamma
    # when new cases are positive
    df = pd.read_csv('../../../fixed_revisions/derived_data/dataset/dataset.csv')
    for days_infectious in range(5, 10 + 1):
        mask = (df['new_cases'] > 0.0) & df['gr_infected_{}'.format(days_infectious)] <= -1 / float(days_infectious)
        assert_allclose(df.loc[mask,].shape[0], 0)

def tests_growth_rates_zero_new_cases():
    # Additional check for growth rates
    # with zero new cases
    df = pd.read_csv('../../../fixed_revisions/derived_data/dataset/dataset.csv')
    mask = df['new_cases'] == 0.0
    for days_infectious in range(5, 10 + 1):
        df['temp'] = df['infected_{}'.format(days_infectious)] / df.groupby('Country/Region').shift()['infected_{}'.format(days_infectious)] - 1
        df_temp = df.loc[:, ['temp', 'gr_infected_{}'.format(days_infectious)]].dropna()
        assert_allclose(df_temp['temp'], df_temp['gr_infected_{}'.format(days_infectious)])

reason: import numpy,from numpy


text:
from timeit import default_timer as timer
import numpy as np
from .writer import ExperimentWriter
from .experiment import Experiment

class SingleEnvExperiment(Experiment):
    '''An Experiment object for training and testing agents that interact with one environment at a time.'''
    def __init__(
            self,
            agent,
            env,
            logdir='runs',
            quiet=False,
            render=False,
            write_loss=True
    ):
        super().__init__(self._make_writer(logdir, agent.__name__, env.name, write_loss), quiet)
        self._agent = agent(env, self._writer)
        self._env = env
        self._render = render
        self._frame = 1
        self._episode = 1

        if render:
            self._env.render(mode="human")

    @property
    def frame(self):
        return self._frame

    @property
    def episode(self):
        return self._episode

    def train(self, frames=np.inf, episodes=np.inf):
        while not self._done(frames, episodes):
            self._run_training_episode()

    def test(self, episodes=100):
        returns = []
        for episode in range(episodes):
            episode_return = self._run_test_episode()
            returns.append(episode_return)
            self._log_test_episode(episode, episode_return)
        self._log_test(returns)
        return returns

    def _run_training_episode(self):
        # initialize timer
        start_time = timer()
        start_frame = self._frame

        # initialize the episode
        state = self._env.reset()
        action = self._agent.act(state)
        returns = 0

        # loop until the episode is finished
        while not state.done:
            if self._render:
                self._env.render()
            state = self._env.step(action)
            action = self._agent.act(state)
            returns += state.reward
            self._frame += 1

        # stop the timer
        end_time = timer()
        fps = (self._frame - start_frame) / (end_time - start_time)

        # log the results
        self._log_training_episode(returns, fps)

        # update experiment state
        self._episode += 1

    def _run_test_episode(self):
        # initialize the episode
        state = self._env.reset()
        action = self._agent.eval(state)
        returns = 0

        # loop until the episode is finished
        while not state.done:
            if self._render:
                self._env.render()
            state = self._env.step(action)
            action = self._agent.eval(state)
            returns += state.reward

        return returns

    def _done(self, frames, episodes):
        return self._frame > frames or self._episode > episodes

    def _make_writer(self, logdir, agent_name, env_name, write_loss):
        return ExperimentWriter(self, agent_name, env_name, loss=write_loss, logdir=logdir)

max_stars_repo_name: ezelikman/autonomous-learning-library

max_issues_repo_name: sadiakhaf/autonomous-learning-library

max_forks_repo_name: sadiakhaf/autonomous-learning-library


text:
import numpy as np
from motpy.metrics import angular_similarity, calculate_iou

from utils import assert_almost_equal, current_milli_time

def test_ious():
    b1 = [[10, 20]]
    b2 = [[10, 21], [30, 40]]
    iou_1d = calculate_iou(b1, b2, dim=1)
    assert_almost_equal(iou_1d, np.array([[0.9091, 0.]]))

    b1 = [[20.1, 20.1, 30.1, 30.1], [15, 15, 25, 25]]
    b2 = [[10, 10, 20, 20]]
    iou_2d = calculate_iou(b1, b2, dim=2)
    assert_almost_equal(iou_2d, [[0], [0.1429]])

    b1 = [[10, 10, 10, 20, 20, 20]]
    b2 = [[10, 11, 10.2, 21, 19.9, 20.3],
          [30, 30, 30, 90, 90, 90]]
    iou_3d = calculate_iou(b1, b2, dim=3)
    assert_almost_equal(iou_3d, [[0.7811, 0]])

def test_ious_large():
    # benchmarking test
    t = current_milli_time()
    for _ in range(100):
        b1 = np.random.randn(30, 4)
        b2 = np.random.randn(20, 4)
        iou = calculate_iou(b1, b2, dim=2)
        assert iou.shape == (30, 20)

    tdiff = current_milli_time() - t
    assert tdiff < 100

def test_angular_similarity():
    a = [[10, 10], [0, 100]]
    b = [[0, 10], [-10, -10]]

    mat = angular_similarity(a, b)
    mat_exp = np.array([[0.8536, 0.],
                        [1., 0.1464]])

    assert_almost_equal(mat, mat_exp)

max_stars_repo_name: brdg-ai/motpy

max_issues_repo_name: brdg-ai/motpy

max_forks_repo_name: brdg-ai/motpy


text:
import logging
import os

import aircv as ac
import cv2
import numpy as np
from PIL import Image

from constant import PC_PROJECT_ROOT, SCREEN_PATH, PC_RECONGNIZE_TEST, getRecongnizeTarget,getTargetImgs

from logger import logger as logging

def matchImg(src_path, obj_path, confidencevalue=0.8):  # imgsrc=原始图像，imgobj=待查找的图片
    imsrc = cv2.imdecode(np.fromfile(src_path,dtype=np.uint8),-1)
    imobj = cv2.imdecode(np.fromfile(obj_path,dtype=np.uint8),-1)
    match_result = ac.find_template(imsrc, imobj, confidencevalue)
    if match_result is not None:
        match_result['shape'] = (imsrc.shape[1], imsrc.shape[0])  # 0为高，1为宽
    return match_result

def find_img_position(debugEveryOne=False):
    if debugEveryOne and not os.path.exists(PC_RECONGNIZE_TEST):
        os.makedirs(PC_RECONGNIZE_TEST)

    for template in getTargetImgs():
        res = matchImg('{}{}'.format(PC_PROJECT_ROOT,SCREEN_PATH),template,confidencevalue=0.9)
        # {'confidence': 0.5435812473297119, 'rectangle': ((394, 384), (394, 416), (450, 384), (450, 416)), 'result': (422.0, 400.0)
        # confidence：匹配相似率
        #
        # rectangle：匹配图片在原始图像上四边形的坐标
        #
        # result：匹配图片在原始图片上的中心坐标点，也就是我们要找的点击点
        # 如果结果匹配到的confidence小于入参传递的相似度confidence，则会返回None，不返回字典
        if res is not None:
            x, y = res['result']
            logging.debug('match result =  {}'.format(res['confidence']))
            logging.debug('match position =  ({},{})'.format(x, y))

            if debugEveryOne is True:
                logging.debug('match Img :{}'.format(template))
                # show the rect of find subImage
                rect = res['rectangle']
                crop_file = template.replace(getRecongnizeTarget(),PC_RECONGNIZE_TEST)
                if os.path.exists(crop_file):
                    os.remove(crop_file)
                Image.open('{}{}'.format(PC_PROJECT_ROOT,SCREEN_PATH)).crop((rect[0][0],
                                                                    rect[0][1],
                                                                     rect[3][0],
                                                                     rect[3][1],)).save(crop_file)
            else:
                return template,x, y
    return None

max_stars_repo_name: senjoeson/kog-money

max_issues_repo_name: senjoeson/kog-money

max_forks_repo_name: senjoeson/kog-money


text:
"""
Traffic Flow Prediction with Neural Networks(SAEs、LSTM、GRU).
"""
import math
import warnings
import numpy as np
import pandas as pd
from data.data import process_data, process_data_combine
from keras.models import load_model
from keras.utils.vis_utils import plot_model
import sklearn.metrics as metrics
import matplotlib
import matplotlib.pyplot as plt

warnings.filterwarnings("ignore")

def MAPE(y_true, y_pred):
    """Mean Absolute Percentage Error
    Calculate the mape.

    # Arguments
        y_true: List/ndarray, ture data.
        y_pred: List/ndarray, predicted data.
    # Returns
        mape: Double, result data for train.
    """

    y = [x for x in y_true if x > 0]
    y_pred = [y_pred[i] for i in range(len(y_true)) if y_true[i] > 0]

    num = len(y_pred)
    sums = 0

    for i in range(num):
        tmp = abs(y[i] - y_pred[i]) / y[i]
        sums += tmp

    mape = sums * (100 / num)

    return mape

def eva_regress(y_true, y_pred):
    """Evaluation
    evaluate the predicted resul.

    # Arguments
        y_true: List/ndarray, ture data.
        y_pred: List/ndarray, predicted data.
    """

    mape = MAPE(y_true, y_pred)
    vs = metrics.explained_variance_score(y_true, y_pred)
    mae = metrics.mean_absolute_error(y_true, y_pred)
    mse = metrics.mean_squared_error(y_true, y_pred)
    r2 = metrics.r2_score(y_true, y_pred)
    print('explained_variance_score:%f' % vs)
    print('mape:%f%%' % mape)
    print('mae:%f' % mae)
    print('mse:%f' % mse)
    print('rmse:%f' % math.sqrt(mse))
    print('r2:%f' % r2)
    return vs, mape, mae, mse, math.sqrt(mse), r2

def plot_results(y_true, y_preds, names):
    """Plot
    Plot the true data and predicted data.

    # Arguments
        y_true: List/ndarray, ture data.
        y_pred: List/ndarray, predicted data.
        names: List, Method names.
    """
    # d = '2016-3-4 00:00'
    # x = pd.date_range(d, periods=288, freq='5min')

    x = np.array(range(288))
    time_list = ['00:00', '02:00', '04:00', '06:00', '08:00',
                 '12:00', '14:00', '16:00', '18:00', '20:00', '22:00', '24:00']
    aa = np.linspace(0, 288, 12)

    fig = plt.figure()

    ax = fig.add_subplot(111)

    # ax.plot(x, y_true, label='True Data')
    for name, y_pred in zip(names, y_preds):
        ax.plot(x, y_pred, label=name)

    plt.legend()
    plt.grid(True)

    plt.xticks(aa, time_list, rotation=0)

    plt.xlabel('Time of Day')
    plt.ylabel('Flow')

    fig.autofmt_xdate()

    plt.show()

def main():
    lstm = load_model('model/lstm.h5')
    # lstm = load_model('model/lstm-no-drop.h5')
    gru = load_model('model/gru.h5')
    saes = load_model('model/saes.h5')
    models = [lstm, gru, saes]
    names = ['LSTM', 'gru', 'SAEs']

    # # 12 24 46尺度的LSTM
    # lstm12 = load_model('model/lstm_12.h5')
    # lstm24 = load_model('model/lstm_24.h5')
    # lstm36 = load_model('model/lstm_36.h5')

    # models = [lstm12, lstm24, lstm36]
    # names = ['LSTM_12', 'LSTM_24' ,'LSTM_36']

    lag = 24
    file1 = 'data/train.csv'
    file2 = 'data/test.csv'
    _, _, X_test12, y_test12, _ = process_data_combine(file1, file2, 12, 0)
    _, _, X_test24, y_test24, scaler = process_data_combine(file1, file2, 24, 0)
    _, _, X_test36, y_test36, _ = process_data_combine(file1, file2, 36, 0)

    num_lane = len(X_test12)

    X_test = X_test24
    y_test = y_test24

    for a in range(1):
        # 融合模型
        evs, mapes, maes, mses, rmses, r2s = np.zeros(num_lane), np.zeros(num_lane), np.zeros(num_lane), \
                                             np.zeros(num_lane), np.zeros(num_lane), \
                                             np.zeros(num_lane)  # explained_variance_score
        # y_preds = []
        for j in range(num_lane):
            stack_predict = []
            # stack_predict = np.zeros(288)

            for name, model in zip(names, models):

                if name == 'SAEs':
                    X_test[j] = np.reshape(X_test[j], (X_test[j].shape[0], X_test[j].shape[1]))
                else:
                    X_test[j] = np.reshape(X_test[j], (X_test[j].shape[0], X_test[j].shape[1], 1))

                predicted = model.predict(X_test[j])
                predicted = scaler.inverse_transform(predicted.reshape(-1, 1)).reshape(1, -1)[0]
                stack_predict.append(predicted)
                # y_preds.append(predicted[:288])
                # print(name)
            predicted = (7 * stack_predict[0] + 0 * stack_predict[1] + 3 * stack_predict[2]) / 10
            # predicted = (10 * stack_predict[0] + 0 * stack_predict[1] + 0 * stack_predict[2]) / 10
            y_test[j] = scaler.inverse_transform(y_test[j].reshape(-1, 1)).reshape(1, -1)[0]
            stack_predict.append(predicted)
            stack_predict.pop(1)

            name2 = ['LSTM', 'SAEs', 'LSTM-SAEs']
            plot_results(y_test[: 288], stack_predict, name2)

            vs, mape, mae, mse, rmse, r2 = eva_regress(y_test[j], predicted)
            evs[j] = vs
            mapes[j] = mape
            maes[j] = mae
            mses[j] = mse
            rmses[j] = rmse
            r2s[j] = r2
            # print("完成了" + str(j))

        print(str(a), str(10 - a), evs.mean(), mapes.mean(), maes.mean(), mses.mean(), rmses.mean(), r2s.mean(),
              r2s.mean() / evs.mean())

if __name__ == '__main__':
    main()