text:
import numpy as np
import autograd.numpy as anp
#from scipy.special.orthogonal import _IntegerType
from .pyerrors import *
from .fits import standard_fit
from .linalg import *
from .roots import find_root
from matplotlib import pyplot as plt
from matplotlib.ticker import NullFormatter  # useful for `logit` scale
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
#import PySimpleGUI as sg
import matplotlib

class Corr:
    """The class for a correlator (time dependent sequence of pe.Obs).

    Everything, this class does, can be achieved using lists or arrays of Obs.
    But it is simply more convenient to have a dedicated object for correlators.
    One often wants to add or multiply correlators of the same length at every timeslice and it is inconvinient
    to iterate over all timeslices for every operation. This is especially true, when dealing with smearing matrices.

    The correlator can have two types of content: An Obs at every timeslice OR a GEVP
    smearing matrix at every timeslice. Other dependency (eg. spacial) are not supported.

    """

    def __init__(self, data_input, padding_front=0, padding_back=0,prange=None):
        #All data_input should be a list of things at different timeslices. This needs to be verified

        if not (isinstance(data_input, list)):
            raise TypeError('Corr__init__ expects a list of timeslices.')
        # data_input can have multiple shapes. The simplest one is a list of Obs.
        #We check, if this is the case
        if all([isinstance(item, Obs) for item in data_input]):
            self.content=[np.asarray([item]) for item in data_input]
            #Wrapping the Obs in an array ensures that the data structure is consistent with smearing matrices.
            self.N = 1 # number of smearings

        #data_input in the form [np.array(Obs,NxN)]
        elif all([isinstance(item,np.ndarray) or item==None for item in data_input]) and any([isinstance(item,np.ndarray)for item in data_input]):
            self.content = data_input

            noNull=[a for a in self.content if not (a is None)] #To check if the matrices are correct for all undefined elements
            self.N = noNull[0].shape[0]
            # The checks are now identical to the case above
            if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]:
                raise Exception("Smearing matrices are not NxN")
            if  (not all([item.shape == noNull[0].shape for item in noNull])):
                raise Exception("Items in data_input are not of identical shape." + str(noNull))
        else: # In case its a list of something else.
            raise Exception ("data_input contains item of wrong type")

        self.tag = None

        #We now apply some padding to our list. In case that our list represents a correlator of length T but is not defined at every value.
        #An undefined timeslice is represented by the None object
        self.content = [None] * padding_front + self.content + [None] * padding_back
        self.T = len(self.content) #for convenience: will be used a lot

        #The attribute "prange" [start,end] marks a range of two timeslices.
        #This is useful for keeping track of plateaus and fitranges.
        #The prange can be inherited from other Corrs, if the operation should not alter a chosen range eg. multiplication with a constant.
        if not prange is None:
            self.prange=prange

        self.gamma_method()

    def gamma_method(self):
        for item in self.content:
            if not(item is None):
                if self.N == 1:
                    item[0].gamma_method()
                else:
                    for i in range(self.N):
                        for j in range(self.N):
                            item[i,j].gamma_method()

    #We need to project the Correlator with a Vector to get a single value at each timeslice.
    #The method can use one or two vectors.
    #If two are specified it returns v1@G@v2 (the order might be very important.)
    #By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to
    def projected(self, vector_l=None, vector_r=None):
        if self.N == 1:
            raise Exception("Trying to project a Corr, that already has N=1.")
            #This Exception is in no way necessary. One could just return self
            #But there is no scenario, where a user would want that to happen and the error message might be more informative.

        self.gamma_method()

        if vector_l is None:
            vector_l,vector_r=np.asarray([1.] + (self.N - 1) * [0.]),np.asarray([1.] + (self.N - 1) * [0.])
        elif(vector_r is None):
            vector_r=vector_l

        if not vector_l.shape == vector_r.shape == (self.N,):
            raise Exception("Vectors are of wrong shape!")

        #We always normalize before projecting! But we only raise a warning, when it is clear, they where not meant to be normalized.
        if (not(0.95 < vector_r@vector_r < 1.05)) or (not(0.95 < vector_l@vector_l < 1.05)):
            print("Vectors are normalized before projection!")

        vector_l,vector_r = vector_l / np.sqrt((vector_l@vector_l)), vector_r / np.sqrt(vector_r@vector_r)

        newcontent = [None if (item is None) else np.asarray([vector_l.T@item@vector_r]) for item in self.content]
        return Corr(newcontent)

    def sum(self):
        return np.sqrt(self.N)*self.projected(np.ones(self.N))
    #For purposes of debugging and verification, one might want to see a single smearing level. smearing will return a Corr at the specified i,j. where both are integers 0<=i,j<N.
    def smearing(self, i, j):
        if self.N == 1:
            raise Exception("Trying to pick smearing from projected Corr")
        newcontent=[None if(item is None) else item[i, j] for item in self.content]
        return Corr(newcontent)

    #Obs and Matplotlib do not play nicely
    #We often want to retrieve x,y,y_err as lists to pass them to something like pyplot.errorbar
    def plottable(self):
        if self.N != 1:
            raise Exception("Can only make Corr[N=1] plottable") #We could also autoproject to the groundstate or expect vectors, but this is supposed to be a super simple function.
        x_list = [x for x in range(self.T) if (not self.content[x] is None)]
        y_list = [y[0].value for y in self.content if (not y is None)]
        y_err_list = [y[0].dvalue for y in self.content if (not y is None)]

        return x_list, y_list, y_err_list

    #symmetric returns a Corr, that has been symmetrized.
    #A symmetry checker is still to be implemented
    #The method will not delete any redundant timeslices (Bad for memory, Great for convenience)
    def symmetric(self):

        if self.T%2 != 0:
            raise Exception("Can not symmetrize odd T")

        newcontent = [self.content[0]]
        for t in range(1, self.T):
            if (self.content[t] is None) or (self.content[self.T - t] is None):
                newcontent.append(None)
            else:
                newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
        if(all([x is None for x in newcontent])):
            raise Exception("Corr could not be symmetrized: No redundant values")

        return Corr(newcontent, prange=self.prange if hasattr(self,"prange") else None)

    def anti_symmetric(self):

        if self.T%2 != 0:
            raise Exception("Can not symmetrize odd T")

        newcontent=[self.content[0]]
        for t in range(1, self.T):
            if (self.content[t] is None) or (self.content[self.T - t] is None):
                newcontent.append(None)
            else:
                newcontent.append(0.5 * (self.content[t] - self.content[self.T - t]))
        if(all([x is None for x in newcontent])):
            raise Exception("Corr could not be symmetrized: No redundant values")
        return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)

    #This method will symmetrice the matrices and therefore make them positive definit.
    def smearing_symmetric(self):
        if self.N > 1:
            transposed = [None if (G is None) else G.T for G in self.content]
            return 0.5 * (Corr(transposed)+self)
        if self.N == 1:
            raise Exception("Trying to symmetrize a smearing matrix, that already has N=1.")

    #We also include a simple GEVP method based on Scipy.linalg
    def GEVP(self, t0, ts,state=1):
        if (self.content[t0] is None) or (self.content[ts] is None):
            raise Exception("Corr not defined at t0/ts")
        G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
        for i in range(self.N):
            for j in range(self.N):
                G0[i, j] = self.content[t0][i, j].value
                Gt[i, j] = self.content[ts][i, j].value

        sp_val,sp_vec = scipy.linalg.eig(Gt, G0)
        sp_vec=sp_vec[:,np.argsort(sp_val)[-state]] #we only want the eigenvector belonging to the selected state
        sp_vec = sp_vec/np.sqrt(sp_vec@sp_vec)
        return sp_vec

    def Eigenvalue(self,t0,state=1):
        G=self.smearing_symmetric()
        G0=G.content[t0]
        L = mat_mat_op(anp.linalg.cholesky, G0)
        Li = mat_mat_op(anp.linalg.inv, L)
        LT=L.T
        LTi=mat_mat_op(anp.linalg.inv, LT)
        newcontent=[]
        for t in range(self.T):
            Gt=G.content[t]
            M=Li@Gt@LTi
            eigenvalues = eigh(M)[0]
            #print(eigenvalues)
            eigenvalue=eigenvalues[-state]
            newcontent.append(eigenvalue)
        return Corr(newcontent)

    def roll(self, dt):
        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))

    def deriv(self, symmetric=True): #Defaults to symmetric derivative
        if not symmetric:
            newcontent = []
            for t in range(self.T - 1):
                if (self.content[t] is None) or (self.content[t+1] is None):
                    newcontent.append(None)
                else:
                    newcontent.append(self.content[t + 1] - self.content[t])
            if(all([x is None for x in newcontent])):
                raise Exception("Derivative is undefined at all timeslices")
            return Corr(newcontent, padding_back=1)
        if symmetric:
            newcontent = []
            for t in range(1, self.T-1):
                if (self.content[t-1] is None) or (self.content[t+1] is None):
                    newcontent.append(None)
                else:
                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
            if(all([x is None for x in newcontent])):
                raise Exception('Derivative is undefined at all timeslices')
            return Corr(newcontent, padding_back=1, padding_front=1)

    def second_deriv(self):
        newcontent = []
        for t in range(1, self.T-1):
            if (self.content[t-1] is None) or (self.content[t+1] is None):
                newcontent.append(None)
            else:
                newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
        if(all([x is None for x in newcontent])):
            raise Exception("Derivative is undefined at all timeslices")
        return Corr(newcontent, padding_back=1, padding_front=1)

    def m_eff(self, variant='log', guess=1.0):
        """Returns the effective mass of the correlator as correlator object

        Parameters
        ----------
        variant -- log: uses the standard effective mass log(C(t) / C(t+1))
                   periodic : Solves C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m. See, e.g., arXiv:1205.5380
        guess -- guess for the root finder, only relevant for the root variant
        """
        if self.N != 1:
            raise Exception('Correlator must be projected before getting m_eff')
        if variant is 'log':
            newcontent = []
            for t in range(self.T - 1):
                if (self.content[t] is None) or (self.content[t + 1] is None):
                    newcontent.append(None)
                else:
                    newcontent.append(self.content[t] / self.content[t + 1])
            if(all([x is None for x in newcontent])):
                raise Exception('m_eff is undefined at all timeslices')

            return np.log(Corr(newcontent, padding_back=1))

        elif variant is 'periodic':
            newcontent = []
            for t in range(self.T - 1):
                if (self.content[t] is None) or (self.content[t + 1] is None):
                    newcontent.append(None)
                else:
                    func = lambda x, d : anp.cosh(x * (t - self.T / 2)) / anp.cosh(x * (t + 1 - self.T / 2)) - d
                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], func, guess=guess)))
            if(all([x is None for x in newcontent])):
                raise Exception('m_eff is undefined at all timeslices')

            return Corr(newcontent, padding_back=1)
        elif variant is 'arccosh':
            newcontent=[]
            for t in range(1,self.T-1):
                if (self.content[t] is None) or (self.content[t+1] is None)or (self.content[t-1] is None):
                    newcontent.append(None)
                else:
                    newcontent.append((self.content[t+1]+self.content[t-1])/(2*self.content[t]))
            if(all([x is None for x in newcontent])):
                raise Exception("m_eff is undefined at all timeslices")
            return np.arccosh(Corr(newcontent,padding_back=1,padding_front=1))

        else:
            raise Exception('Unkown variant.')

    #We want to apply a pe.standard_fit directly to the Corr using an arbitrary function and range.
    def fit(self, function, fitrange=None, silent=False, **kwargs):
        if self.N != 1:
            raise Exception("Correlator must be projected before fitting")

        #The default behaviour is:
            #1 use explicit fitrange
            # if none is provided, use the range of the corr
            # if this is also not set, use the whole length of the corr (This could come with a warning!)

        if fitrange is None:
            if hasattr(self,"prange"):
                fitrange=self.prange
            else:
                fitrange=[0, self.T]

        xs = [x for x in range(fitrange[0], fitrange[1]) if not self.content[x] is None]
        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1]) if not self.content[x] is None]
        result = standard_fit(xs, ys, function, silent=silent, **kwargs)
        if isinstance(result, list):
            [item.gamma_method() for item in result if isinstance(item,Obs)]
        elif isinstance(result, dict):
            [item.gamma_method() for item in result['fit_parameters'] if isinstance(item,Obs)]
        else:
            raise Exception('Unexpected fit result.')
        return result

    #we want to quickly get a plateau
    def plateau(self, plateau_range=None, method="fit"):
        if not plateau_range:
            if hasattr(self,"prange"):
                plateau_range=self.prange
            else:
                raise Exception("no plateau range provided")
        if self.N != 1:
            raise Exception("Correlator must be projected before getting a plateau.")
        if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1])])):
                raise Exception("plateau is undefined at all timeslices in plateaurange.")
        if method == "fit":
            def const_func(a, t):
                return a[0] # At some point pe.standard fit had an issue with single parameter fits. Being careful does not hurt
            return self.fit(const_func,plateau_range)[0]
        elif method in ["avg","average","mean"]:
            returnvalue= np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1]+1] if not item is None])
            returnvalue.gamma_method()
            return returnvalue

        else:
            raise Exception("Unsupported plateau method: " + method)

    def set_prange(self, prange):
        if not len(prange)==2:
            raise Exception("prange must be a list or array with two values")
        if not ((isinstance(prange[0],int)) and (isinstance(prange[1],int))):
            raise Exception("Start and end point must be integers")
        if not (0<=prange[0]<=self.T and 0<=prange[1]<=self.T and prange[0]<prange[1]  ):
            raise Exception("Start and end point must define a range in the interval 0,T")

        self.prange = prange
        return

    # Plotting routine for correlator
    def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, save=None, ylabel=None):
        """Plots the correlator, uses tag as label if available.

        Parameters
        ----------
        x_range -- list of two values, determining the range of the x-axis e.g. [4, 8]
        comp -- Correlator or list of correlators which are plotted for comparison.
        logscale -- Sets y-axis to logscale
        save -- path to file in which the figure should be saved
        """
        if self.N!=1:
            raise Exception("Correlator must be projected before plotting")
        if  x_range is None:
            x_range = [0, self.T]

        fig = plt.figure()
        ax1 = fig.add_subplot(111)

        x,y,y_err=self.plottable()
        ax1.errorbar(x, y, y_err, label=self.tag)
        if logscale:
            ax1.set_yscale('log')
        else:
            # we generate ylim instead of using autoscaling.
            if y_range is None:
                try:
                    y_min=min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]:x_range[1]] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
                    y_max=max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]:x_range[1]] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
                except:
                    pass
            else:
                ax1.set_ylim(y_range)
        if comp:
            if isinstance(comp, Corr) or isinstance(comp, list):
                for corr in comp if isinstance(comp, list) else [comp]:
                    x,y,y_err=corr.plottable()
                    plt.errorbar(x, y, y_err, label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
            else:
                raise Exception('comp must be a correlator or a list of correlators.')

        if plateau:
            if isinstance(plateau, Obs):
                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=plateau.__repr__()[4:-1])
                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
            else:
                raise Exception('plateau must be an Obs')
        if hasattr(self,"prange"):
            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')

        if fit_res:
            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
            ax1.plot(x_samples,
                    fit_res['fit_function']([o.value for o in fit_res['fit_parameters']], x_samples)
                    , ls='-', marker=',', lw=2)

        ax1.set_xlabel(r'$x_0 / a$')
        if ylabel:
            ax1.set_ylabel(ylabel)
        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])

        handles, labels = ax1.get_legend_handles_labels()
        if labels:
            legend = ax1.legend()
        plt.draw()

        if save:
            if isinstance(save, str):
                fig.savefig(save)
            else:
                raise Exception('safe has to be a string.')

        return

    def dump(self,filename):
        dump_object(self,filename)
        return

    def print(self, range=[0, None]):
        print(self.__repr__(range))

    def __repr__(self, range=[0, None]):
        if range[1]:
            range[1] += 1
        content_string = 'x0/a\tCorr(x0/a)\n------------------\n'
        for i, sub_corr in enumerate(self.content[range[0]:range[1]]):
            if sub_corr is None:
                content_string += str(i + range[0]) + '\n'
            else:
                content_string += str(i + range[0])
                for element in sub_corr:
                    content_string += '\t' + element.__repr__()[4:-1]
                content_string += '\n'
        return content_string
    def __str__(self):
        return self.__repr__()
        #return ("Corr[T="+str(self.T)+" , N="+str(self.N)+" , content="+str([o[0] for o in [o for o in self.content]])+"]")

    #We define the basic operations, that can be performed with correlators.
    #While */+- get defined here, they only work for Corr*Obs and not Obs*Corr.
    #This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception.
    #One could try and tell Obs to check if the y in __mul__ is a Corr and

    def __add__(self, y):
        if isinstance(y, Corr):
            if ((self.N!=y.N) or (self.T!=y.T) ):
                raise Exception("Addition of Corrs with different shape")
            newcontent=[]
            for t in range(self.T):
                if (self.content[t] is None) or (y.content[t] is None):
                    newcontent.append(None)
                else:
                    newcontent.append(self.content[t]+y.content[t])
            return Corr(newcontent)

        elif isinstance(y, Obs) or isinstance(y, int) or isinstance(y,float):
            newcontent=[]
            for t in range(self.T):
                if (self.content[t] is None):
                    newcontent.append(None)
                else:
                    newcontent.append(self.content[t]+y)
            return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
        else:
            raise TypeError("Corr + wrong type")

    def __mul__(self,y):
        if isinstance(y,Corr):
            if not((self.N==1 or y.N==1 or self.N==y.N) and self.T==y.T):
                raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T")
            newcontent=[]
            for t in range(self.T):
                if (self.content[t] is None) or (y.content[t] is None):
                    newcontent.append(None)
                else:
                    newcontent.append(self.content[t]*y.content[t])
            return Corr(newcontent)

        elif isinstance(y, Obs) or isinstance(y, int) or isinstance(y,float):
            newcontent=[]
            for t in range(self.T):
                if (self.content[t] is None):
                    newcontent.append(None)
                else:
                    newcontent.append(self.content[t]*y)
            return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
        else:
            raise TypeError("Corr * wrong type")

    def __truediv__(self,y):
        if isinstance(y,Corr):
            if not((self.N==1 or y.N==1 or self.N==y.N) and self.T==y.T):
                raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T")
            newcontent=[]
            for t in range(self.T):
                if (self.content[t] is None) or (y.content[t] is None):
                    newcontent.append(None)
                else:
                    newcontent.append(self.content[t]/y.content[t])
            #Here we set the entire timeslice to undefined, if one of the smearings has encountered an division by zero.
            #While this might throw away perfectly good values in other smearings, we will never have to check, if all values in our matrix are defined
            for t in range(self.T):
                if newcontent[t] is None:
                    continue
                if np.isnan(np.sum(newcontent[t]).value):
                    newcontent[t]=None

            if all([item is None for item in newcontent]):
                raise Exception("Division returns completely undefined correlator")

            return Corr(newcontent)

        elif isinstance(y, Obs):
            if y.value==0:
                raise Exception("Division by Zero will return undefined correlator")
            newcontent=[]
            for t in range(self.T):
                if (self.content[t] is None):
                    newcontent.append(None)
                else:
                    newcontent.append(self.content[t]/y)
            if hasattr(self,"prange"):
                newrange=self.prange
            else:
                newrange=None
            return Corr(newcontent,prange=newrange)

        elif isinstance(y, int) or isinstance(y,float):
            if y==0:
                raise Exception("Division by Zero will return undefined correlator")
            newcontent=[]
            for t in range(self.T):
                if (self.content[t] is None):
                    newcontent.append(None)
                else:
                    newcontent.append(self.content[t]/y)
            return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
        else:
            raise TypeError("Corr / wrong type")

    def __neg__(self):
        newcontent=[None if (item is None) else -1.*item for item in self.content]
        if hasattr(self,"prange"):
            newrange=self.prange
        else:
            newrange=None
        return Corr(newcontent,prange=newrange)

    def __sub__(self,y):
        return self +(-y)

    def __pow__(self, y):
        if isinstance(y, Obs) or isinstance(y,int) or isinstance(y,float):
            newcontent=[None if (item is None) else item**y for item in self.content]
            return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
        else:
            raise TypeError("type of exponent not supported")

    def __abs__(self):
        newcontent=[None if (item is None) else np.abs(item) for item in self.content]
        return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)

#The numpy functions:
    def sqrt(self):
        return self**0.5

    def log(self):
        newcontent=[None if (item is None) else np.log(item) for item in self.content]
        return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)

    def exp(self):
        newcontent=[None if (item is None) else np.exp(item) for item in self.content]
        return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)

    def sin(self):
        newcontent=[None if (item is None) else np.sin(item) for item in self.content]
        return Corr(newcontent)

    def cos(self):
        newcontent=[None if (item is None) else np.cos(item) for item in self.content]
        return Corr(newcontent)

    def tan(self):
        newcontent=[None if (item is None) else np.tan(item) for item in self.content]

        for t in range(self.T):
                if newcontent[t] is None:
                    continue
                if np.isnan(np.sum(newcontent[t]).value):
                    newcontent[t]=None
        if all([item is None for item in newcontent]):
                raise Exception("Operation returns completely undefined correlator")

        return Corr(newcontent)

    def sinh(self):
        newcontent=[None if (item is None) else np.sinh(item) for item in self.content]
        return Corr(newcontent)

    def cosh(self):
        newcontent=[None if (item is None) else np.cosh(item) for item in self.content]
        return Corr(newcontent)

    def tanh(self):
        newcontent=[None if (item is None) else np.tanh(item) for item in self.content]
        for t in range(self.T):
                if newcontent[t] is None:
                    continue
                if np.isnan(np.sum(newcontent[t]).value):
                    newcontent[t]=None
        if all([item is None for item in newcontent]):
                raise Exception("Operation returns completely undefined correlator")
        return Corr(newcontent)

    def arcsin(self):
        newcontent=[None if (item is None) else np.arcsin(item) for item in self.content]
        for t in range(self.T):
                if newcontent[t] is None:
                    continue
                if np.isnan(np.sum(newcontent[t]).value):
                    newcontent[t]=None
        if all([item is None for item in newcontent]):
                raise Exception("Operation returns completely undefined correlator")
        return Corr(newcontent)

    def arccos(self):
        newcontent=[None if (item is None) else np.arccos(item) for item in self.content]
        for t in range(self.T):
                if newcontent[t] is None:
                    continue
                if np.isnan(np.sum(newcontent[t]).value):
                    newcontent[t]=None
        if all([item is None for item in newcontent]):
                raise Exception("Operation returns completely undefined correlator")
        return Corr(newcontent)

    def arctan(self):
        newcontent=[None if (item is None) else np.arctan(item) for item in self.content]
        for t in range(self.T):
                if newcontent[t] is None:
                    continue
                if np.isnan(np.sum(newcontent[t]).value):
                    newcontent[t]=None
        if all([item is None for item in newcontent]):
                raise Exception("Operation returns completely undefined correlator")
        return Corr(newcontent)

    def arcsinh(self):
        newcontent=[None if (item is None) else np.arcsinh(item) for item in self.content]
        for t in range(self.T):
                if newcontent[t] is None:
                    continue
                if np.isnan(np.sum(newcontent[t]).value):
                    newcontent[t]=None
        if all([item is None for item in newcontent]):
                raise Exception("Operation returns completely undefined correlator")
        return Corr(newcontent)

    def arccosh(self):
        newcontent=[None if (item is None) else np.arccosh(item) for item in self.content]
        for t in range(self.T):
                if newcontent[t] is None:
                    continue
                if np.isnan(np.sum(newcontent[t]).value):
                    newcontent[t]=None
        if all([item is None for item in newcontent]):
                raise Exception("Operation returns completely undefined correlator")
        return Corr(newcontent)

    def arctanh(self):
        newcontent=[None if (item is None) else np.arctanh(item) for item in self.content]
        for t in range(self.T):
                if newcontent[t] is None:
                    continue
                if np.isnan(np.sum(newcontent[t]).value):
                    newcontent[t]=None
        if all([item is None for item in newcontent]):
                raise Exception("Operation returns completely undefined correlator")
        return Corr(newcontent)

    #right hand side operations (require tweak in main module to work)
    def __rsub__(self,y):
        return -self+y
    def __rmul__(self, y):
        return self * y
    def __radd__(self,y):
        return self + y

max_stars_repo_name: s-kuberski/pyerrors

max_issues_repo_name: s-kuberski/pyerrors

max_forks_repo_name: s-kuberski/pyerrors

reason: import numpy,from scipy


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

r'''This module contains a base class for modeling computation graphs.

Neural networks are really just a concise, computational way of describing a
mathematical model of a computation graph that operates on a particular set of
data.

At a high level, a neural network is a computation graph that describes a
parametric mapping

.. math::
   F_\theta: \mathcal{S} \to \mathcal{T}

between a source space :math:`\mathcal{S}` and a target space
:math:`\mathcal{T}`, using parameters :math:`\theta`. For example, suppose we
are processing vectors representing the MNIST handwritten digits. We could think
of :math:`\mathcal{S} = \mathbb{R}^{28 \times 28} = \mathbb{R}^{784}` (i.e., the
space of all 28×28 images), and for classifying the MNIST digits we could think
of :math:`\mathcal{T} = \mathbb{R}^{10}`.

This mapping is assumed to be fairly complex. If it were not -- if you could
capture the mapping using a simple expression like :math:`F_{\{a\}}(x) = ax^2`
-- then we would just use the expression directly and not need to deal with an
entire network. So if the mapping is complex, we will do a couple of things to
make our problem tractable. First, we will assume some structure for
:math:`F_\theta`. Second, we will fit our model to some set of data that we have
obtained, so that our parameters :math:`\theta` are tuned to the problem at
hand.

Graph structure
---------------

.. image:: _static/feedforward_layers.svg

The mapping :math:`F_\theta` is implemented in neural networks by assuming a
specific, layered form. Computation nodes -- also called units or (sometimes)
neurons -- are arranged in a :math:`k+1` partite graph, with layer :math:`k`
containing :math:`n_k` nodes. The number of input nodes in the graph is referred
to as :math:`n_0`.

Most layers are connected together using a set of weights. A **weight matrix**
:math:`W^k \in \mathbb{R}^{n_{k-1} \times n_k}` specifies the strength of the
connection between nodes in layer :math:`k` and those in layer :math:`k-1` --
all other pairs of nodes are typically not connected. Each layer of nodes also
typically has a **bias vector** that determines the offset of each node from the
origin. Together, the parameters :math:`\theta` of the model are these :math:`k`
weight matrices and :math:`k` bias vectors (there are no weights or biases for
the input nodes in the graph).
'''

import climate
import fnmatch
import gzip
import hashlib
import numpy as np
import pickle
import theano
import theano.tensor as TT
import pdb

from . import layers

logging = climate.get_logger(__name__)

class Network(object):
    '''The network class encapsulates a network computation graph.

    In addition to defining standard functionality for common types of
    feedforward nets, there are also many options for specifying topology and
    regularization, several of which must be provided to the constructor at
    initialization time.

    Parameters
    ----------
    layers : sequence of int, tuple, dict, or :class:`Layer <layers.Layer>`
        A sequence of values specifying the layer configuration for the network.
        For more information, please see :ref:`creating-specifying-layers`.

    weighted : bool, optional
        If True, the network will require an additional input during training
        that provides weights for the target outputs of the network; the weights
        will be the last input argument to the network, and they must be the
        same shape as the target output.

        This can be particularly useful for recurrent networks, where the length
        of each input sequence in a minibatch is not necessarily the same number
        of time steps, or for classifier networks where the prior proabibility
        of one class is significantly different than another. The default is not
        to use weighted outputs.

    sparse_input : bool
        If True, create an input variable that can hold a sparse matrix.
        Defaults to False, which assumes all arrays are dense.

    Attributes
    ----------
    inputs : list of theano variables
        A list of the symbolic inputs this network expects during training.
    layers : list of :class:`Layer <layers.Layer>`
        A list of the layers in this network model.
    weighted : bool
        True iff this network expects target weight inputs during training.
    '''

    def __init__(self, layers, weighted=False, encdec = False,sparse_input=False):
        self._graphs = {}     # cache of symbolic computation graphs
        self._functions = {}  # cache of callable feedforward functions
        self.weighted = weighted
        self.encdec = encdec
        self.inputs = list(self._setup_vars(sparse_input))
        self.layers = []
        for i, layer in enumerate(layers):
            self.add_layer(layer, is_output=i == len(layers) - 1)
        logging.info('network has %d total parameters', self.num_params)

    def _setup_vars(self, sparse_input):
        '''Setup Theano variables required by our network.

        Subclasses must implement this method to specify variables that are
        required for training the model. For example, a supervised model might
        specify a variable that represents the target output for a particular
        input.

        Parameters
        ----------
        sparse_input : bool
            If True, create an input variable that can hold a sparse matrix.
            Defaults to False, which assumes all arrays are dense.

        Returns
        -------
        vars : sequence of theano variables
            The variables that this network requires as inputs during training.
        '''
        raise NotImplementedError

    def add_layer(self, layer, is_output=False):
        '''Add a layer to our network graph.

        Parameters
        ----------
        layer : int, tuple, dict, or :class:`Layer <layers.Layer>`
            A value specifying the layer to add. For more information, please
            see :ref:`creating-specifying-layers`.
        is_output : bool, optional
            True iff this is the output layer for the graph. This influences the
            default activation function used for the layer: output layers in
            most models have a linear activation, while output layers in
            classifier networks default to a softmax activation.
        '''
        # if the given layer is a Layer instance, just add it and move on.
        if isinstance(layer, layers.Layer):
            self.layers.append(layer)
            return

        # for the first layer, create an 'input' layer.
        if len(self.layers) == 0:
            if isinstance(layer, dict):
                src_size = layer['src_size']
                dst_size = layer['dst_size']
                emb_size = layer['emb_size']
                self.layers.append(layers.build('inputencdec', src_size = src_size,
                    dst_size = dst_size,
                    emb_size = emb_size,
                    name = 'in'))
            else:
                assert isinstance(layer, int), 'first layer must be an int'
                self.layers.append(layers.build('input', size=layer, name='in'))
            return

        # here we set up some defaults for constructing a new layer.
        act = getattr(self, 'DEFAULT_OUTPUT_ACTIVATION', 'linear')
        form = 'feedforward'
        kwargs = dict(
            name='out' if is_output else 'hid{}'.format(len(self.layers)),
            activation=act if is_output else 'logistic',
            inputs={self.layers[-1].output_name(): self.layers[-1].size},
            size=layer,
        )

        # if layer is a tuple, assume that it contains one or more of the following:
        # - a layers.Layer subclass to construct (type)
        # - the name of a layers.Layer class (str)
        # - the name of an activation function (str)
        # - the number of units in the layer (int)
        if isinstance(layer, (tuple, list)):
            for el in layer:
                try:
                    if issubclass(el, layers.Layer):
                        form = el.__name__
                except TypeError:
                    pass
                if isinstance(el, str):
                    if layers.Layer.is_registered(el):
                        form = el
                    else:
                        kwargs['activation'] = el
                if isinstance(el, int):
                    kwargs['size'] = el

        # if layer is a dictionary, try to extract a form for the layer, and
        # override our default keyword arguments with the rest.
        if isinstance(layer, dict):
            layer = dict(layer)
            if 'form' in layer:
                form = layer.pop('form').lower()
            kwargs.update(layer)

        if isinstance(form, str) and form.lower() == 'bidirectional':
            if not (isinstance(layer, dict) and 'name' in layer):
                kwargs['name'] = 'bd{}{}'.format(
                    kwargs.get('worker', 'rnn'), len(self.layers))

        if isinstance(form, str) and form.lower() == 'tied':
            partner = kwargs.get('partner')
            if isinstance(partner, str):
                # if the partner is named, just get that layer.
                partner = [l for l in self.layers if l.name == partner][0]
            else:
                # otherwise, we look backwards through our list of layers.
                # any "tied" layer that we find increases a counter by one,
                # and any "untied" layer decreases the counter by one. our
                # partner is the first layer we find with count zero.
                #
                # this is intended to handle the hopefully common case of a
                # (possibly deep) tied-weights autoencoder.
                tied = 1
                partner = None
                for l in self.layers[::-1]:
                    tied += 1 if isinstance(l, layers.Tied) else -1
                    if tied == 0:
                        partner = l
                        break
                assert partner is not None, \
                    'could not find tied layer partner for {} in {}'.format(
                        layer, self.layers)
            kwargs['partner'] = partner

        self.layers.append(layers.build(form, **kwargs))

    def error(self, outputs):
        '''Build a theano expression for computing the network error.

        Parameters
        ----------
        outputs : dict mapping str to theano expression
            A dictionary of all outputs generated by the layers in this network.

        Returns
        -------
        error : theano expression
            A theano expression representing the network error.
        '''
        err = outputs[self.output_name()] - self.x
        if self.weighted:
            return (self.weights * err * err).sum() / self.weights.sum()
        return (err * err).mean()

    def _hash(self, **kwargs):
        '''Construct a string key for representing a computation graph.

        This key will be unique for a given network topology and set of keyword
        arguments.

        Returns
        -------
        key : str
            A hash representing the computation graph for the current network.
        '''
        def add(s):
            h.update(str(s).encode('utf-8'))
        h = hashlib.md5()
        add(kwargs)
        for l in self.layers:
            add('{}{}{}'.format(l.__class__.__name__, l.name, l.size))
        return h.hexdigest()

    def build_graph(self, **kwargs):
        '''Connect the layers in this network to form a computation graph.

        Parameters
        ----------
        noise : dict mapping str to float, optional
            A dictionary that maps layer output names to standard deviation
            values. For an output "layer:output" in the graph, white noise with
            the given standard deviation will be added to the output. Defaults
            to 0 for all layer outputs.
        dropout : dict mapping str to float in [0, 1], optional
            A dictionary that maps layer output names to dropout values. For an
            output "layer:output" in the graph, the given fraction of units in
            the output will be randomly set to 0. Default to 0 for all layer
            outputs.

        Returns
        -------
        outputs : list of theano variables
            A list of expressions giving the output of each layer in the graph.
        updates : list of update tuples
            A list of updates that should be performed by a theano function that
            computes something using this graph.
        '''
        key = self._hash(**kwargs)
        if key not in self._graphs:
            noise = kwargs.get('noise')
            if noise is None:
                noise = {}
                for i, l in enumerate(self.layers):
                    which = 'hidden_noise'
                    if i == 0:
                        which = 'input_noise'
                    if i == len(self.layers) - 1:
                        which = 'output_noise'
                    noise[l.output_name()] = kwargs.get(which, 0)
            dropout = kwargs.get('dropout')
            if dropout is None:
                dropout = {}
                for i, l in enumerate(self.layers):
                    which = 'hidden_dropouts'
                    if i == 0:
                        which = 'input_dropouts'
                    if i == len(self.layers) - 1:
                        which = 'output_dropouts'
                    dropout[l.output_name()] = kwargs.get(which, 0)
            if self.encdec:
                outputs, updates = dict(src=self.src, src_mask = self.src_mask,
                        dst = self.dst), []
            else:
                outputs, updates = dict(x=self.x), []
            for i, layer in enumerate(self.layers):
                out, upd = layer.connect(outputs, noise, dropout)
                outputs.update(out)
                updates.extend(upd)
                outputs['out'] = outputs[layer.output_name()]
            self._graphs[key] = outputs, updates
        return self._graphs[key]

    @property
    def params(self):
        '''A list of the learnable theano parameters for this network.'''
        return [p for l in self.layers for p in l.params]

    @property
    def num_params(self):
        '''Number of parameters in the entire network model.'''
        return sum(l.num_params for l in self.layers)

    def output_name(self):
        '''Get the fully-scoped name of the output for this network model.

        Returns
        -------
        name : str
            A name of the form "{layer}:{output}" that indicates the overall
            output for the network.
        '''
        return self.layers[-1].output_name()

    def find(self, layer, param):
        '''Get a parameter from a layer in the network.

        Parameters
        ----------
        layer : int or str
            The layer that owns the parameter to return.

            If this is an integer, then 0 refers to the input layer, 1 refers
            to the first hidden layer, 2 to the second, and so on.

            If this is a string, the layer with the corresponding name, if any,
            will be used.

        param : int or str
            Name of the parameter to retrieve from the specified layer, or its
            index in the parameter list of the layer.

        Raises
        ------
        KeyError
            If there is no such layer, or if there is no such parameter in the
            specified layer.

        Returns
        -------
        param : theano shared variable
            A shared parameter variable from the indicated layer.
        '''
        for i, l in enumerate(self.layers):
            if layer == i or layer == l.name:
                return l.find(param)
        raise KeyError(layer)

    def feed_forward(self, x, **kwargs):
        '''Compute a forward pass of all layers from the given input.

        All keyword arguments are passed directly to :func:`build_graph`.

        Parameters
        ----------
        x : ndarray (num-examples, num-variables)
            An array containing data to be fed into the network. Multiple
            examples are arranged as rows in this array, with columns containing
            the variables for each example.

        Returns
        -------
        layers : list of ndarray (num-examples, num-units)
            The activation values of each layer in the the network when given
            input `x`. For each of the hidden layers, an array is returned
            containing one row per input example; the columns of each array
            correspond to units in the respective layer. The "output" of the
            network is the last element of this list.
        '''
        key = self._hash(**kwargs)
        if key not in self._functions:
            outputs, updates = self.build_graph(**kwargs)
            labels, exprs = list(outputs.keys()), list(outputs.values())
            self._functions[key] = (
                labels,
                theano.function([self.x], exprs, updates=updates),
            )
        labels, f = self._functions[key]
        return dict(zip(labels, f(x)))

    def predict(self, x):
        '''Compute a forward pass of the inputs, returning the network output.

        Parameters
        ----------
        x : ndarray (num-examples, num-variables)
            An array containing data to be fed into the network. Multiple
            examples are arranged as rows in this array, with columns containing
            the variables for each example.

        Returns
        -------
        y : ndarray (num-examples, num-variables)
            Returns the values of the network output units when given input `x`.
            Rows in this array correspond to examples, and columns to output
            variables.
        '''
        return self.feed_forward(x)[self.output_name()]

    def score(self, x, y, w=None):
        '''Compute R^2 coefficient of determination for a given labeled input.

        Parameters
        ----------
        x : ndarray (num-examples, num-inputs)
            An array containing data to be fed into the network. Multiple
            examples are arranged as rows in this array, with columns containing
            the variables for each example.
        y : ndarray (num-examples, num-outputs)
            An array containing expected target data for the network. Multiple
            examples are arranged as rows in this array, with columns containing
            the variables for each example.

        Returns
        -------
        r2 : float
            The R^2 correlation between the prediction of this netork and its
            target output.
        '''
        u = y - self.predict(x)
        v = y - y.mean()
        if w is None:
            w = np.ones_like(u)
        return 1 - (w * u * u).sum() / (w * v * v).sum()

    def save(self, filename):
        '''Save the state of this network to a pickle file on disk.

        Parameters
        ----------
        filename : str
            Save the state of this network to a pickle file at the named path.
            If this name ends in ".gz" then the output will automatically be
            gzipped; otherwise the output will be a "raw" pickle.
        '''
        opener = gzip.open if filename.lower().endswith('.gz') else open
        handle = opener(filename, 'wb')
        pickle.dump(self, handle, -1)
        handle.close()
        logging.info('%s: saved model', filename)

    @classmethod
    def load(cls, filename):
        '''Load a saved network from disk.

        Parameters
        ----------
        filename : str
            Load the state of a network from a pickle file at the named path. If
            this name ends in ".gz" then the input will automatically be
            gunzipped; otherwise the input will be treated as a "raw" pickle.
        '''
        opener = gzip.open if filename.lower().endswith('.gz') else open
        handle = opener(filename, 'rb')
        model = pickle.load(handle)
        handle.close()
        logging.info('%s: loaded model', filename)
        return model

    def loss(self, **kwargs):
        '''Return a variable representing the loss for this network.

        The loss includes both the error for the network as well as any
        regularizers that are in place.

        Parameters
        ----------
        weight_l1 : float, optional
            Regularize the L1 norm of unit connection weights by this constant.
        weight_l2 : float, optional
            Regularize the L2 norm of unit connection weights by this constant.
        hidden_l1 : float, optional
            Regularize the L1 norm of hidden unit activations by this constant.
        hidden_l2 : float, optional
            Regularize the L2 norm of hidden unit activations by this constant.
        contractive : float, optional
            Regularize model using the Frobenius norm of the hidden Jacobian.
        noise : float, optional
            Standard deviation of desired noise to inject into input.
        dropout : float in [0, 1], optional
            Proportion of input units to randomly set to 0.

        Returns
        -------
        loss : theano expression
            A theano expression representing the loss of this network.
        '''
        outputs, _ = self.build_graph(**kwargs)
        hiddens = [outputs[l.output_name()] for l in self.layers[1:-1]]
        regularizers = dict(
            weight_l1=(abs(w).mean() for l in self.layers
                       for w in l.params if w.ndim > 1),
            weight_l2=((w * w).mean() for l in self.layers
                       for w in l.params if w.ndim > 1),
            hidden_l1=(abs(h).mean() for h in hiddens),
            hidden_l2=((h * h).mean() for h in hiddens),
            contractive=(TT.sqr(TT.grad(h.mean(), self.x)).mean()
                         for h in hiddens),
        )
        return self.error(outputs) + sum(
            kwargs[weight] * sum(expr)
            for weight, expr in regularizers.items()
            if kwargs.get(weight, 0) > 0)

    def monitors(self, **kwargs):
        '''Return expressions that should be computed to monitor training.

        Returns
        -------
        monitors : list of (name, expression) pairs
            A list of named monitor expressions to compute for this network.
        '''
        outputs, _ = self.build_graph(**kwargs)
        monitors = [('err', self.error(outputs))]

        def parse_pattern(pattern):
            '''Yield graph expressions that match the given pattern.'''
            for name, expr in outputs.items():
                if fnmatch.fnmatch(name, pattern):
                    yield name, expr
            for l in self.layers:
                for p in l.params:
                    if fnmatch.fnmatch(p.name, pattern):
                        yield p.name, p

        def parse_levels(levels):
            '''Yield named monitor callables.'''
            if isinstance(levels, dict):
                levels = levels.items()
            if isinstance(levels, (int, float)):
                levels = [levels]
            for level in levels:
                if isinstance(level, (tuple, list)):
                    label, call = level
                    yield ':{}'.format(label), call
                if isinstance(level, (int, float)):
                    def call(expr):
                        return (expr < level).mean()
                    yield '<{}'.format(level), call

        inputs = kwargs.get('monitors', {})
        if isinstance(inputs, dict):
            inputs = inputs.items()
        for pattern, levels in inputs:
            for name, expr in parse_pattern(pattern):
                for key, value in parse_levels(levels):
                    monitors.append(('{}{}'.format(name, key), value(expr)))

        return monitors

    def updates(self, **kwargs):
        '''Return expressions to run as updates during network training.

        Returns
        -------
        updates : list of (parameter, expression) pairs
            A list of named parameter update expressions for this network.
        '''
        _, updates = self.build_graph(**kwargs)
        return updates

ext: py

reason: import numpy,import theano


text:
import pandas as pd
import numpy as np
import anomaly_detection as ad
import pystan
import fbprophet

def test_pystan():
    """Check pystan for build errors"""
    import pystan

    model_code = "parameters {real y;} model {y ~ normal(0,1);}"
    model = pystan.StanModel(model_code=model_code)  # this will take a minute
    y = model.sampling(n_jobs=1).extract()["y"]
    assert np.round(y.mean(), 0) == 0

def test__format_dataframe():
    prophet_example = ad.utilities._get_prophet_example()
    example = ad.utilities._get_test_example()

    # should pass
    detector = ad.AnomalyDetector(data=prophet_example)
    detector = ad.AnomalyDetector(
        data=example, datetime_column="datetime", target="sales_float"
    )
    detector = ad.AnomalyDetector(
        data=example.set_index("datetime"), target="sales_float"
    )

    # should fail
    try:
        detector = ad.AnomalyDetector(data=example)
    except AssertionError:
        pass

    try:
        detector = ad.AnomalyDetector(data=example, target="sales_float")
    except AssertionError:
        pass

    assert set(("sales_float", "datetime")).issubset(set(example))
    assert set(("y", "ds")).issubset(set(detector.data.columns))

def test_fit():
    detector = ad.utilities._get_detector_example()
    detector.fit()

    assert detector.model

def test_predict():
    detector = ad.utilities._get_detector_example()
    detector.predict()

    forecasts = detector.results

    assert (forecasts["yhat"] >= 100000).all() & (forecasts["yhat"] <= 518253).all()

def test_detect_anomalies():
    detector = ad.utilities._get_detector_example()
    detector.detect_anomalies()

    total_rows = len(detector.results)

    null_count = sum(detector.results["anomaly_score"] == 0)

    # As a rough guideline, there shouldn't be that many anomalies
    assert (null_count / total_rows) >= 0.92

def test_get_results():
    example = ad.utilities._get_test_example()
    prophet_example = ad.utilities._get_prophet_example()

    detector = ad.utilities._get_detector_example(type_="test")
    prophet_detector = ad.utilities._get_detector_example(type_="prophet")

    for data, detect in [(example, detector), (prophet_example, prophet_detector)]:

        detect.predict()

        detect.detect_anomalies()

        final_output = detect.get_results()

        # check that changepoints are assigned correctly
        output_changepoints = set(
            final_output.loc[
                final_output["changepoint_flag"] == 1, detect.datetime_column
            ]
        )
        model_changepoints = set(detect.model.changepoints)
        print(detect.results)
        print(final_output)
        assert output_changepoints == model_changepoints

        # check that column / date order hasn't changed
        final_output_set = set(final_output[[detect.target, detect.datetime_column]])
        data_set = set(data[[detect.target, detect.datetime_column]])
        assert final_output_set == data_set

        # check that correct columns are included
        assert set(("anomaly_score", "changepoint_flag")).issubset(
            set(final_output.columns)
        )

def test_plot_forecasts():
    detector = ad.utilities._get_detector_example()
    results = detector.detect_anomalies()

    print(detector.plot_forecasts())

def test_plot_components():
    detector = ad.utilities._get_detector_example()
    results = detector.detect_anomalies()

    print(detector.plot_components())

def test_plot_anomalies():
    detector = ad.utilities._get_detector_example()
    results = detector.detect_anomalies()

    print(detector.plot_anomalies())

if __name__ == "__main__":
    test_pystan()
    test__format_dataframe()
    test_fit()
    test_predict()
    test_detect_anomalies()
    test_get_results()
    test_plot_forecasts()
    test_plot_components()
    test_plot_anomalies()
    print("Passed!")

max_stars_repo_name: ntlind/anomaly-detection

max_issues_repo_name: ntlind/anomaly-detection

max_forks_repo_name: ntlind/anomaly-detection

reason: import numpy


text:
'''
ECE276A WI20 HW1
Stop Sign Detector
'''

import os, cv2
import numpy as np
import skimage

class StopSignDetector():
    def __init__(self):
        '''
            Initilize your stop sign detector with the attributes you need,
            e.g., parameters of your classifier
        '''
        self.weights = np.array([-1198.70681915, 2975.83862542, -19176.06813495, 9368.96722851]) # removed a third of the red pixels so the entire dataset would be evenly split 50/50

        self.ARC_COEFF = 0.01
        self.MIN_PTS_IN_CONTOUR = 8
        self.MAX_PTS_IN_CONTOUR = 12
        self.MIN_AREA_RATIO = 0.00015
        self.MAX_ECCENTRICITY = 0.61
        self.MAX_SIGNS_PER_IMG = 2
        self.KERNEL_SIZE = 3

    def segment_image(self, img):
        '''
            Obtain a segmented image using a color classifier,
            Logistic Regression

            Inputs:
                img - original image
            Outputs:
                mask_img - a binary image with 1 if the pixel in the original image is red and 0 otherwise
        '''

        x = img.flatten().reshape(img.shape[0] * img.shape[1], img.shape[2])

        intercept = np.ones((x.shape[0], 1))
        x = np.concatenate((intercept, x), axis = 1) # add bias column

        mask_img = np.matmul(x, self.weights)
        mask_img = mask_img.reshape(img.shape[0], img.shape[1])
        mask_img = 1.0 * (mask_img > 0) # binarize image
        mask_img = skimage.img_as_ubyte(mask_img) # convert to uint8

        return mask_img

    def get_bounding_box(self, img):
        '''
            Find the bounding box of the stop sign
            call other functions in this class if needed

            Inputs:
                img - original image
            Outputs:
                boxes - a list of lists of bounding boxes. Each nested list is a bounding box in the form of [x1, y1, x2, y2]
                where (x1, y1) and (x2, y2) are the top left and bottom right coordinate respectively. The order of bounding boxes in the list
                is from left to right in the image.

            Our solution uses xy-coordinate instead of rc-coordinate. More information: http://scikit-image.org/docs/dev/user_guide/numpy_images.html#coordinate-conventions
        '''
        boxes = []
        mask_img = self.segment_image(img)
        img_area = img.shape[0] * img.shape[1]

        # get thresholded image and blur. Canny, Sobel, Laplacian all had worse performance
        _, thresh_img = cv2.threshold(mask_img, 0, 255, cv2.THRESH_BINARY)

        # Gaussian blur and erosion for connected, smooth contours
        kernel = (self.KERNEL_SIZE, self.KERNEL_SIZE)
        smoothed_img = cv2.GaussianBlur(thresh_img, kernel, cv2.BORDER_DEFAULT)
        eroded_img = cv2.erode(smoothed_img, kernel, iterations = 1)

        # find largest contours
        contours, _ = cv2.findContours(eroded_img, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)[-2:]
        contours = sorted(contours, key = cv2.contourArea, reverse = True)[:self.MAX_SIGNS_PER_IMG]

        for contour in contours:
            perimeter = self.ARC_COEFF * cv2.arcLength(contour, True)
            poly_approx = cv2.approxPolyDP(contour, perimeter, True)

            if (len(poly_approx) >= self.MIN_PTS_IN_CONTOUR):
                _, axes, _ = cv2.fitEllipse(poly_approx)

                eccentricity = np.sqrt(1 - np.square(min(axes) / max(axes)))

                if (len(poly_approx) <= self.MAX_PTS_IN_CONTOUR or eccentricity <= self.MAX_ECCENTRICITY):
                    (x, y, w, h) = cv2.boundingRect(poly_approx)
                    area = w * h

                    if ((area / img_area) >= self.MIN_AREA_RATIO):
                        bounds = [x, img.shape[0] - y - h, x + w, img.shape[0] - y] # convert coordinates to project specification
                        boxes.append(bounds)

#            out = cv2.drawContours(img, [contour], -1, (0, 255, 0), 3)
        print(boxes)
        return boxes

if __name__ == '__main__':
    folder = "trainset"
    my_detector = StopSignDetector()
    for filename in os.listdir(folder):
        # read one test image
        img = cv2.imread(os.path.join(folder, filename))
        cv2.imshow('image', img)
        cv2.waitKey(0)
        cv2.destroyAllWindows()

    # Display results:
    # (1) Segmented images
    #	 mask_img = my_detector.segment_image(img)
    # (2) Stop sign bounding box
    #    boxes = my_detector.get_bounding_box(img)
    # The autograder checks your answers to the functions segment_image() and get_bounding_box()
    # Make sure your code runs as expected on the testset before submitting to Gradescope


text:
# This is the Method 1 for MaxPooling with NumPy

# Import the necessary packages
import numpy as np

#generating a 4x4 matrix with random elements ranging from 0 to 20
large_matrix = np.random.randint(0,20,(4,4))
print(large_matrix)

#This functioin is defined to split the 4x4 matrix to 2x2 smaller matrices
def split(matrix, n_rows, n_cols):
  r, c = matrix.shape
  return (matrix.reshape(c//n_rows, n_rows, -1, n_cols).swapaxes(1,2).reshape(-1,n_rows, n_cols))

#calling the function and storing the matrices to A, B, C, and D
A, B, C, D = split(large_matrix,2,2)
print("\n")
print("A=",A)
print("B=", B)
print("C=", C)
print("D=", D)

#Retrieving the highest value elements from the generated matrices
final = np.array([A.max(), B.max(),C.max(), D.max()])

#generating the matrix with the maximum values
final_matrix = final.reshape(2,2)
print("final_matrix= \n ",final_matrix)

max_stars_repo_name: the-engrossed/Maxpool-with-NumPy

max_issues_repo_name: the-engrossed/Maxpool-with-NumPy

max_forks_repo_name: the-engrossed/Maxpool-with-NumPy


text:
#%%

# import modules
import pandas as pd
import uuid
import re
import math
import networkx as nx
import json

#%%

class Bibliography:
    """
    a bibliography containing bibliographic data
    """

    def __init__(self, path_list, db_src):
        self.bib_df = _read_bib(path_list, db_src)
        self.bib_dict = _convert_bib_to_dict(self.bib_df, db_src)
        self.cit_rel = _extract_cit_rel(self.bib_dict)

    def export_bib_dict(self, path):
        """
        write bibliographic dictionary to a json file
        """
        with open(path, 'w') as json_file:
            json.dump(self.bib_dict, json_file)

#%%

class CitationNetwork:
    """
    a citation network in the form of a weighted DAG
    """

    def __init__(self, BibObject, keywords, w_method):
        self.keywords = keywords
        self.cit_net = _create_citation_network(BibObject, keywords, w_method)
        self.out_deg_dist = _extract_out_deg_dist(self.cit_net)

    def add_enclosing(self):
        """
        add sink and source to citation network
        """
        self.cit_net.add_edges_from(_enclose_cit_rel(self.cit_net))

#%%

class MainPathNetwork:
    """
    a network of main paths constructed from a citation network
    """

    def __init__(self, CitNetObject, mc=1, iterations=1):
        self.cit_net_raw = CitNetObject
        self.main_paths, self.main_net = _main_path_analysis(CitNetObject, mc, iterations)
        self.main_path_net = nx.DiGraph(kw=CitNetObject.keywords)

    def create_network(self, mw=1):
        """
        create a main path network
        contains all edges that are part of a main path
        network is pruned with some options
        """
        self.main_path_net.add_weighted_edges_from(self.main_net)
        _rm_list_e = []
        for e in self.main_path_net.edges:
            if self.main_path_net.get_edge_data(e[0], e[1])['weight'] < mw:
                _rm_list_e.append(e)
        self.main_path_net.remove_edges_from(_rm_list_e)
        _rm_list_n = []
        for n in self.main_path_net.nodes:
            if (self.main_path_net.in_degree(n) == 0 and self.main_path_net.out_degree(n) == 0):
                _rm_list_n.append(n)
        self.main_path_net.remove_nodes_from(_rm_list_n)

#%%

def _read_bib(path_list, db_src):
    """
    read output data from bibliographic database
    args:
        path_list: a list of paths to read from
        db_src: database source of the bibliographic data
    returns:
        bib_df: a bibliographic dataframe
    """

    # get database specifics
    db_specs = _switch_case_db(db_src)

    # read file(s)
    if len(path_list) == 1:
        bib_df = pd.read_csv(
            path_list[0], usecols=db_specs['columns'],
            dtype=db_specs['dtype'], sep=db_specs['sep'],
            index_col=False
        )
    else:
        bib_df = pd.concat(
            [
                pd.read_csv(
                    path, usecols=db_specs['columns'],
                    dtype=db_specs['dtype'], sep=db_specs['sep'],
                    index_col=False
                ) for path in path_list
            ],
            ignore_index=True
        )

    # some more formatting
    bib_df['db_src'] = db_src
    bib_df = bib_df.rename(columns=db_specs['new_cols'])

    return bib_df

#%%

def _convert_bib_to_dict(bib_df, db_src):
    """
    convert bibliographic dataframe into python dictionary
    args:
        bib_df: a bibliographic dataframe
        db_src: database source of the bibliographic data
    returns:
        bib_dict: a bibliographic dictionary
    """

    # get database specifics
    db_specs = _switch_case_db(db_src)

    # extract and reformat data
    keys = bib_df.columns
    bib_dict = {}

    # create a dictionary entries
    for j in range(bib_df.shape[0]):
        entry = {}
        for key in keys:
            entry[key] = bib_df.loc[j, key]
        entry = _split_columns(entry, db_specs['to_split'])  # split non-atomic columns
        bib_dict[str(uuid.uuid4())] = entry

    return bib_dict

#%%

# switch scopus specifics
def _switch_scopus():
    """
    in case db_src='scopus', use these specifications
    args:
        None
    returns:
        wos: the scopus specifics
    """

    # define scopus specifics
    scopus = {
        'columns': [
            'Title', 'Authors', 'Author(s) ID', 'Source title',
            'Year', 'Volume', 'Issue', 'Cited by',
            'References', 'DOI', 'ISBN', 'ISSN',
            'CODEN', 'PubMed ID', 'EID'
        ],
        'dtype': {
            'Title': 'str', 'Authors': 'str',
            'Author(s) ID': 'str', 'Source title': 'str',
            'Year': 'float64', 'Volume': 'str',
            'Issue': 'str', 'Cited by': 'float64',
            'References': 'str', 'DOI': 'str',
            'ISBN': 'str', 'ISSN': 'str',
            'CODEN': 'str', 'PubMed ID': 'str',
            'EID': 'str'
        },
        'sep': ',',
        'new_cols': {
            'Title': 'title', 'Authors': 'authors',
            'Author(s) ID': 'scopus_author_id', 'Source title': 'source',
            'Year': 'year', 'Volume': 'vol',
            'Issue': 'issue', 'Cited by': 'cit_count',
            'References': 'refs', 'DOI': 'doi',
            'ISBN': 'isbn', 'ISSN': 'issn',
            'CODEN': 'coden', 'PubMed ID': 'pmid',
            'EID': 'eid'
        },
        'to_split': {
            'authors': ',', 'scopus_author_id': ';',
            'refs': ';', 'isbn': ';'
        }
    }

    return scopus

#%%

# switch web of science specifics
def _switch_wos():
    """
    in case db_src='wos', use these specifications
    args:
        None
    returns:
        wos: the wos specifics
    """

    # define web of science specifics
    wos = {
        'columns': [
            'TI', 'AU', 'AF', 'RI', 'OI', 'SO', 'PY', 'VL',
            'IS', 'TC', 'Z9', 'CR', 'U1', 'U2', 'DI', 'D2',
            'BN', 'SN', 'PM', 'J9', 'JI', 'BP', 'EP'
        ],
        'dtype': {
            'TI': 'str', 'AU': 'str', 'AF': 'str', 'RI': 'str',
            'OI': 'str', 'SO': 'str', 'PY': 'float64', 'VL': 'str',
            'IS': 'str', 'TC': 'float64', 'Z9': 'float64', 'CR': 'str',
            'U1': 'float64', 'U2': 'float64', 'DI': 'str', 'D2': 'str',
            'BN': 'str', 'SN': 'str', 'PM': 'str', 'J9': 'str',
            'JI': 'str', 'BP': 'str', 'EP': 'str'
        },
        'sep': '\t',
        'new_cols': {
            'TI': 'title', 'AU': 'authors', 'AF': 'authors_full',
            'RI': 'researcher_id', 'OI': 'orcid', 'SO': 'source',
            'PY': 'year', 'VL': 'vol', 'IS': 'issue',
            'TC': 'cit_count', 'Z9': 'cit_count_z9', 'CR': 'refs',
            'U1': 'usage_count_180d', 'U2': 'usage_count_2013', 'DI': 'doi',
            'D2': 'book_doi', 'BN': 'isbn', 'SN': 'issn',
            'PM': 'pmid', 'J9': 'src_abb_29', 'JI': 'src_abb_iso',
            'BP': 'start_page', 'EP': 'end_page'
        },
        'to_split': {
            'authors': ';', 'authors_full': ';', 'researcher_id': ';',
            'orcid': ';', 'refs': ';', 'isbn': ';'
        }
    }

    return wos

#%%

# switch case replacement
# adapted from: https://www.pydanny.com/why-doesnt-python-have-switch-case.html
def _switch_case_db(arg):
    """
    replacement for switch case
    args:
        arg: the case to execute
    returns:
        func(): the switch function
    """

    # dictionary referring to switch statements
    switcher = {
        'scopus': _switch_scopus,
        'wos': _switch_wos
    }

    # print status if exception occurs
    if arg not in ['scopus', 'wos']:
        print(
            'Expection: Unknown db_source ' +
            str(arg) +
            '\nOutput might not be in proper format'
        )

    # switch function
    func = switcher.get(arg, lambda: [[], {}])

    return func()

#%%

# split non-atomic columns
def _split_columns(entry, split_list):
    """
    split pre-defined dictionary entries along pre-defined separators
    args:
        entry: a bibliography entry
        split_list: a pre-defined list of columns to separate
    returns:
        entry: the split entry
    """

    # function to strip trailing spaces
    space_stripper = lambda x: x.rstrip().lstrip()

    # iterate over pre-defined split_list
    for label, separator in split_list.items():
        try:
            entry[label] = space_stripper(entry[label]).rstrip(';')
            entry[label] = entry[label].split(separator)
            entry[label] = list(
                map(space_stripper, entry[label])
            )
        except:
            continue

    return entry

#%%

def _extract_cit_rel(bib_dict):
    """
    extract the citation relation
    args:
        bib_dict: a bibliographic dictionary
    returns:
        cit_rel: a list of tuples (x, y) with 'x cites y'
    """

    # initialize citation relation
    cit_rel = []

    # iterate over all bibliography entries
    for key in bib_dict.keys():
        doi_at_key = bib_dict[key]['doi']
        # check if a doi is available, if not: continue
        if len(str(doi_at_key)) > 8:
            pass
        else:
            continue
        refs_at_key = bib_dict[key]['refs']
        # try to extract doi and append to citation relation
        try:
            for ref_idx in range(len(refs_at_key)):
                ref_doi = _extract_doi(
                    refs_at_key[ref_idx], bib_dict[key]['db_src']
                )
                if ref_doi == 'NO_DOI':
                    continue
                if ref_doi != doi_at_key:
                    cit_rel.append((ref_doi, doi_at_key))
        except:
            continue

    return cit_rel

#%%

# extract doi from reference string
# with a little help from: https://www.crossref.org/blog/dois-and-matching-regular-expressions/
def _extract_doi(ref_elem, db_src='wos'):
    """
    extract doi from reference (CAUTION: only works for db_src == 'wos' so far!)
    args:
        ref_elem: a reference
        db_src: database source of the bibliographic data
    returns:
        doi if doi is extractable
        'NO_DOI' if no doi is extractable
    """

    # currently: only works for web of science
    if db_src == 'wos':

        regex_doi = re.compile('10.\d{4,9}/[-._;()/:A-Za-z0-9]+$')  # define regex for doi

        # try searching for doi in reference string
        try:
            res = regex_doi.search(ref_elem)
            doi = res.group(0)
            return doi
        # if not successful: return value for missing doi
        except:
            return 'NO_DOI'
    # <CODE FOR SCOPUS HERE> -> elif db_src == 'scopus': <CODE>
    else:
        return 'NO_DOI'

#%%

#%%

def _create_citation_network(BibObject, keywords=['citation network'], w_method='nppc'):
    """
    create a citation network
    args:
        BibObject: Bibliography object
        keywords: list of keywords as graph attributes (for later reference)
        w_method: weighting method (so far only 'nppc')
    returns:
        net: a citation network
    """

    # create directed graph
    net = nx.DiGraph(kw=keywords)
    net.add_edges_from(BibObject.cit_rel)
    net = _break_cycles(net)
    net = _compute_edge_weights(net, w_method)

    return net

#%%

def _compute_edge_weights(net, method='nppc'):
    """
    compute the edge weights of a citation network
    args:
        net: a citation network
        method: a weighting method (so far only 'nppc')
    returns:
        net: the net with weighted edges
    """

    if method == 'nppc':  # node pair projection count
        # extract all pairs of connected nodes
        con_nodes = _find_connected_nodes(net)
        # compute nppc weights
        for sub in [nx.all_simple_paths(
                net, source=pair[0], target=pair[1]
            ) for pair in con_nodes]:
            tmp = net.subgraph({node for path in sub for node in path})
            for edge in tmp.edges:
                try:
                    net[edge[0]][edge[1]]['weight'] += 1
                except:
                    net[edge[0]][edge[1]]['weight'] = 1
    else:
        print('This does not seem to be a valid weighting method.')

    return net

#%%

def _find_connected_nodes(net):
    """
    find all connected nodes (a, b) where there is a path from a to b
    args:
        net: a citation network
    returns:
        con_nodes: a list of tuples with all connected nodes
    """

    con_nodes = [
        (struc[0], i)
        for node in net.nodes
        for struc in nx.bfs_successors(net, node)
        for i in struc[1]
    ]

    return(con_nodes)

#%%

def _extract_out_deg_dist(net):
    """
    extract the distribution of out-degrees in the citation network for later pruning
    args:
        net: a citation graph
    returns:
        a pandas Dataframe with the frequency distribution of out-degrees
    """

    out_deg_list = [net.out_degree(node) for node in net.nodes]
    freq = [0] * (max(out_deg_list) + 1)
    for i in out_deg_list:
        freq[i] += 1
    data = {
        'out_degree': list(range(max(out_deg_list) + 1)),
        'freq': freq
    }

    return(pd.DataFrame(data))

#%%

def _enclose_cit_rel(net):
    """
    add sink and source to a citation network
    args:
        net: an unweighted citation network
    return:
        enclosing: list of edges to add sink and source
    """

    try:
        source_edges = [('source', node, {'weight': 1}) for node in net.nodes if net.in_degree(node) == 0]
        sink_edges = [(node, 'sink', {'weight': 1}) for node in net.nodes if net.out_degree(node) == 0]
        return(source_edges + sink_edges)
    except:
        return(False)

#%%

def _break_cycles(net):
    """
    breaks potential cycles in a citation net
    args:
        net: a citation net
    returns:
        net: the acyclic citation net
    """

    flag = True
    counter = 0
    while(flag):
        try:
            counter += 1
            cycles = nx.find_cycle(net)
            net.remove_edge(cycles[-1][0], cycles[-1][1])
        except:
            if counter == 1:
                print('NOTE: {} edge was removed to break cycles'.format(counter))
            else:
                print('NOTE: {} edges were removed to break cycles'.format(counter))
            flag = False
        if counter >= 100:
            print('NOTE: There are oddly many cycles. Please check your data to avoid problems in the further process.')
            flag = False

    return(net)

#%%

def _main_path_analysis(CitNetObject, mc, iterations):
    """
    conduct main path analysis on a citation network
    args:
        CitNetObject: a CitationNetwork object
        mc: minimum citations for start nodes
        iterations: number of iterations to conduct
    returns:
        main_paths: a list of all mined main paths
        main_net: the edges of all main paths in one list
    """

    # setup
    CitNetObject.cit_net.remove_nodes_from(['source', 'sink'])
    print('NOTE: source and sink have been removed from your citation network')
    init_net = CitNetObject.cit_net

    # iterations
    for i in range(iterations):

        # initialize
        start_nodes = _retrieve_start_nodes(init_net, mc)
        main_paths = []
        main_net = []

        for node in start_nodes:

            # loop setup
            counter = 0  # if paths cross lengths of 100, the while loop is interrupted
            flag = True
            error = False
            mp = []
            cur_nodes = [node]

            # find main path
            while(flag):

                # create list containing all lists of out_edges of the current nodes
                candidates = [
                    [
                        (e[0], e[1], init_net.edges[e[0], e[1]]['weight'])
                        for e in init_net.out_edges(cur_node)
                    ] for cur_node in cur_nodes
                ]

                # extract weights from candidate edges
                weights = [[t[2] for t in l] for l in candidates]

                # determine maximum and prepare next iteration
                mw = [max(w) for w in weights]
                max_idx = [
                    [
                        i for i, j in enumerate(weights[m]) if j == mw[m]
                    ] for m in range(len(mw))
                ]

                # update current nodes and extend current main path
                cur_nodes.clear()
                for i, mi in enumerate(max_idx):
                    next_edges = [candidates[i][j] for j in mi]
                    mp.extend(next_edges)
                    cur_nodes.extend([e[1] for e in next_edges])
                cur_nodes = list(dict.fromkeys(cur_nodes))

                # remove node from current nodes if end of path is reached
                rm_idx = []
                for i in range(len(cur_nodes)):
                    if init_net.out_degree(cur_nodes[i]) == 0:
                        rm_idx.append(i)
                for idx in sorted(rm_idx, reverse=True):
                    del cur_nodes[idx]

                counter += 1

                # stop criteria
                if not cur_nodes:
                    flag = False
                if counter >= 100:
                    print('This takes oddly long. Something must have gone wrong.')
                    error = True
                    flag = False

            # append extracted main path to main path collection and extend main path network
            mp = list(dict.fromkeys(mp))
            main_paths.append(mp)
            main_net.extend(mp)

            # report potential error
            if error:
                print('An error occurred.')
                break

        init_net = nx.DiGraph()
        init_net.add_weighted_edges_from(main_net)

    return(main_paths, main_net)

#%%

def _retrieve_start_nodes(net, mc):
    """
    retrieve nodes with in_degree == 0 as starting points for main path analysis
    args:
        net: a CitationNetwork object
        mc: minimum citations
    returns:
        list of start nodes in net
    """
    return [node for node in net.nodes if (net.in_degree(node) == 0 and net.out_degree(node) > mc)]

max_stars_repo_name: JohannesNakayama/CitationNetworks

max_issues_repo_name: JohannesNakayama/CitationNetworks

max_forks_repo_name: JohannesNakayama/CitationNetworks

reason: import networkx


text:
#!/usr/bin/env python
import numpy
import math
import ant
import types
import random
from random import shuffle

from multiprocessing import Process, Queue

# Parameters used -
# index - stores the city number
# x - x coordinate for the city i
# y - y coordinate for the city i
class City:
    def __init__(self,i=0, x_cord=0, y_cord=0):
        self.index = i
        self.x = x_cord
        self.y = y_cord

# Parameters used -
# num_cities - number of cities for travel
# initial_pheromone - initial value of pheromone deposited
# alpha - weight parameter for pheromones
# beta - weight parameter for desirability (attractiveness), which is inverse of distance
# pheromone_deposit - amount of pheromones which can be deposited
# evaporation_constant - amount of pheromone which will evaporate after a cycle
class ACO:
    def __init__(self, num_cities, initial_pheromone=1, alpha=1, beta=3,
                pheromone_deposit=1, evaporation_constant=0.6):
        self.cities = []
        self.shortest_paths = []
        self.shortest_paths_lens = []
        self.shortest_path_len = -1
        self.evaporationConst = evaporation_constant
        self.pheromone_deposit = pheromone_deposit
        self.pheromone = numpy.full((num_cities,num_cities),initial_pheromone)
        self.alpha = alpha
        self.beta = beta
        self.attractiveness = numpy.zeros((num_cities,num_cities))
        self.routing_table = numpy.full((num_cities,num_cities),(1.00/(num_cities-1)))

    # Adding cities present in the dataset
    def add_cities(self, city):
        if isinstance(city, list):
            self.cities.extend(city)
        else:
            self.cities.append(city)

    # Use euclidean formula to calculate distance between two cities.
    # Cities class has x and y coordinates
    def euclidean_distance(self,city1,city2):
        return  (math.sqrt(math.pow((city1.x-city2.x),2)+math.pow((city1.y-city2.y),2)))

    # Attractiveness is reciprocal of distance. This function calculates
    # attractiveness between cities.
    def calc_attraction(self):
        city_list=self.cities
        for i in range(len(city_list)):
            for j in range(len(city_list)):
                distance = self.euclidean_distance(city_list[i], city_list[j])
                if distance > 0:
                    self.attractiveness[i][j] = 1/distance
                else:
                    self.attractiveness[i][j] = 0.00

    # a - a single ant
    # pheromone deposit on each path is updated by the ant visiting that path
    # we have deposited the phermone uniformly across the path
    # additionally, evaporation part is removed from the pheromone
    def update_pheromone(self,a):
        for i in range(0,len(a.path)-1):
            try:
                curr_pher = self.pheromone[a.path[i].index][a.path[i+1].index]
                self.pheromone[a.path[i].index][a.path[i+1].index] = curr_pher + self.pheromone_deposit/a.path_length
            except:
                break
        self.pheromone = self.pheromone*(1-self.evaporationConst)

    # get the pheromone value for the path between city i and j
    def get_pheromone(self,i,j):
        return self.pheromone[i][j]

    # Calculation for numerator part of probability distribution for the next city selection
    def city_sum(self, city_cur,city_next):
        return ((math.pow(self.pheromone[city_cur.index][city_next.index], self.alpha))*(math.pow(self.attractiveness[city_cur.index][city_next.index],self.beta)))

    # Calculating probability density for the next city selection,
    # and storing that value in the routing table
    def update_routing_table(self, a):
        denom = 0.0
        for c in a.path:
            temp_cities = list(a.path)
            temp_cities.remove(c)
            for valid in temp_cities:
                denom += self.city_sum(c, valid)

            for valid in temp_cities:
                numerator = self.city_sum(c, valid)
                if denom > 0:
                    self.routing_table[c.index][valid.index] = numerator/denom
                else:
                    self.routing_table[c.index][valid.index] = 0

    def mp_tour(self,ant,q):
        ant.reset_ant(self)
        while(ant.unvisited):
            if random.random()<0:
                next_city = ant.unvisited.pop(random.randint(0,len(ant.unvisited)-1))
                ant.path.append(next_city)
                ant.currCity = next_city

            else:
                for c in ant.unvisited:
                    ant.transition_probs.append(ant.get_transition_prob(self,c))

                ant.transition_probs = ant.transition_probs/sum(ant.transition_probs)
                selection = numpy.random.choice(ant.unvisited, 1, p = ant.transition_probs)
                next_city = selection[0]
                ant.path.append(next_city)
                ant.currCity = next_city
                ant.unvisited.pop(ant.unvisited.index(next_city))
            ant.transition_probs = []
        q.put(ant)

    # Getting the best path by traversing the number of ants num_steps times
    def get_best_path(self, num_ants=2, num_steps=3):
        ants = []
        self.calc_attraction()
        shuffle(self.cities)

        for i in range(0, num_ants):
            ants.append(ant.Ant(i, self))

        for step in range(0, num_steps):
            print("Step: {} of {}".format(step+1, num_steps))
            path_lens = []
            paths = []
            procs = []

            q = Queue()

            for a in ants:
                p = Process(target = self.mp_tour, args = (a, q,))
                procs.append(p)
                p.start()

            for p in procs:
                p.join()

            ants = []
            while q.empty() == False:
                ants.append(q.get())

            for a in ants:
                path_len = a.calc_path_length()
                path_lens.append(path_len)
                paths.append(a.path)
                self.update_pheromone(a)
                self.update_routing_table(a)

            best_path_len = min(path_lens)
            best_path = paths[path_lens.index(best_path_len)]

            print("Step best path: {} Step: {}".format(best_path_len, step+1))

            self.shortest_paths.append(best_path)
            self.shortest_paths_lens.append(best_path_len)

        output_index = self.shortest_paths_lens.index(min(self.shortest_paths_lens))
        output_path = self.shortest_paths[output_index]
        self.shortest_path_len = self.shortest_paths_lens[output_index]
        self.shortest_paths = []
        self.shortest_paths_lens = []

        return output_path


text:
import torch
import numpy as np

class Instance(object):
    def __init__(self, words, label):
        self.words = words     # 词序列
        self.label = label     # 标签值

def load_dataset(path):
    insts = []
    with open(path, 'r', encoding='utf-8') as fin:
        for line in fin:
            lbl, raw_data = line.strip().split('|||')
            words = raw_data.split()
            insts.append(Instance(words, lbl))
    np.random.shuffle(insts)
    return insts

# 封装Batch
class Batch(object):
    def __init__(self, src, target):
        self.src = src          # 数据
        self.target = target    # 标签

def batch_iter(dataset, batch_size, wd_vocab, shuffle=True, device=torch.device('cpu')):
    if shuffle:
        np.random.shuffle(dataset)

    # batch数量
    nb_batch = int(np.ceil(len(dataset) / batch_size))

    for i in range(nb_batch):
        batch_data = dataset[i*batch_size: (i+1)*batch_size]
        if shuffle:
            np.random.shuffle(batch_data)

        yield batch_gen(batch_data, wd_vocab, device)

# batch数据变量化，提取数据和标签
def batch_gen(batch_data, wd_vocab, device=torch.device('cpu')):
    batch_size = len(batch_data)
    # 求最长的序列长度
    max_seq_len = max([len(inst.words) for inst in batch_data])

    wd_idx = torch.zeros((batch_size, max_seq_len), dtype=torch.long).to(device)
    lbl_idx = torch.zeros(batch_size, dtype=torch.long).to(device)
    for i, inst in enumerate(batch_data):
        seq_len = len(inst.words)
        wd_idx[i, :seq_len] = torch.tensor(wd_vocab.extwd2idx(inst.words))
        lbl_idx[i] = torch.tensor(wd_vocab.label2index(inst.label))

    return Batch(wd_idx, lbl_idx)

# class BatchIter(object):
#     def __init__(self, dataset):
#         self._dataset = dataset
#         self.x = None
#         self.y = None
#
#     def __iter__(self):
#         return self
#
#     def __next__(self):
#         pass