text:
import numpy, scipy, scipy.optimize, scipy.interpolate
import os, os.path, time, copy, pylab, operator
from scipy import interp
import matplotlib.pyplot as plt

import numpy as np

def myeval(c):
    if c=='None':
        c=None
    elif c=='nan' or c=='NaN':
        c=numpy.nan
    else:
        temp=c.lstrip('0')
        if (temp=='' or temp=='.') and '0' in c:
            c=0
        else:
            c=eval(temp)
    return c

def get_dict_item_keylist(d, keylist):
    return reduce(lambda d, k: d[k], keylist, d)

def removeoutliers_meanstd(arr, nptsoneside, nsig, gapptsoneside=0): #avrages maximum of 2*nptoneside points and usees distance from mean scaled by std compared to nsig to determine if the value should be replaced by the mean. if gapptsoneside>0, will do this leaving a gap around the point in question and using nptsoneside-gaps points for the mean and std
    if nptsoneside==1 and gapptsoneside==0:
        return removesinglepixoutliers(arr, critratiotoneighbors=nsig)
    nsig=max(nsig, 1.)
    nptsoneside=max(nptsoneside, 2.)
    gapptsoneside=min(gapptsoneside, nptsoneside-2.)
    for gap in range(gapptsoneside+1):
        starti=numpy.uint32([max(i-(nptsoneside-gap), 0) for i in range(len(arr))])
        stopi=numpy.uint32([min(i+(nptsoneside-gap)+1, len(arr)) for i in range(len(arr))])
        #print [numpy.append(arr[i0:i], arr[i+1:i1]) for i, i0, i1 in zip(range(len(arr)), starti, stopi)][8]
        #print [(((numpy.append(arr[i0:i], arr[i+1:i1]).mean()-arr[i]))**2, (numpy.append(arr[i0:i], arr[i+1:i1]).std()*nsig)**2) for i, i0, i1 in zip(range(len(arr)), starti, stopi)][8]
        arr=numpy.array([(((numpy.append(arr[i0:i], arr[i+1:i1]).mean()-arr[i]))**2<(numpy.append(arr[i0:i], arr[i+1:i1]).std()*nsig)**2 and (arr[i],) or (numpy.append(arr[i0:i], arr[i+1:i1]).mean(),))[0] for i, i0, i1 in zip(range(len(arr)), starti, stopi)], dtype=arr.dtype)
    return arr

def CalcArrSS(x, WeightExp=1., TestPts=10):
    p=WeightExp
    i=TestPts
    s0=x[:i].std()/i**p+1
    while x[:i].std()/i**p<s0 and i<len(x):
        s0=x[:i].std()/i**p
        i+=TestPts
    return x[:i].mean()

def xx(x, WeightExp=1., TestPts=10):
    p=WeightExp
    i=TestPts
    s0=x[:i].std()/i**p+1
    while x[:i].std()/i**p<s0 and i<len(x):
        s0=x[:i].std()/i**p
        i+=TestPts
    return i

def concat_extrap_ends(x, npts, polyorder=1, lowside=True, highside=True):
    i=numpy.arange(npts, dtype='float64')
    if lowside:
        ans=scipy.polyfit(-1*(i+1.), x[:npts], polyorder)
        x=numpy.concatenate([scipy.polyval(list(ans), i[::-1]), x])
    if highside:
        ans=scipy.polyfit(-1*(i[::-1]-1.), x[-1*npts:], polyorder)
        x=numpy.concatenate([x, scipy.polyval(list(ans), i)])
    return x

def lininterpbetweenregularpoints(existy, interval):
    existy=numpy.array(existy)
    x=numpy.arange(interval,dtype='float32')/interval
    diff=existy[1:]-existy[:-1]
    o=numpy.outer(diff,x)
    return numpy.concatenate([arr+start for arr,start in zip(o,existy[:-1])]+[existy[-1:]])

def interpwithinarr(existind, existy, order=3, interpplotax=None, interpcols=['k', 'r']):
    if order==1:
        existind=numpy.array(existind)
        diff=existind[1:]-existind[:-1]
        if numpy.all(diff==diff[0]):
            return lininterpbetweenregularpoints(existy, diff[0])
    interind=sorted(list(set(numpy.arange(max(existind)+1))-set(existind)))
    yfull=numpy.zeros(max(existind)+1, existy.dtype)
    yfull[existind]=existy[:]
    yfull[interind]=scipy.interpolate.spline(existind, existy, interind, order=order)
    if not interpplotax is None:
        interpplotax.plot(existind, existy, interpcols[0])
        interpplotax.plot(interind,  yfull[interind], interpcols[1])
    return yfull

def savgolsmooth(x, nptsoneside=7, order = 4, dx=1.0, deriv=0, binprior=0): #based on scipy cookbook. x is 1-d array, window is the number of points used to smooth the data, order is the order of the smoothing polynomial, will return the smoothed "deriv"th derivative of x
    if nptsoneside<=1:
        return x
    if binprior>1:
        origlen=len(x)
        x=numpy.array([x[i*binprior:(i+1)*binprior].mean() for i in range(origlen//binprior)])
        dx*=binprior
    side=numpy.uint16(max(nptsoneside, numpy.ceil(order/2.)))
    s=numpy.r_[2*x[0]-x[side:0:-1],x,2*x[-1]-x[-2:-1*side-2:-1]]
    # a second order polynomal has 3 coefficients
    b = numpy.mat([[k**i for i in range(order+1)] for k in range(-1*side, side+1)])
    m = numpy.linalg.pinv(b).A[deriv] #this gives the dth ? of the base array (.A) of the pseudoinverse of b

    # precompute the offset values for better performance
    offsets = range(-1*side, side+1)
    offset_data = zip(offsets, m)

    smooth_data=[numpy.array([(weight * s[i + offset]) for offset, weight in offset_data]).sum() for i in xrange(side, len(s) - side)]
    smooth_data=numpy.array(smooth_data)/(dx**deriv)

    if binprior>1:
        ia=numpy.arange(binprior, dtype='float32')/binprior
        xr=numpy.concatenate([ia*(b-a)+a for a, b in zip(smooth_data[:-1], smooth_data[1:])])
        xr=numpy.concatenate([(smooth_data[1]-smooth_data[0])*ia[:binprior//2]+smooth_data[0], xr, (smooth_data[-1]-smooth_data[-2])*ia[:binprior//2]+smooth_data[-1]])
        smooth_data=numpy.concatenate([xr, (smooth_data[-1]-smooth_data[-2])*ia[:origlen-len(xr)]+smooth_data[-1]])

    return smooth_data

class fitfcns: #datatuples are x1,x2,...,y
    #.finalparams .sigmas .parnames useful, returns fitfcn(x)
    def genfit(self, fcn, initparams, datatuple, markstr='unspecified', parnames=[], interaction=0,  maxfev=2000, weights=None, optimizerfcn=None):
        self.maxfev=maxfev
        self.performfit=True
        self.initparams=initparams
        self.sigmas=scipy.zeros(len(initparams))
        self.parnames=parnames
        self.finalparams=initparams
        self.error=False
        if weights is None:
            def wts(x):
                return 1.
        elif weights=='parabolic':
            a=(datatuple[0][0]+datatuple[0][-1])/2.0
            b=(datatuple[0][-1]-datatuple[0][0])/2.0
            def wts(x):
                return 1.0+((x-a)/b)**2

        def res1(p, x1, y):
            return (y-fcn(p, x1))*wts(x1)

        def res2(p, x1,x2,y):
            return y-fcn(p, x1, x2)

        def res3(p, x1,x2,x3, y):
            return y-fcn(p, x1, x2, x3)

        def res4(p, x1,x2,x3,x4,  y):
            return y-fcn(p, x1, x2, x3, x4)

        resdic={1:res1,  2:res2,  3:res3,  4:res4}
        self.resfcn=resdic[len(datatuple)-1]
        i=0
        for arr in datatuple:  #if the numerical data is given as a list or tuple then convert to arrays. regardless convert to float64 because leastsq REQUIRES THIS
            datatuple=datatuple[0:i]+tuple([numpy.float64(arr)])+datatuple[i+1:]
            i=i+1
        while self.performfit:
            self.sigmas=scipy.zeros(len(self.finalparams))
            if not optimizerfcn is None:
                try:
                    self.finalparams=optimizerfcn(self.resfcn,self.initparams, args=datatuple, maxfun=self.maxfev, xtol=1.e-10, ftol=1.e-10)
                    self.error=0
                except:
                    self.error=1
            else:
                fitout = scipy.optimize.leastsq(self.resfcn,self.initparams, args=datatuple, maxfev=self.maxfev, full_output=1)#, warning=False)
                self.performfit=False
                self.finalparams=fitout[0]
                if not fitout[4] in [1, 2]:
                    print 'Fitting Error ', fitout[4], ' at ', markstr,': ', fitout[3]
                    self.error=True
                else:
                    #self.finalparams=fitout[0]
                    self.covmat=fitout[1]
                    try:
                        self.sigmas=scipy.array([self.covmat[i, i] for i in range(len(self.sigmas))])
                    except:
                        pass

        def fitfcn(x):
            return fcn(self.finalparams, x)
        return fitfcn

    def poly(self, p, x):#both must be numpy arrays
        return numpy.array([p[i]*(x**i) for i in range(p.size)]).sum(0)

    def polyfit(self, datatuple, initparams, markstr='unspecified', interaction=0,  maxfev=2000, weights=None):
        #initparams can be an array of coefficients [constant,lin term, quad term,...] or an integer indicating the order of the polynomial
        if isinstance(initparams, int):
            initparams=numpy.ones(initparams+1)
        else:
            initparams=numpy.float64(initparams)
        parnames=[]
        i=0
        for par in initparams:
            parnames+=[''.join(('coef', `i`))]
            i+=1

        return self.genfit(self.poly, initparams, datatuple, markstr, parnames, interaction, maxfev, weights=weights)

    def gaussianfit(self, datatuple, initparams=scipy.array([1, 0, 1]), markstr='unspecified', interaction=0, showplot=True, maxfev=2000, weights=None):
        return self.genfit(self.gaussian, initparams, datatuple, markstr, parnames=['coef', 'center', 'sigma'], interaction=interaction, maxfev=maxfev, weights=weights)

    def gaussian(self, p, x):
        return p[0]*scipy.exp(-0.5*((x-p[1])/p[2])**2)

    def lorentzianfit(self, datatuple, initparams=scipy.array([1, 0, 1]), markstr='unspecified', interaction=0, showplot=True, maxfev=2000, weights=None):
        return self.genfit(self, self.lorentzian, initparams, datatuple, markstr, parnames=['coef', 'center', 'gamma'], interaction=interaction, maxfev=maxfev, weights=weights)

    def lorentzian(self, p, x):
        return (p[0]/scipy.pi)*p[2]/((x-p[1])**2+p[2]**2)

def Gaussian(pars, x):
    return pars[2]*numpy.exp(-0.5*((x-pars[0])/pars[1])**2)

def Lorentzian(pars, x):#defined in nontraditional way so that pars[2] is the peak height
    return pars[2]/(1+((x-pars[0])/pars[1])**2)

def GaussLorentz(pars, x):
    gw=min(max(pars[3], 0.), 1.)
    return gw*Gaussian(pars, x)+(1.-gw)*Lorentzian(pars, x)

def GaussHalfLorentz(pars, x):
    return .5*Gaussian(pars, x)+.5*Lorentzian(pars, x)

PeakFcnLibrary={'Gaussian':Gaussian, 'Lorentzian':Lorentzian, 'GaussHalfLorentz':GaussHalfLorentz}

def fitpeakset(X, Y, initpars, peakfcn, negpeaks=True, optimizerfcn=None, nsigcut=3.):#peak function must be a function that accepts a list of 3 parameters (the reshape 3 needs to be changed if num params differs)
    numgauss=len(initpars)
    if numgauss==0:
        return (numpy.float32([]), numpy.float32([]), 0.)
    if nsigcut is None:
        imin=0
        imax=len(X)
    else:
        xmin=initpars[0][0]
        xmax=initpars[0][0]
        for p, w, h in initpars:
            xmin=min(xmin, p-w*nsigcut)
            xmax=max(xmax, p+w*nsigcut)
        imin=numpy.argmin((X-xmin)**2)
        imax=numpy.argmin((X-xmax)**2)

    zeroedpeakinds=[]
    repeatwithpkremoved=True #peaks are removed if their fitted height is <0. At the end, these peaks are added to the fit parameter list with 0 height and 0 error
    while repeatwithpkremoved:
        initparscpy=copy.copy(list(initpars))
        for pkind in reversed(zeroedpeakinds):#reverse so opo gets the right index
            initparscpy.pop(pkind)
        if len(initparscpy)==0:
            break
        initparsflat=numpy.float64(initparscpy).flatten()
        def fitfcn(p, x):
            allpars=numpy.reshape(p, (p.size//initpars.shape[1], initpars.shape[1]))
            if isinstance(x, numpy.ndarray):
                val=numpy.zeros(x.size, dtype='float32')
            else:
                val=0.0
            for pars in allpars:
                val+=peakfcn(pars, x)
            return val
#        def residfcn(p, x, y):
#            err=y-fitfcn(p, x)
#            return err
        Ya=numpy.float64(Y[imin:imax])
        Xa=numpy.float64(X[imin:imax])

        #if not optimizerfcn is None:
        ff=fitfcns()
        ff.genfit(fitfcn, initparsflat, (Xa, Ya), optimizerfcn=optimizerfcn)
        finalparams=ff.finalparams
#        else:
#            fitout=scipy.optimize.leastsq(residfcn, initparsflat, args=(X, Y), full_output=1)
#            if not (fitout[4] in [1, 2]):
#                print 'Fitting Error', fitout[4],': ', fitout[3]
#            finalparams=numpy.float32(fitout[0])
        finalparamsshaped=numpy.reshape(finalparams, (len(finalparams)//initpars.shape[1], initpars.shape[1]))
        if negpeaks:
            repeatwithpkremoved=False
        else:
            negpeakinds=numpy.where(finalparamsshaped[:, 2]<0)[0]
            zeroedpeakinds+=list(negpeakinds)
            zeroedpeakinds.sort()
            repeatwithpkremoved=len(negpeakinds)>0
#        print '^^^^^^^^^^^^^^^'
#        print initparsflat
#        print finalparamsshaped
#        pylab.plot(X, Y, 'b.')
#        pylab.show()
#    if not (fitout[1] is None):
#        covmat=fitout[1]
#        sigmas=numpy.float32([covmat[i, i] for i in range(len(finalparams))])
#    else:
#        print 'COVARIANCE NOT CALCULATED:', fitout[4],': ', fitout[3]
#        sigmas=numpy.zeros(len(finalparams), dtype='float32')
        sigmas=ff.sigmas
    finalresid=numpy.sqrt((ff.resfcn(finalparams, X, Y)**2).sum())
    #pylab.plot(X, Y, 'k.', X, fitfcn(finalparams, X), 'r-')

    sigmashaped=numpy.reshape(sigmas, (len(finalparams)//initpars.shape[1], initpars.shape[1]))
    for pkind in zeroedpeakinds:
        finalparamsshaped=list(finalparamsshaped)
        sigmashaped=list(sigmashaped)
        temp=copy.copy(initpars[pkind][:])
        temp[2]=0.#zero the height
        finalparamsshaped.insert(pkind, temp)
        sigmashaped.insert(pkind, numpy.zeros(initpars.shape[1], dtype='float64'))
        finalparamsshaped=numpy.float64(finalparamsshaped)
        sigmashaped=numpy.float64(sigmashaped)
    return (finalparamsshaped, sigmashaped, finalresid)

def arrayzeroind1d(arr, postoneg=False, negtopos=True):
    sarr=numpy.sign(arr)
    if postoneg:
        zeroind=numpy.where(sarr[:-1]>sarr[1:])[0]
        if negtopos:
            zeroind=numpy.append(zeroind, numpy.where(sarr[:-1]*sarr[1:]<=0)[0])
    else:#assume that if not postoneg then negtopos
        zeroind=numpy.where(sarr[:-1]*sarr[1:]<=0)[0]
    return (1.0*zeroind*arr[(zeroind+1,)]-(zeroind+1)*arr[(zeroind,)])/(arr[(zeroind+1,)]-arr[(zeroind,)]) #returns array of the floating point "index" linear interpolation between 2 indeces

def clustercoordsbymax1d(arr, pkind, critsepind):#results will be sorted. wherever there are peak indeces too close together. the peak index next to the peak index with highest arr value gets removed
    pkind.sort()
    indindslow=numpy.where((pkind[1:]-pkind[:-1])<critsepind)[0]
    indindshigh=indindslow+1
    while indindslow.size>0:
        maxindindindlow=numpy.nanargmax(arr[pkind[(indindslow,)]])
        maxindindindhigh=numpy.nanargmax(arr[pkind[(indindshigh,)]])
        if arr[pkind[indindslow[maxindindindlow]]]>arr[pkind[indindshigh[maxindindindhigh]]]:
            pkind=numpy.delete(pkind, indindshigh[maxindindindlow])
        else:
            pkind=numpy.delete(pkind, indindslow[maxindindindhigh])

        indindslow=numpy.where((pkind[1:]-pkind[:-1])<critsepind)[0]
        indindshigh=indindslow+1
    return pkind

def peaksearch1dSG(x, dx=1., critpeakheight=10, critsepind=5, critcurve=None, firstdernpts=7, firstderorder=1, secdernpts=14, secderorder=1, pospeaks=True, negpeaks=True):
    #dx is delta q for one index. zeros of the first derivative of inn are grouped together if within critsepind. only negative slope in the firstder is used so no secder is necessary unless specify a critical curvature in count nm^2
    if not (pospeaks or negpeaks):
        return numpy.float32([])
    ifirstder=savgolsmooth(x, nptsoneside=firstdernpts, order=firstderorder, dx=dx, deriv=1)
    fullpkind=numpy.float32([])
    if pospeaks:
        zeroind=arrayzeroind1d(ifirstder, postoneg=True, negtopos=False)
        temp=numpy.where(x[(numpy.uint32(numpy.round(zeroind)),)]>critpeakheight)
        fullpkind=numpy.append(fullpkind, zeroind[temp])
    if negpeaks:
        zeroind=arrayzeroind1d(ifirstder, postoneg=False, negtopos=True)
        temp=numpy.where(x[(numpy.uint32(numpy.round(zeroind)),)]<(-1*critpeakheight))
        fullpkind=numpy.append(fullpkind, zeroind[temp])

    if fullpkind.size==0:
        return fullpkind
    pkind=clustercoordsbymax1d(x, numpy.uint32(numpy.round(fullpkind)), critsepind)
    if critcurve is not None:
        isecder=savgolsmooth(x, nptsoneside=secdernpts, order=secderorder, dx=dx, deriv=2)
        temp=numpy.where(numpy.abs(isecder[(numpy.uint32(numpy.round(pkind)),)])>(critcurve))
        pkind=numpy.array(pkind)[temp]
#    pkind=list(pkind)
#    pkind.reverse()#highest to smallest for pairing below
    return numpy.array(pkind, dtype=numpy.float32)

def removeoutliers_meanstd(arr, nptsoneside, nsig, gapptsoneside=0): #avrages maximum of 2*nptoneside points and usees distance from mean scaled by std compared to nsig to determine if the value should be replaced by the mean. if gapptsoneside>0, will do this leaving a gap around the point in question and using nptsoneside-gaps points for the mean and std
    if nptsoneside==1 and gapptsoneside==0:
        return removesinglepixoutliers(arr, critratiotoneighbors=nsig)
    nsig=max(nsig, 1.)
    nptsoneside=max(nptsoneside, 2.)
    gapptsoneside=min(gapptsoneside, nptsoneside-2.)
    for gap in range(gapptsoneside+1):
        starti=numpy.uint32([max(i-(nptsoneside-gap), 0) for i in range(len(arr))])
        stopi=numpy.uint32([min(i+(nptsoneside-gap)+1, len(arr)) for i in range(len(arr))])
        #print [numpy.append(arr[i0:i], arr[i+1:i1]) for i, i0, i1 in zip(range(len(arr)), starti, stopi)][8]
        #print [(((numpy.append(arr[i0:i], arr[i+1:i1]).mean()-arr[i]))**2, (numpy.append(arr[i0:i], arr[i+1:i1]).std()*nsig)**2) for i, i0, i1 in zip(range(len(arr)), starti, stopi)][8]
        arr=numpy.array([(((numpy.append(arr[i0:i], arr[i+1:i1]).mean()-arr[i]))**2<(numpy.append(arr[i0:i], arr[i+1:i1]).std()*nsig)**2 and (arr[i],) or (numpy.append(arr[i0:i], arr[i+1:i1]).mean(),))[0] for i, i0, i1 in zip(range(len(arr)), starti, stopi)], dtype=arr.dtype)
    return arr

def removesinglepixoutliers(arr,critratiotoneighbors=1.5):
    c=numpy.where((arr[1:-1]>(critratiotoneighbors*arr[:-2]))*(arr[1:-1]>(critratiotoneighbors*arr[2:])))
    c0=c[0]+1
    #print len(c0), ' pixels being replaced'
    arr[c0]=(arr[c0-1]+arr[c0+1])/2
    return arr

def clustercoordsbymax1d(arr, pkind, critqsepind):#results will be sorted. wherever there are peak indeces too close together. the peak index next to the peak index with highest arr value gets removed
    pkind.sort()
    indindslow=numpy.where((pkind[1:]-pkind[:-1])<critqsepind)[0]
    indindshigh=indindslow+1
    while indindslow.size>0:
        maxindindindlow=numpy.argmax(arr[pkind[(indindslow,)]])
        maxindindindhigh=numpy.argmax(arr[pkind[(indindshigh,)]])
        if arr[pkind[indindslow[maxindindindlow]]]>arr[pkind[indindshigh[maxindindindhigh]]]:
            pkind=numpy.delete(pkind, indindshigh[maxindindindlow])
        else:
            pkind=numpy.delete(pkind, indindslow[maxindindindhigh])

        indindslow=numpy.where((pkind[1:]-pkind[:-1])<critqsepind)[0]
        indindshigh=indindslow+1
    return pkind

def findlinearsegs(y, dydev_frac,  dydev_nout, dn_segstart,  SGnpts=10, plotbool=False, dx=1., dydev_abs=0., maxfracoutliers=.5, critdy_fracmaxdy=None, critdy_abs=None, npts_SGconstdy=2):
    if 2*npts_SGconstdy+dydev_nout>=len(y):
        print 'array not long enough to find linear segments'
        return [], [], [], [], []
    dy=savgolsmooth(y, nptsoneside=SGnpts, order = 2, dx=dx, deriv=1)
#    lenconstdy=numpy.array([(dy[i]==0. and (0, ) or \
#    (numpy.all((numpy.abs((dy[i:]-dy[i])/dy[i])<=dydev_frac)|(numpy.abs(dy[i:]-dy[i])<=dydev_abs)) and (len(dy)-i, ) or \
#    (numpy.where(numpy.logical_not((numpy.abs((dy[i:]-dy[i])/dy[i])<=dydev_frac)|(numpy.abs(dy[i:]-dy[i])<=dydev_abs)))[0][:dydev_nout][-1],)))[0] for i in range(len(dy)-dydev_nout)])
    lenconstdy=numpy.array([(dy[i]==0. and (0, ) or \
    (numpy.all(numpy.abs(dy[i:]-dy[i])<max(numpy.abs(dy[i]*dydev_frac), dydev_abs)) and (len(dy)-i, ) or \
    (numpy.where(numpy.logical_not(numpy.abs(dy[i:]-dy[i])<max(numpy.abs(dy[i]*dydev_frac), dydev_abs)))[0][:dydev_nout][-1],)))[0] for i in range(len(dy)-dydev_nout)])
    #print len(y), len(lenconstdy), lenconstdy.max()

    if len(lenconstdy)==0:
        len_segs=[]
        istart_segs=[]
    else:
        lendn=savgolsmooth(numpy.float32(lenconstdy), nptsoneside=npts_SGconstdy, order = 1, dx=1.0, deriv=1, binprior=0)
        if plotbool:
            pylab.figure()
            pylab.plot(lenconstdy)
            pylab.ylabel('initial len of consecutive "constant dy" points')
            pylab.figure()
            pylab.plot(lendn)
            pylab.ylabel('lendn = deriv of len of const dy pts')
        istart_segs=numpy.where((lendn[:-1]>0)&(lendn[1:]<0))[0]
        if numpy.any(lenconstdy[:npts_SGconstdy+1]>=lenconstdy[npts_SGconstdy+1]):
            itemp=numpy.argmax(lenconstdy[:npts_SGconstdy])
            if not itemp in istart_segs:
                istart_segs=numpy.append(itemp, istart_segs)
        istart_segs[(istart_segs<npts_SGconstdy*2)]=npts_SGconstdy*2
        istart_segs[(istart_segs>len(y)-1-npts_SGconstdy*2)]=len(y)-1-npts_SGconstdy*2
        istart_segs+=numpy.array([numpy.argmax(lenconstdy[i-npts_SGconstdy*2:i+npts_SGconstdy*2]) for i in istart_segs])-npts_SGconstdy*2
        istart_segs=numpy.array(clustercoordsbymax1d(lenconstdy, istart_segs, dn_segstart))
        istart_segs=istart_segs[lenconstdy[istart_segs]>=dydev_nout/maxfracoutliers]

        #
        #istart=numpy.array([i for i in istart if (numpy.abs(i-istart)<dn_segstart).sum()==1 or numpy.all(lenconstdy[numpy.abs(i-istart)<dn_segstart]<=lenconstdy[i])])
        #istart=numpy.array([((numpy.abs(i-istart)<dn_segstart).sum()==1 and (i,) or (numpy.median(, ))[0] for i in istart)
        len_segs=lenconstdy[istart_segs]

    if not critdy_abs is None:
        critdy_fracmaxdy=critdy_abs/numpy.abs(dy).max()
    if not critdy_fracmaxdy is None:
        istart_constsegs, len_constsegs, garb, garb=findzerosegs(dy, critdy_fracmaxdy,  dydev_nout, dn_segstart,  SGnpts=SGnpts, plotbool=plotbool, dx=1., maxfracoutliers=maxfracoutliers)
#        if len(istart_constsegs)>0:
#            print istart_constsegs, len_constsegs
        temp=[[i, l] for i, l in zip(istart_constsegs, len_constsegs) if numpy.min((istart_segs-i)**2)>dn_segstart**2]
        if len(temp)>0:
            istart_constsegs, len_constsegs=numpy.array(temp).T
            istart_segs=numpy.append(istart_segs, istart_constsegs)
            len_segs=numpy.append(len_segs, len_constsegs)
    if len(istart_segs)==0:
        return numpy.array([]), numpy.array([]), numpy.array([]), numpy.array([]), dy

    #dy_segs=numpy.array([dy[i:i+l].mean() for i, l in zip(istart_segs, len_segs)])
    #interc_segs=numpy.array([(y[i:i+l]-(i+numpy.arange(l))*d).mean() for i, l, d in zip(istart_segs, len_segs, dy_segs)])

    fitdy_segs, fitinterc_segs=numpy.array([numpy.polyfit(dx*(i+numpy.arange(l)), y[i:i+l], 1) for i, l in zip(istart_segs, len_segs)]).T

    if plotbool:
        pylab.figure()
        pylab.plot(y)
        cols=['g', 'c', 'y', 'm']
        for count, (i, l, fd, fy0) in enumerate(zip(istart_segs, len_segs, fitdy_segs, fitinterc_segs)):
            if l==max(len_segs):
                c='r'
            else:
                c=cols[count%len(cols)]
            pylab.plot(i, fy0+i*fd*dx, 'g*')
            #pylab.plot(i+numpy.arange(l),numpy.arange(l)*d+y0, 'g-')
            pylab.plot(i+numpy.arange(l),(i+numpy.arange(l))*fd*dx+fy0, '-', color=c, lw=2)

        pylab.figure()
        pylab.plot(dy)
        pylab.ylabel('dy')
        for count, (i, l, fd) in enumerate(zip(istart_segs, len_segs, fitdy_segs)):
            if l==max(len_segs):
                c='r'
            else:
                c=cols[count%len(cols)]
            pylab.plot(i, fd, 'g*')
            pylab.plot(i+numpy.arange(l),numpy.ones(l)*fd, '-', color=c, lw=2)
        pylab.figure()
        pylab.plot(lenconstdy, 'ko')
        pylab.plot(istart_segs, len_segs, 'r.')
        pylab.ylabel('len of consecutive "constant dy" points ')
    return istart_segs, len_segs, fitdy_segs, fitinterc_segs, dy # fit intercept is wrt the beginning of the array, index=0 not x=0

def findmatchinglinearsegs(y, dydev_frac,  dydev_nout, dn_segstart,  SGnpts=10, plotbool=False, dx=1., dydev_abs=0., maxfracoutliers=.5, critdy_fracmaxdy=None, critdy_abs=None):
    dy=savgolsmooth(y, nptsoneside=SGnpts, order = 2, dx=dx, deriv=1)
#    lenconstdy=numpy.array([(numpy.any(dyall[:,i]==0.) and (0,) or \
#    (numpy.all([(numpy.abs((dy[i:]-dyall[:,i].mean())/dyall[:,i].mean())<=dydev_frac)|(numpy.abs(dy[i:]-dy[i])<=dydev_abs)) and (len(dy)-i, ) or \
#    (numpy.where(numpy.logical_not((numpy.abs((dy[i:]-dy[i])/dy[i])<=dydev_frac)|(numpy.abs(dy[i:]-dy[i])<=dydev_abs)))[0][:dydev_nout][-1],)))[0] for i in range(len(dy)-dydev_nout)])

    nptstemp=2
    lendn=savgolsmooth(numpy.float32(lenconstdy), nptsoneside=nptstemp, order = 1, dx=1.0, deriv=1, binprior=0)
    istart_segs=numpy.where((lendn[:-1]>0)&(lendn[1:]<0))[0]
    istart_segs[(istart_segs<nptstemp*2)]=nptstemp*2
    istart_segs[(istart_segs>len(y)-1-nptstemp*2)]=len(y)-1-nptstemp*2
    istart_segs+=numpy.array([numpy.argmax(lenconstdy[i-nptstemp*2:i+nptstemp*2]) for i in istart_segs])-nptstemp*2
    istart_segs=numpy.array(clustercoordsbymax1d(lenconstdy, istart_segs, dn_segstart))
    istart_segs=istart_segs[lenconstdy[istart_segs]>=dydev_nout/maxfracoutliers]

    #
    #istart=numpy.array([i for i in istart if (numpy.abs(i-istart)<dn_segstart).sum()==1 or numpy.all(lenconstdy[numpy.abs(i-istart)<dn_segstart]<=lenconstdy[i])])
    #istart=numpy.array([((numpy.abs(i-istart)<dn_segstart).sum()==1 and (i,) or (numpy.median(, ))[0] for i in istart)
    len_segs=lenconstdy[istart_segs]

    if not critdy_abs is None:
        critdy_fracmaxdy=critdy_abs/numpy.abs(dy).max()
    if not critdy_fracmaxdy is None:
        istart_constsegs, len_constsegs, garb, garb=findzerosegs(dy, critdy_fracmaxdy,  dydev_nout, dn_segstart,  SGnpts=SGnpts, plotbool=plotbool, dx=1., maxfracoutliers=maxfracoutliers)
#        if len(istart_constsegs)>0:
#            print istart_constsegs, len_constsegs
        temp=[[i, l] for i, l in zip(istart_constsegs, len_constsegs) if numpy.min((istart_segs-i)**2)>dn_segstart**2]
        if len(temp)>0:
            istart_constsegs, len_constsegs=numpy.array(temp).T
            istart_segs=numpy.append(istart_segs, istart_constsegs)
            len_segs=numpy.append(len_segs, len_constsegs)
    if len(istart_segs)==0:
        return numpy.array([]), numpy.array([]), numpy.array([]), numpy.array([]), dy

    #dy_segs=numpy.array([dy[i:i+l].mean() for i, l in zip(istart_segs, len_segs)])
    #interc_segs=numpy.array([(y[i:i+l]-(i+numpy.arange(l))*d).mean() for i, l, d in zip(istart_segs, len_segs, dy_segs)])

    fitdy_segs, fitinterc_segs=numpy.array([numpy.polyfit(dx*(i+numpy.arange(l)), y[i:i+l], 1) for i, l in zip(istart_segs, len_segs)]).T

    if plotbool:
        pylab.figure()
        pylab.plot(y)
        cols=['g', 'c', 'y', 'm']
        for count, (i, l, fd, fy0) in enumerate(zip(istart_segs, len_segs, fitdy_segs, fitinterc_segs)):
            if l==max(len_segs):
                c='r'
            else:
                c=cols[count%len(cols)]
            pylab.plot(i, fy0+i*fd*dx, 'g*')
            #pylab.plot(i+numpy.arange(l),numpy.arange(l)*d+y0, 'g-')
            pylab.plot(i+numpy.arange(l),(i+numpy.arange(l))*fd*dx+fy0, '-', color=c, lw=2)

        pylab.figure()
        pylab.plot(dy)
        pylab.ylabel('dy')
        for count, (i, l, fd) in enumerate(zip(istart_segs, len_segs, fitdy_segs)):
            if l==max(len_segs):
                c='r'
            else:
                c=cols[count%len(cols)]
            pylab.plot(i, fd, 'g*')
            pylab.plot(i+numpy.arange(l),numpy.ones(l)*fd, '-', color=c, lw=2)
        pylab.figure()
        pylab.plot(lenconstdy, 'ko')
        pylab.plot(istart_segs, len_segs, 'r.')
        pylab.ylabel('len of consecutive "constant dy" points ')
    return istart_segs, len_segs, fitdy_segs, fitinterc_segs, dy # fit intercept is wrt the beginning of the array, index=0 not x=0

def findzerosegs(y, yzero_maxfrac,  ydev_nout, dn_segstart,  SGnpts=10, plotbool=False, dx=1., maxfracoutliers=.5):#ydev_nout is number of outliers allowed in segment, dn_segstart is how close to each other the segments are allowed to start
    if ydev_nout>=len(y):
        print 'array not long enough to find zero segments'
        return [], [], [], []
    y=savgolsmooth(y, nptsoneside=SGnpts, order = 2)
    yzero_maxfrac=numpy.abs(y).max()*yzero_maxfrac
    lenzeroy=numpy.array([\
    (numpy.all(numpy.abs(y[i:])<=yzero_maxfrac) and (len(y)-i, ) or \
    (numpy.where(numpy.abs(y[i:])>yzero_maxfrac)[0][:ydev_nout][-1],))[0] for i in range(len(y)-ydev_nout)])

    nptstemp=2
    lendn=savgolsmooth(numpy.float32(lenzeroy), nptsoneside=nptstemp, order = 1, dx=1.0, deriv=1, binprior=0)

    istart_segs=numpy.where((lendn[:-1]>0)&(lendn[1:]<0))[0]
    if numpy.any(lenzeroy[:nptstemp+1]>=lenzeroy[nptstemp+1]):
        itemp=numpy.argmax(lenzeroy[:nptstemp])
        if not itemp in istart_segs:
            istart_segs=numpy.append(itemp, istart_segs)
    istart_segs[(istart_segs<nptstemp*2)]=nptstemp*2
    istart_segs[(istart_segs>len(y)-1-nptstemp*2)]=len(y)-1-nptstemp*2
    istart_segs+=numpy.array([numpy.argmax(lenzeroy[i-nptstemp*2:i+nptstemp*2]) for i in istart_segs])-nptstemp*2
    istart_segs=numpy.array(clustercoordsbymax1d(lenzeroy, istart_segs, dn_segstart))
    istart_segs=istart_segs[lenzeroy[istart_segs]>ydev_nout/maxfracoutliers]

    if len(istart_segs)==0:
        return numpy.array([]), numpy.array([]), numpy.array([]), numpy.array([])

    #istart=numpy.array([i for i in istart if (numpy.abs(i-istart)<dn_segstart).sum()==1 or numpy.all(lenconstdy[numpy.abs(i-istart)<dn_segstart]<=lenconstdy[i])])
    #istart=numpy.array([((numpy.abs(i-istart)<dn_segstart).sum()==1 and (i,) or (numpy.median(, ))[0] for i in istart)
    len_segs=lenzeroy[istart_segs]

    #dy_segs=numpy.array([dy[i:i+l].mean() for i, l in zip(istart_segs, len_segs)])
    #interc_segs=numpy.array([(y[i:i+l]-(i+numpy.arange(l))*d).mean() for i, l, d in zip(istart_segs, len_segs, dy_segs)])

    fitdy_segs, fitinterc_segs=numpy.array([numpy.polyfit(dx*(i+numpy.arange(l)), y[i:i+l], 1) for i, l in zip(istart_segs, len_segs)]).T

    if plotbool:
        pylab.figure()
        pylab.plot(y)
        cols=['g', 'c', 'y', 'm']
        for count, (i, l, fd, fy0) in enumerate(zip(istart_segs, len_segs, fitdy_segs, fitinterc_segs)):
            if l==max(len_segs):
                c='r'
            else:
                c=cols[count%len(cols)]
            pylab.plot(i, fy0+i*fd*dx, 'g*')
            #pylab.plot(i+numpy.arange(l),numpy.arange(l)*d+y0, 'g-')
            pylab.plot(i+numpy.arange(l),(i+numpy.arange(l))*fd*dx+fy0, '-', color=c, lw=2)

        for count, (i, l, fd) in enumerate(zip(istart_segs, len_segs, fitdy_segs)):
            if l==max(len_segs):
                c='r'
            else:
                c=cols[count%len(cols)]
            pylab.plot(i, fd, 'g*')
            pylab.plot(i+numpy.arange(l),numpy.ones(l)*fd, '-', color=c, lw=2)
        pylab.figure()
        pylab.plot(lenzeroy, 'ko')
        pylab.plot(istart_segs, lenzeroy[istart_segs], 'r.')
        pylab.ylabel('len of consecutive "zero" points')
    return istart_segs, len_segs, fitdy_segs, fitinterc_segs # fit intercept is wrt the beginning of the array, index=0 not x=0

def reggrid(x, y):
    if not numpy.all((x[:-1]-x[1:])==(x[1]-x[0])):
        xg=numpy.linspace(x.min(), x.max(), len(x))
        yg=interp(xg, x, y)
        return xg, yg
    else:
        return x, y

def calcmeandt_dlist(dlist):
    for d in dlist:
        d['dt']=(d['t(s)'][1:]-d['t(s)'][:-1]).mean()

def calcmeandEdt_dlist(dlist):
    for d in dlist:
        d['dt']=(d['t(s)'][1:]-d['t(s)'][:-1]).mean()
        d['dE']=(numpy.abs(d['Ewe(V)'][1:]-d['Ewe(V)'][:-1])).mean()
        d['dEdt']=d['dE']/d['dt']

def calcsegind_dlist(dlist, SGpts=8):
    if not 'dt' in dlist[0].keys():
        calcmeandt_dlist(dlist)
    for d in dlist:
        d['Ewe(V)_dtSG']=savgolsmooth(d['Ewe(V)'], nptsoneside=SGpts, order = 1, dx=d['dt'], deriv=1, binprior=0)
        rising=d['Ewe(V)_dtSG']>=0
        d['segind']=numpy.empty(len(rising), dtype='uint16')
        endptcorr=int(1.5*SGpts)#remove possibiliies of segment breaks within 1.5SGpts of the edges
        rising[:endptcorr+1]=rising[endptcorr+2]
        rising[-endptcorr-1:]=rising[-endptcorr-2]
        inds=numpy.where(rising[:-1]!=rising[1:])[0]
        inds=numpy.concatenate([[0], inds, [len(rising)]])
        for count, (i0, i1) in enumerate(zip(inds[:-1], inds[1:])):
            d['segind'][i0:i1]=count
        d['segprops_dlist']=[]
        for si, i0 in zip(range(d['segind'].max()+1), inds[:-1]):
            d['segprops_dlist']+=[{}]
            d['segprops_dlist'][-1]['rising']=rising[i0+1]
            inds2=numpy.where(d['segind']==si)[0]
            d['segprops_dlist'][-1]['inds']=inds2
            d['segprops_dlist'][-1]['npts']=len(d['segprops_dlist'][-1]['inds'])
            d['segprops_dlist'][-1]['dEdt']=d['Ewe(V)_dtSG'][inds2][SGpts:-SGpts].mean()

def manualsegind_dlist(dlist, istart=[0], SGpts=8):
    if not 'dt' in dlist[0].keys():
        calcmeandt_dlist(dlist)

    for d in dlist:
        d['Ewe(V)_dtSG']=savgolsmooth(d['Ewe(V)'], nptsoneside=SGpts, order = 1, dx=d['dt'], deriv=1, binprior=0)
        d['segprops_dlist']=[]
        d['segind']=numpy.empty(len(d['Ewe(V)']), dtype='uint16')
        iend=istart[1:]+[len(d['Ewe(V)'])]
        for count, (i0, i1) in enumerate(zip(istart, iend)):
            d['segind'][i0:i1]=count
            d['segprops_dlist']+=[{}]
            inds2=numpy.arange(i0, i1)
            d['segprops_dlist'][-1]['inds']=inds2
            d['segprops_dlist'][-1]['npts']=i1-i0
            d['segprops_dlist'][-1]['dEdt']=d['Ewe(V)_dtSG'][inds2][SGpts:-SGpts].mean()
            d['segprops_dlist'][-1]['rising']=d['segprops_dlist'][-1]['dEdt']>0.

def SegSG_dlist(dlist, SGpts=10, order=1, k='I(A)'):
    kSG=k+'_SG'
    for d in dlist:
        d[kSG]=numpy.zeros(d[k].shape, dtype='float32')
        for segd in d['segprops_dlist']:
            inds=segd['inds']
            d[kSG][inds]=numpy.float32(savgolsmooth(d[k][inds], nptsoneside=SGpts, order = order, deriv=0, binprior=0))

def SegdtSG_dlist(dlist, SGpts=10, order=1, k='I(A)', dxk='dt'):
    kSG=k+'_dtSG'
    for d in dlist:
        if not kSG in d.keys() and dxk in d.keys():
            continue
        d[kSG]=numpy.zeros(d[k].shape, dtype='float32')
        for segd in d['segprops_dlist']:
            inds=segd['inds']
            d[kSG][inds]=numpy.float32(savgolsmooth(d[k][inds], nptsoneside=SGpts, order = order, deriv=1, binprior=0, dx=d[dxk]))

def calcIllum_LVsettings(d, lv4tuple, savekey='Illumcalc'):
    lv_dark, lv_ill, lv_duty, lv_period=lv4tuple
    illum=numpy.zeros(len(d['t(s)']), dtype='bool')
    indsill=numpy.where((d['t(s)']>lv_dark)&(d['t(s)']<=lv_ill))[0]
    till=d['t(s)'][indsill]
    till-=till[0]
    cycfrac=(till%lv_period)/lv_period
    illum[indsill[cycfrac<=lv_duty]]=1
    d[savekey]=illum

def calcdiff_stepill(d, interd, ikey='Illum', ykeys=['Ewe(V)'], xkeys=['t(s)', 'I(A)'], illfracrange=(.4, .95), darkfracrange=(.4, .95)):
    if isinstance(ikey, list) or isinstance(ikey, numpy.ndarray):
        calcIllum_LVsettings(d, ikey, savekey='Illumcalc')
        illum=d['Illumcalc']
    else:
        illum=d[ikey]!=0
    istart_len_calc=lambda startind, endind, fracrange: (startind+numpy.floor(fracrange[0]*(endind-startind)), numpy.ceil((fracrange[1]-fracrange[0])*(endind-startind)))
    riseinds=numpy.where(illum[1:]&numpy.logical_not(illum[:-1]))[0]+1
    fallinds=numpy.where(numpy.logical_not(illum[1:])&illum[:-1])[0]+1

    if len(fallinds)==0 and len(riseinds)==0:
        print 'insufficient light cycles'
        return 1
    if illum[0]:
        illstart=0
        illend=fallinds[0]
        darkstart=fallinds[0]
        if len(riseinds)==0:
            darkend=len(illum)
        else:
            darkend=riseinds[0]
    else:
        darkstart=0
        darkend=riseinds[0]
        illstart=riseinds[0]
        if len(fallinds)==0:
            illend=len(illum)
        else:
            illend=fallinds[0]

    ill_istart, ill_len=istart_len_calc(illstart, illend, illfracrange)
    dark_istart, dark_len=istart_len_calc(darkstart, darkend, darkfracrange)

    inds_ill=[range(int(ill_istart), int(ill_istart+ill_len))]
    inds_dark=[range(int(dark_istart), int(dark_istart+dark_len))]

    d['inds_ill']=inds_ill
    d['inds_dark']=inds_dark

    getillvals=lambda arr:numpy.array([arr[inds].mean() for inds in inds_ill])
    getdarkvals=lambda arr:numpy.array([arr[inds].mean() for inds in inds_dark])

    interd['rawselectinds']=numpy.int32(numpy.round(getdarkvals(numpy.arange(len(d[ikey])))))
    for k in xkeys+ykeys:#these are len 1 arrays that will get saved as itnermediate data
        interd[k+'_ill']=getillvals(d[k])
        interd[k+'_dark']=getdarkvals(d[k])
    for k in ykeys:
        #d[k+'_illdiffmean']=d[k+'_ill'][0]-d[k+'_dark'][0]
        interd[k+'_illdiff']=numpy.array([d[k+'_ill'][0]-d[k+'_dark'][0]])
    return 0

def calcdiff_choppedill(d, interd, ikey='Illum', ykeys=['I(A)'], xkeys=['t(s)', 'Ewe(V)'], illfracrange=(.4, .95), darkfracrange=(.4, .95)):
    if isinstance(ikey, list) or isinstance(ikey, numpy.ndarray):
        calcIllum_LVsettings(d, ikey, savekey='Illumcalc')
        illum=d['Illumcalc']
    else:
        illum=d[ikey]!=0
    istart_len_calc=lambda startind, endind, fracrange: (startind+numpy.floor(fracrange[0]*(endind-startind)), numpy.ceil((fracrange[1]-fracrange[0])*(endind-startind)))
    riseinds=numpy.where(illum[1:]&numpy.logical_not(illum[:-1]))[0]+1
    fallinds=numpy.where(numpy.logical_not(illum[1:])&illum[:-1])[0]+1
    if len(fallinds)==0 or len(riseinds)==0:
        print 'insufficient light cycles'
        return 1
    riseinds=riseinds[riseinds<fallinds[-1]]#only consider illum if there is a dark before and after
    fallinds=fallinds[fallinds>riseinds[0]]
    if len(fallinds)==0 or len(riseinds)==0:
        print 'insufficient light cycles'
        return 1
    ill_istart, ill_len=istart_len_calc(riseinds, fallinds, illfracrange)
    darkstart, darkend=numpy.where(numpy.logical_not(illum))[0][[0, -1]]
    dark_istart, dark_len=istart_len_calc(numpy.concatenate([[darkstart], fallinds]), numpy.concatenate([riseinds, [darkend]]), darkfracrange)

    #inds_ill=[range(int(i0), int(i0+ilen)) for i0, ilen in zip(ill_istart, ill_len)]
    #inds_dark=[range(int(i0), int(i0+ilen)) for i0, ilen in zip(dark_istart, dark_len)]

    indstemp=[(range(int(i0ill), int(i0ill+ilenill)), range(int(i0dark), int(i0dark+ilendark))) for i0ill, ilenill, i0dark, ilendark in zip(ill_istart, ill_len, dark_istart, dark_len) if ilenill>0 and ilendark>0]
    inds_ill=map(operator.itemgetter(0), indstemp)
    inds_dark=map(operator.itemgetter(1), indstemp)
    if dark_len[-1]>0:
        inds_dark+=[range(int(dark_istart[-1]), int(dark_istart[-1]+dark_len[-1]))]
    else:
        inds_ill=inds_ill[:-1]

    d['inds_ill']=inds_ill
    d['inds_dark']=inds_dark

    getillvals=lambda arr:numpy.array([arr[inds].mean() for inds in inds_ill])
    getdarkvals=lambda arr:numpy.array([arr[inds].mean() for inds in inds_dark])

    #use the average dark inds as rawselectinds - ideally this would be the il inds but there are one less ill and don't want to endfill the selectinds with nan
    interd['rawselectinds']=numpy.int32(numpy.round(getdarkvals(numpy.arange(len(d[ikey])))))

    for k in xkeys+ykeys:
        interd[k+'_ill']=numpy.append(getillvals(d[k]), numpy.nan)
        interd[k+'_dark']=getdarkvals(d[k])
    for k in ykeys:
        interd[k+'_illdiff']=interd[k+'_ill']-0.5*(interd[k+'_dark']+numpy.append(interd[k+'_dark'], numpy.nan)[1:])
#        d[k+'_illdiffmean']=numpy.mean(d[k+'_illdiff'])
#        d[k+'_illdiffstd']=numpy.std(d[k+'_illdiff'])
    return 0
def calcdiff_ill_caller(d, interd, ikey='Illum', thresh=0, **kwargs):
    if isinstance(ikey, list) or isinstance(ikey, numpy.ndarray):
        calcIllum_LVsettings(d, ikey, savekey='Illumcalc')
        illum=d['Illumcalc']
        ikey='Illumcalc'
    else:
        illum=d[ikey]>thresh
    riseinds=numpy.where(illum[1:]&numpy.logical_not(illum[:-1]))[0]+1
    fallinds=numpy.where(numpy.logical_not(illum[1:])&illum[:-1])[0]+1
    d['IllumBool']=illum
    try:
        if len(riseinds)==0 or len(fallinds)==0 or (len(riseinds)==1 and len(fallinds)==1 and riseinds[0]>fallinds[0]):
            err=calcdiff_stepill(d, interd, ikey=ikey, **kwargs)
        else:
            err=calcdiff_choppedill(d, interd, ikey='IllumBool', **kwargs)
    except:
        return 1
    return err

#def illumtimeshift(illum, t, tshift):
#    tmod=t-tshift
#    inds=[numpy.argmin((tv-tmod)**2) for tv in t]
#    return illum[inds]

def illumtimeshift(d, ikey, tkey, tshift):
    tmod=d[tkey]-tshift
    inds=[numpy.argmin((t-tmod)**2) for t in d[tkey]]
    d[ikey]=d[ikey][inds]

def calc_comps_multi_element_inks(platemapdlist, cels_set_ordered, conc_el_chan, key_append_conc='.Loading', key_append_atfrac='.AtFrac', tot_conc_label=None):
    #platemap should have standard A-H channels, at least as many as conc_el_chan.shape[1], conc_el_chan is by "new" element and by platemap channel, make the key_append a string to update dlist or None to skip
    inkchannelsused=['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H'][:conc_el_chan.shape[1]]#grab the channels used
    #with this update the platemap is in units of concentration instead of fraction of printer channel loading
    if key_append_conc is None:
        key_append_conc='TEMP'
    cels_set_ordered_conc=[k+key_append_conc for k in cels_set_ordered]
    [d.update(\
                    zip(cels_set_ordered_conc, (numpy.float32([d[k] for k in inkchannelsused])[numpy.newaxis, :]*conc_el_chan).sum(axis=1)))\
                              for d in platemapdlist]
    if key_append_atfrac is None:
        cels_set_ordered_atfrac=[]
    else:
        cels_set_ordered_atfrac=[k+key_append_atfrac for k in cels_set_ordered]
        concarr=numpy.array([[d[k] for k in cels_set_ordered_conc] for d in platemapdlist])
        if not tot_conc_label is None:
            [d.update([(tot_conc_label, carr.sum())]) for d, carr in  zip(platemapdlist, concarr)]
        comparr2d=concarr/concarr.sum(axis=1)[:, numpy.newaxis]
        [d.update(zip(cels_set_ordered_atfrac, comparr)) for d, comparr in  zip(platemapdlist, comparr2d)]
    if key_append_conc=='TEMP':
        for d in platemapdlist:
            for k in cels_set_ordered_conc:
                del d[k]
        cels_set_ordered_conc=[]
    return cels_set_ordered_conc, cels_set_ordered_atfrac

def ternary_comp_to_cart(a, b, c):
    return 1.-a-b/2., numpy.sqrt(3)*b/2.0

def quaternary_comp_to_cart(a, b, c, d):
    return 1.-a-b/2.-d/2., b/2.*(3.**.5)+d/2./(3.**.5), d*(2.**.5)/(3.**.5)

def filterbydistancefromline(arr_of_xyz, xyz1, xyz2, critdist, betweenpoints=True, invlogic=False, returnonlyinds=False, is_composition=False): #if N points in D dimmensions arr_ is NxD, xyz1 and 2 are D, critdist is D-dimmension Euclidean, lineparameter is 0 to 1 from xyz 1 to xyz2 and can be out of that range
    #see http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
    if is_composition:#probably cheaper to not do coordinate transform and just multiply by the distance ratio for the ndim, but this allows quat. to be done with 3D math
        arr_of_xyz/=arr_of_xyz.sum(axis=1)[:,numpy.newaxis]
        xyz1/=xyz1.sum()
        xyz2/=xyz2.sum()

        if xyz1.size==1:
            print 'unary compositions do not make sense'
        elif xyz1.size==2:
            arr_of_xyz=arr_of_xyz[:, 0]
            xyz1=xyz1[:1]
            xyz2=xyz2[:1]
        elif xyz1.size==3:
            arr_of_xyz=numpy.array(ternary_comp_to_cart(*arr_of_xyz.T)).T
            xyz1=numpy.array(ternary_comp_to_cart(*xyz1))
            xyz2=numpy.array(ternary_comp_to_cart(*xyz2))
        elif xyz1.size==4:
            arr_of_xyz=numpy.array(quaternary_comp_to_cart(*arr_of_xyz.T)).T
            xyz1=numpy.array(quaternary_comp_to_cart(*xyz1))
            xyz2=numpy.array(quaternary_comp_to_cart(*xyz2))
        else:
            print 'quinary or higher compositions not supported'
            return None
    else:
        if xyz1.size>3:
            print 'Only 1D,2D,3D supported'
            return None
    notnoninds=numpy.where(numpy.logical_not(numpy.isnan(arr_of_xyz)).prod(axis=1))[0]
    numnansamples=len(arr_of_xyz)-len(notnoninds)
    if numnansamples>0:
        arr_of_xyz=arr_of_xyz[notnoninds]
    distfromlin=numpy.array([numpy.linalg.norm(numpy.cross(xyz2-xyz1, xyz1-xyz))/numpy.linalg.norm((xyz2-xyz1)) for xyz in arr_of_xyz])
    if betweenpoints:
        lineparameter=numpy.array([-numpy.inner(xyz2-xyz1, xyz1-xyz)/numpy.linalg.norm((xyz2-xyz1))**2 for xyz in arr_of_xyz])
        if invlogic:
            inds=numpy.where(numpy.logical_not((distfromlin<=critdist) & (lineparameter>=0) & (lineparameter<=1)))[0]
        else:
            inds=numpy.where((distfromlin<=critdist) & (lineparameter>=0) & (lineparameter<=1))[0]
    else:
        lineparameter=numpy.array([-numpy.inner(xyz2-xyz1, xyz1-xyz)/numpy.linalg.norm((xyz2-xyz1))**2 for xyz in arr_of_xyz])
        if invlogic:
            inds=numpy.where(numpy.logical_not(distfromlin<=critdist))[0]
        else:
            inds=numpy.where(distfromlin<=critdist)[0]
    if returnonlyinds:
        return notnoninds[inds]
    else:
        if numnansamples>0:
            temp=numpy.ones(len(notnoninds)+numnansamples, dtype='float64')*numpy.nan
            temp[notnoninds]=distfromlin
            distfromlin=temp

            temp=numpy.ones(len(notnoninds)+numnansamples, dtype='float64')*numpy.nan
            temp[notnoninds]=lineparameter
            lineparameter=temp
            inds=notnoninds[inds]
        return {'select_inds':inds, 'dist_from_line':distfromlin, 'norm_dist_along_line':lineparameter}

def twotheta_q(q, wl=0.15418, units='deg'): #q is scattering vector in 1/nm.  wl is wavelength in nm. return float32 2Theta value in degrees unless units='rad'
    if units=='rad':
        return numpy.float32(2.0*numpy.arcsin(wl*q/(4.0*numpy.pi)))
    else:
        return numpy.float32(2.0*numpy.arcsin(wl*q/(4.0*numpy.pi))*180.0/numpy.pi)

def d_q(q): #q is scattering vector in 1/nm.  wl is wavelength in nm. return float32 2Theta value in degrees unless units='rad'
    return numpy.float32(2.0*numpy.pi/q)

def q_twotheta(twotheta, wl=0.15418, units='deg'): #units are those of twotheta, wl in nm, q in 1/nm
    if units=='deg':
        twotheta*=(numpy.pi/180.0)
    return 4*numpy.pi*numpy.sin(twotheta/2.0)/wl

def integrate_spectrum_with_piecewise_weights(spectrumx, spectrumy, weightsx, weightsy, below_val=0, above_val=0, return_completed_weights=False):
    #assume spectrum x and weights x are in increasing order
    completed_weights=numpy.zeros(len(spectrumy), dtype='float64')
    completed_weights[spectrumx<min(weightsx)]=below_val
    completed_weights[spectrumx>max(weightsx)]=above_val

    for x0, x1, y0, y1 in zip(weightsx[:-1], weightsx[1:],weightsy[:-1], weightsy[1:]):
        inds=numpy.where((spectrumx>=x0)&(spectrumx<=x1))[0]
        xarr=spectrumx[inds]
        yarr=y0+(y1-y0)*(xarr-x0)/(x1-x0)
        completed_weights[inds]=yarr
    newy=spectrumy*completed_weights
    if 0:
        plt.figure()
        plt.plot(spectrumx, spectrumy, 'k-')
        plt.plot(spectrumx, newy, 'r-')
        twax=plt.twinx()
        twax.plot(spectrumx, completed_weights, 'b-')
        twax.plot(weightsx, weightsy, 'bo')
        plt.show()
    integrated_val=numpy.trapz(newy, x=spectrumx)
    if return_completed_weights:
        return integrated_val, completed_weights
    return integrated_val

def linear_2d_extrapolation(xy, z, xy_to_extrapolate, num_pts_in_lin_fit=6):#returns array of 2d interp values using linear interp of num_pts_in_lin_fit (up to num_pts_in_lin_fit x 2 if need be). Arguments mimic griddata
    x, y=xy.T
    new_z=[]
    for xi, yi in xy_to_extrapolate:
        dist_tups=sorted([((xv-xi)**2+(yv-yi)**2,count) for count,(xv,yv,zv) in enumerate(zip(x,y,z)) if not np.isnan(zv)])
        for extra_pts in range(num_pts_in_lin_fit+1):#try up to twice as many points if necessary to get unique x and y values
            fitinds=[tup[1] for tup in dist_tups[:num_pts_in_lin_fit+extra_pts]]
            xf=x[fitinds]
            yf=y[fitinds]
            zf=z[fitinds]
            if ((xf[0]==xf).prod()+(yf[0]==yf).prod())==0:#if more than 1 x value and also more than 1 y value
                break
        if extra_pts==num_pts_in_lin_fit:#cannot perform 2d extrap here
            print 'extrap failed because locally data is in 1d'
            new_z+=[np.nan]
            continue
        A = np.vstack([xf, yf, np.ones(len(xf))]).T
        fitcoefs=np.linalg.lstsq(A, zf)[0]
        new_z+=[(fitcoefs*np.array([xi,yi,1.])).sum()]
    return np.float64(new_z)

def linear_1d_extrapolation(x, z, x_to_extrapolate, num_pts_in_lin_fit=6):
    new_z=[]
    for xi in x_to_extrapolate:
        dist_tups=sorted([((xv-xi)**2,count) for count,(xv,zv) in enumerate(zip(x,z)) if not np.isnan(zv)])
        for extra_pts in range(num_pts_in_lin_fit+1):#try up to twice as many points if necessary to get unique x and y values
            fitinds=[tup[1] for tup in dist_tups[:num_pts_in_lin_fit+extra_pts]]
            xf=x[fitinds]
            zf=z[fitinds]
            if (xf[0]==xf).prod()==0:#if more than 1 x value
                break
        if extra_pts==num_pts_in_lin_fit:#cannot perform 2d extrap here
            print 'extrap failed because locally data is in 1d'
            new_z+=[np.nan]
            continue
        c=np.polyfit(xf, zf, 1)
        f = np.poly1d(c)
        new_z+=[f(xi)]
    return np.float64(new_z)

max_stars_repo_name: johnmgregoire/JCAPDataProcess

max_issues_repo_name: johnmgregoire/JCAPDataProcess

max_forks_repo_name: johnmgregoire/JCAPDataProcess


text:
import numpy as np
def computeStrainTensorGreen(ref, moving):
    strain = np.zeros((ref.n_cells, 3, 3))
    triangles = ref.faces.reshape((-1, 3))[:, 1:]
    for i, t in enumerate(triangles):
        E = np.pinv(ref.points[t])@moving.points[t]
        strain[i] = (E.T@E - np.eye(3))/2
    return strain

def computeStrainCoefficientAlongDirection(strain, v):
    """
    Compute the green strain along a direction, with the correction to make it comparable to engineering strain.
    """
    return = np.sqrt(1 + 2*np.einsum('ni,nij,nj->n', v, strain, v) -1

ext: py

max_stars_repo_name: gbernardino/rvmep

max_issues_repo_name: gbernardino/rvmep

max_forks_repo_name: gbernardino/rvmep

reason: import numpy


text:
import cv2
import base as bs
import numpy as np

def quant4(_img):
    ret = np.zeros_like(_img, dtype='uint8')

    print(_img.shape)
    for i in range(3):
        ind = _img[..., i] < 64
        print(ind.shape)
        ret[..., i][ind] = 32
        ind = _img[..., i] < 128
        i_inserted = ret[..., i] == 0
        ret[..., i][np.logical_and(ind, i_inserted)] = 96
        ind = _img[..., i] < 192
        i_inserted = ret[..., i] == 0
        ret[..., i][np.logical_and(ind, i_inserted)] = 160
        i_inserted = ret[..., i] == 0
        ret[..., i][i_inserted] = 224

    return ret

def pool2d(_img, w_size=8):
    out = img.copy()

    h, w, c = _img.shape
    nh = int(h / w_size)
    nw = int(w / w_size)
    for y in range(nh):
        for x in range(nw):
            for c_i in range(c):
                out[w_size * y:w_size * (y + 1), w_size * x: w_size * (x + 1), c_i] = np.max(_img[w_size * y:w_size * (y + 1), w_size * x: w_size * (x + 1), c_i])

    return out

img = cv2.imread("imori.jpg")
q_img = pool2d(img)
bs.show_img(q_img)


text:
#
# onnx2ks - convert ONNX to Knossos
#
# References:
# https://github.com/onnx/onnx/blob/master/docs/IR.md
# https://github.com/onnx/onnx/blob/72b701f7a55cafa4b8ab66a21dc22da0905b2f4c/onnx/onnx.in.proto

# TODO:
# - Emit "def" FunctionProtos in prelude
# - Add pre/post comments to expressions, and emit argument names as comments
# - roundtrip test
# - Fix pylint on TensorProto.FLOAT

#%%

from itertools import chain
from types import SimpleNamespace

import sys
import re
import warnings

import numpy as np

import onnx
from onnx import TensorProto, AttributeProto

import onnx.numpy_helper
import onnx.helper
import onnx.optimizer

import onnxruntime.capi.onnxruntime_pybind11_state as ort
from onnxruntime.capi.onnxruntime_pybind11_state.schemadef import OpSchema

from ksc.utils import paren
from ksc.type import Type
from ksc.expr import Expr, Def, EDef, Rule, Const, Var, Lam, Call, Let, If, Assert
from ksc.expr import pystr

# Pretty printing
# Importing prettyprint to get the decorated printers for Expression and Type
import ksc.prettyprint  # pylint: disable=unused-import

# Import the prettyprinter routines we use explicitly in this file
from prettyprinter import cpprint, pprint, pformat

# Needed this in order to see the error messages when pprint fails
import warnings

warnings.filterwarnings("always")

############################################################################
ATTR_TYPE_TO_KS_TYPE = {
    OpSchema.AttrType.FLOAT: Type.Float,
    OpSchema.AttrType.INT: Type.Integer,
    OpSchema.AttrType.STRING: Type.String,
    OpSchema.AttrType.TENSOR: Type.Tensor(-1, Type.Float),
    OpSchema.AttrType.FLOATS: Type.Tensor(1, Type.Float),
    OpSchema.AttrType.INTS: Type.Tensor(1, Type.Integer),
    OpSchema.AttrType.STRINGS: Type.Tensor(1, Type.String),
    OpSchema.AttrType.TENSORS: Type.Tensor(1, Type.Tensor(-1, Type.Float)),
    OpSchema.AttrType.GRAPH: Type.Lam(None, None),
}

def onnxAttrType_to_Type(ty):
    """
    Convert ONNX AttrType to KS Type.
    # TODOS:
    #     ty.SPARSE_TENSOR
    #     ty.SPARSE_TENSORS
    #     ty.GRAPHS
    """
    assert isinstance(ty, OpSchema.AttrType)
    return ATTR_TYPE_TO_KS_TYPE[ty]

# See https://github.com/onnx/onnx/blob/master/onnx/mapping.py
TENSOR_TYPE_TO_KS_TYPE = {
    int(TensorProto.FLOAT): Type.Float,
    int(TensorProto.UINT8): Type.Integer,
    int(TensorProto.UINT16): Type.Integer,
    int(TensorProto.INT8): Type.Integer,
    int(TensorProto.INT16): Type.Integer,
    int(TensorProto.INT32): Type.Integer,
    int(TensorProto.INT64): Type.Integer,
    int(TensorProto.BOOL): Type.Bool,
    int(TensorProto.FLOAT16): Type.Float,
    int(TensorProto.DOUBLE): Type.Float,
    int(TensorProto.COMPLEX64): Type.Tuple(Type.Float, Type.Float),
    int(TensorProto.COMPLEX128): Type.Tuple(Type.Float, Type.Float),
    int(TensorProto.UINT32): Type.Integer,
    int(TensorProto.UINT64): Type.Integer,
    int(TensorProto.STRING): Type.String,
}

def onnxTensorType_to_Type(ety):
    ty = TENSOR_TYPE_TO_KS_TYPE.get(ety)
    if ty is not None:
        return ty
    raise NotImplementedError(f"type {ety}")

def onnxType_to_Type(proto):
    ety = onnxTensorType_to_Type(proto.tensor_type.elem_type)
    rank = len(proto.tensor_type.shape.dim)
    return Type.Tensor(rank, ety)

def get_value(init):
    """
    Get scalar value from a TensorProto
    """
    # https://github.com/onnx/onnx/blob/72b701f7a55cafa4b8ab66a21dc22da0905b2f4c/onnx/onnx.in.proto#L448
    # Find https://github.com/onnx/onnx/blob/72b701f7a55cafa4b8ab66a21dc22da0905b2f4c/onnx/mapping.py#L49

    a = onnx.numpy_helper.to_array(init)

    ksty = onnxTensorType_to_Type(init.data_type)

    if ksty == Type.Bool:
        return Const(bool(a))

    if ksty == Type.Float:
        return Const(float(a))

    if ksty == Type.Integer:
        return Const(int(a))

    field = onnx.mapping.STORAGE_TENSOR_TYPE_TO_FIELD[init.data_type]
    if not field:
        raise NotImplementedError(f"field Type {field} {init.data_type}, {init}")

    value = getattr(init, field)
    if not value:
        raise NotImplementedError(
            f"value {init.int32_data} {field} {value} {init.data_type}\n{init}"
        )

    return Const(value)

def mkVec(val):
    if isinstance(val, (list, np.ndarray)):
        return Call("Vec_init", [Const(v) for v in val])
    else:
        return Call("Vec_init", [Const(val)])

def Expr_from_TensorProto(val, name):
    """
    Make KS values from a TensorProto
    """
    if len(val.dims) == 1 and val.dims[0] < 16:
        # vec constructor
        return mkVec(onnx.numpy_helper.to_array(val))

    if len(val.dims) == 0:
        return get_value(val)

    nptype = onnx.mapping.TENSOR_TYPE_TO_NP_TYPE[val.data_type]
    return Call(
        f"load-from-onnx-{nptype}", [*(Const(x) for x in val.dims), Const(name)]
    )

def Expr_from_AttrVal(val):
    if isinstance(val, onnx.TensorProto):
        return Expr_from_TensorProto(val, "?Expr_from_AttrVal?")

    if isinstance(val, list):
        return mkVec(val)

    if isinstance(val, bytes):
        val = val.decode("ascii")
    return Const(val)

def Expr_from_Attr(attr: AttributeProto):
    # ty = onnxAttrType_to_Type(attr.type)
    a = onnx.helper.get_attribute_value(attr)
    return Expr_from_AttrVal(a)

def emit_inits(inits, body):
    """
    Emit initializers
    """
    for init in reversed(inits):
        var = useVar(init.name)

        value = Expr_from_TensorProto(init, var.name)

        body = Let(var, value, body)
    return body

def get_attribute_default_value(attr):
    if not hasattr(attr, "_default_value") or not attr._default_value:
        return None

    val = onnx.AttributeProto()
    _nbytes = val.ParseFromString(attr._default_value)
    val = onnx.helper.get_attribute_value(val)

    ty = onnxAttrType_to_Type(attr.type)
    if ty.is_scalar:
        return Expr_from_AttrVal(val)

    if ty.is_tensor:
        assert isinstance(val, list)
        if isinstance(val[0], bytes):
            val = [v.decode("ascii") for v in val]
        return mkVec(val)

    return val

def get_default_value(schema, attr):
    val = get_attribute_default_value(attr)
    if val != None:
        return val

    print(
        f"** Attribute {attr.name} of op {schema.name} has no default value -- special casing"
    )

    # TODO: Formalize this, at least into prelude
    if schema.name == "MaxPool":
        if attr.name == "dilations":
            return Expr_from_AttrVal([1, 1])

    if schema.name == "Conv" or schema.name == "ConvTranspose":
        if attr.name == "pads":
            return Expr_from_AttrVal(
                [-1, -1]
            )  # TODO: this should either match dims, or maybe be empty
        if attr.name == "output_shape":
            return Expr_from_AttrVal([-1, -1])

    if schema.name == "SoftmaxCrossEntropyLossGrad" and attr.name == "ignore_index":
        return Const(-1)

    if schema.name == "SoftmaxCrossEntropyLoss" and attr.name == "ignore_index":
        return Const(-1)

    print(schema.doc)

    raise NotImplementedError(f"No default value for {schema.name} attr {attr.name}")

def get_all_schemas():
    # Form schema group, take latest version
    schemas = dict()
    for schema in ort.get_all_operator_schema():
        name = schema.name if schema.domain == "" else schema.domain + "." + schema.name
        if name in schemas and schema.since_version < schemas[name].since_version:
            pass
        schemas[name] = schema
    return schemas

def encode_name(name: str):
    if re.match("^[0-9]", name):
        return "%" + name
    else:
        return name

def useVar(name: str):
    return Var(encode_name(name), None, False)

def defVar(name: str, ty: Type):
    return Var(encode_name(name), ty, True)

def onnx2ks(g):
    """
    Convert ONNX graph g into Knossos Expr.
    """

    # names_to_value_infos = {x.name:x.type for x in g.value_info}

    schemas = get_all_schemas()

    # Innermost body is "outputs"
    output_vars = [useVar(o.name) for o in g.output]
    output_types = tuple(onnxType_to_Type(o.type) for o in g.output)
    if len(g.output) == 1:
        body = output_vars[0]
        body_type = output_types[0]
    else:
        body = Call("tuple", output_vars)
        body_type = Type.Tuple(*output_types)

    # ONNX graph nodes are topologically sorted, so just run through in reverse order
    # making a chain of lets.
    for node in reversed(g.node):
        opname = node.op_type if not node.domain else node.domain + "." + node.op_type
        found_schema = opname in schemas
        if found_schema:
            schema = schemas[opname]
        else:
            # Not found, but let's try to proceed, consing up a
            # fake object with the bare minimum number of fields.
            warnings.warn(f"onnx2ks: Op {opname} not found -- just making a call")
            schema = SimpleNamespace()
            schema.name = opname
            schema.attributes = {}

        name = opname

        # Collect args from input
        args = [useVar(i) for i in node.input]

        # Special cases
        if opname == "Constant":
            # Constant: exactly one of the optionals should be set
            assert len(node.attribute) == 1
            n = node.attribute[0]
            args = [Expr_from_Attr(n)]

        elif opname == "Cast":
            # Cast: output type depends on input value.  We postpend the type to the name
            assert len(node.attribute) == 1
            n = node.attribute[0]
            assert n.type == OpSchema.AttrType.INT
            out_type = OpSchema.AttrType(n.i)
            out_type = str(out_type)
            assert out_type.startswith("AttrType.")
            out_type = out_type[9:]
            name = "Cast_" + out_type

        else:
            # Collect attributes.
            # - Some are set on the node.
            node_attrs = dict()
            for n in node.attribute:
                val = Expr_from_Attr(n)
                node_attrs[n.name] = val
                if n.name not in schema.attributes:
                    warnings.warn(
                        f"Attribute {n.name} not in schema for {opname} -- adding arg anyway"
                    )
                    args += [val]

            # - The rest are defaulted
            for attrname, attr in schema.attributes.items():
                if attrname in node_attrs:
                    args += [node_attrs[attrname]]
                else:
                    args += [get_default_value(schema, attr)]

        # Are we dropping optional outputs?
        # This is typically only known at the call site.
        # Change the function name to reflect that,
        # otherwise type annotation would have to happen
        # at more than one place
        n_outputs = len(node.output)
        if found_schema and n_outputs < len(schema.outputs):
            name = f"take{n_outputs}${name}"

        # Generate the call node
        rhs = Call(name, args)

        # Now bind multiple outputs
        # (let ((out1 ... outN) (op args))
        #   body)
        vars = [useVar(o) for o in node.output]
        if len(vars) == 1:
            vars = vars[0]

        # Rest of the time, wrap previous body in a Let
        body = Let(vars, rhs, body)

    body = Let(Var("$end_of_inits #|End of initializers|# "), Const(99999), body)

    # And now gather the initializers, and place in Lets
    inputs = set([i.name for i in g.input])

    inits = {init.name: init for init in g.initializer}

    internal_inits = []
    for init in g.initializer:
        if init.name not in inputs:
            internal_inits.append(init)

    # value_infos -> types
    args = [defVar(i.name, onnxType_to_Type(i.type)) for i in g.input]

    body = Let(
        Var("$beg_of_internal_inits #|Begin of internal initializers|# "),
        Const(99999),
        body,
    )
    ex = emit_inits(internal_inits, body)

    # Emit a def for the whole graph
    #    (def GraphName None args body)
    decl = Def(g.name, body_type, args, ex)

    # And emit a "main" to call it with initializers
    # (def main None ()
    #     (GraphName initializers.. ) )
    main_body = Call(g.name, [useVar(i.name) for i in g.input])
    input_inits = [inits[i.name] for i in g.input if i.name in inits]
    main_args = [arg for arg in args if arg.name not in inits]
    main = Def("main", body_type, main_args, emit_inits(input_inits, main_body))

    # Collect the two decls
    return [decl, main]

if __name__ == "__main__":
    import os
    from ksc.parse_ks import parse_ks_filename
    from ksc.type_propagate import type_propagate_decls

    argv = sys.argv

    if len(argv) == 1:
        os.makedirs("obj/test/onnx2ks", exist_ok=True)
        argv = [
            "*EXE*",
            "test/onnx2ks/mnist-7.onnx",
            #'test/onnx2k/bert-tiny-uncased_L_3_H_128_A_2_V_30528_S_512_Dp_0.1_optimized_layer_norm_opset12.onnx',
            "obj/",
        ]

    nargs = len(argv) - 1
    if nargs == 0:
        filename = "/dev/stdin"
    else:
        filename = argv[1]

    outbase = None
    if nargs == 2:
        outbase = argv[2]
        if outbase.endswith("/"):
            outbase = outbase + re.sub(r"\.onnx$", "", filename)

    type_propagate = False  # Let's use onnx shape inference

    def save(decls, ofn, msg):
        print(f"onnx2ks: Writing to {ofn}")
        with open(ofn, "w") as out:
            out.write(f";; {msg}\n")
            for decl in decls:
                pprint(decl, width=256, ribbon_width=256, stream=out)
                out.write("\n")

    # Load our file
    print(f"onnx2ks: Reading from {filename}")
    model = onnx.load(filename)

    # Make output dir
    odir = os.path.dirname(outbase)
    print(f"onnx2ks: Making output dir {odir}")
    os.makedirs(odir, exist_ok=True)

    print(f"onnx2ks: Writing graph to {outbase}.onnx.txt")
    with open(outbase + ".onnx.txt", "w") as out:
        out.write(onnx.helper.printable_graph(model.graph))

    # Convert to untyped KS
    model = onnx.shape_inference.infer_shapes(model)

    if False:
        # Trying to fix non-topological sort from
        onnx.optimizer.optimize(model, ["extract_constant_to_initializer"])
        onnx.checker.check_model(model)
    decls = onnx2ks(model.graph)

    if outbase:
        save(decls, outbase + ".untyped.ks", f" untyped KS from {filename}")

    # Check model -- do ORT models need to be toposorted?
    # onnx.checker.check_model(model)

    if type_propagate:
        # Load preludes
        print(f"onnx2ks: Reading prelude")
        symtab = dict()

        prelude_decls = parse_ks_filename("etc/onnx-prelude.ks")
        type_propagate_decls(prelude_decls, symtab)

        prelude_decls = parse_ks_filename("etc/onnx-prelude-autogen.ks")
        type_propagate_decls(prelude_decls, symtab)

        # Apply type annotations
        type_propagate_decls(decls, symtab)

    if outbase:
        # And save
        save(decls, outbase + ".ks", f"AUTOGENERATED FROM {filename}")

        # And load again...
        decls = parse_ks_filename(outbase + ".ks")

    if False:
        # Pretty print
        for decl in decls:
            cpprint(decl, width=132, ribbon_width=132)
            print("")

# %%

max_stars_repo_name: tomjaguarpaw/knossos-ksc

max_issues_repo_name: tomjaguarpaw/knossos-ksc

max_forks_repo_name: tomjaguarpaw/knossos-ksc


text:
########### test in sep on FPN ###########
import numpy as np
import torch
torch.manual_seed(2019)
np.random.seed(2019)

import torch.nn as nn
import torch.nn.parallel
from torch.nn.parameter import Parameter
from torch.autograd import Variable
from torch.nn import functional as F
import torch.nn.functional as F
from config import BATCHNORM_MOMENTUM, BOX_SCALE, IM_SCALE, ModelConfig
from lib.fpn.nms.functions.nms import apply_nms

from lib.fpn.box_utils import bbox_overlaps, center_size, bbox_preds, nms_overlaps
from lib.get_union_boxes import NewUnionBoxesAndFeats, UnionBoxesAndFeats
from lib.fpn.proposal_assignments.new_rel_assignments import rel_assignments
from lib.object_detector import ObjectDetector, gather_res, load_vgg, Result
from lib.pytorch_misc import transpose_packed_sequence_inds, to_onehot, arange, enumerate_by_image, diagonal_inds, Flattener, get_dropout_mask
from lib.sparse_targets import FrequencyBias
from lib.surgery import filter_dets
from lib.word_vectors import obj_edge_vectors
from lib.fpn.roi_align.functions.roi_align import RoIAlignFunction
from lib.lstm.highway_lstm_cuda.alternating_highway_lstm import block_orthogonal
import time
from config import IM_SCALE, ROOT_PATH, CO_OCCOUR_PATH
#from lib.ggnn import GGNNObj
#from torch.distributions.uniform import Uniform
#import torch.distributions

MODES = ('sgdet', 'sgcls', 'predcls')
conf = ModelConfig()

def myNNLinear(input_dim, output_dim, bias=True):
    ret_layer = nn.Linear(input_dim, output_dim, bias=bias)
    ret_layer.weight = torch.nn.init.xavier_normal_(ret_layer.weight, gain=1.0)
    return ret_layer

class DynamicFilterContext(nn.Module):

    def __init__(self, classes, rel_classes, mode='sgdet', use_vision=True,
                 embed_dim=200, hidden_dim=512, obj_dim=4096, pooling_dim=4096,
                 pooling_size=7, dropout_rate=0.2, use_bias=True, use_tanh=True,
                 limit_vision=True, sl_pretrain=False, num_iter=-1, use_resnet=False,
                 reduce_input=False, debug_type=None, post_nms_thresh=0.5,
                 output_dim=512,):
##               time_step_num=3, use_knowledge=True, knowledge_matrix=''):

        super(DynamicFilterContext, self).__init__()
        self.classes = classes
        self.rel_classes = rel_classes
        assert mode in MODES
        self.mode = mode

        self.use_vision = use_vision
        self.use_bias = use_bias
        self.use_tanh = use_tanh
        self.use_highway = True
        self.limit_vision = limit_vision
        self.pooling_dim = pooling_dim
        self.pooling_size = pooling_size
        self.nms_thresh = post_nms_thresh

        self.obj_compress = myNNLinear(self.pooling_dim, self.num_classes, bias=True)

        self.roi_fmap_obj = load_vgg(pretrained=False).classifier

        if self.use_bias:
            self.freq_bias = FrequencyBias()
        self.reduce_dim = 256
        self.reduce_obj_fmaps = nn.Conv2d(512, self.reduce_dim, kernel_size=1)

        similar_fun = [myNNLinear(self.reduce_dim, self.reduce_dim),
##        similar_fun = [myNNLinear(self.reduce_dim*2, self.reduce_dim),
                       nn.ReLU(inplace=True),
                       nn.Dropout(p=0.5),
#                       nn.BatchNorm2d(num_features=self.reduce_dim),
#                      myNNLinear(self.reduce_dim, self.reduce_dim),
             #          nn.ReLU(inplace=True),
                       myNNLinear(self.reduce_dim, 1)]
        self.similar_fun = nn.Sequential(*similar_fun)

        roi_fmap = [Flattener(),
                    nn.Linear(self.reduce_dim*2*self.pooling_size*self.pooling_size, 4096, bias=True),
                    nn.ReLU(inplace=True),
                    nn.Dropout(p=0.5),
                    nn.Linear(4096, 4096, bias=True)]
        self.roi_fmap = nn.Sequential(*roi_fmap)

        self.hidden_dim = hidden_dim
        self.rel_compress = myNNLinear(self.hidden_dim*3, self.num_rels)
##        self.rel_compress = myNNLinear(self.pooling_dim, self.num_rels, bias=True)
        self.post_obj = myNNLinear(self.pooling_dim, self.hidden_dim*2)
        self.mapping_x = myNNLinear(self.hidden_dim*2, self.hidden_dim*3)
        self.reduce_rel_input = myNNLinear(self.pooling_dim, self.hidden_dim*3)

        def obj_feature_map(self, features, rois):
        feature_pool = RoIAlignFunction(self.pooling_size, self.pooling_size, spatial_scale=1 / 16)(
            features, rois)
        feature_pool2 = RoIAlignFunction(self.pooling_size-2, self.pooling_size-2, spatial_scale=1 / 16)(
            features, rois)
        feature_pool3 = RoIAlignFunction(1, 1, spatial_scale=1 / 16)(
            features, rois)
###        return feature_pool
        return feature_pool, feature_pool2, feature_pool3
       # return self.roi_fmap_obj(feature_pool.view(rois.size(0), -1))

    @property
    def num_classes(self):
        return len(self.classes)
    @property
    def num_rels(self):
        return len(self.rel_classes)

    @property
    def is_sgdet(self):
        return self.mode == 'sgdet'

    @property
    def is_sgcls(self):
        return self.mode == 'sgcls'

    def forward(self, *args, **kwargs):

        results = self.base_forward(*args, **kwargs)
        return results

    def base_forward(self, fmaps, obj_logits, im_inds, rel_inds, msg_rel_inds, reward_rel_inds, im_sizes, boxes_priors=None, boxes_deltas=None,$
        assert self.mode == 'sgcls'
        num_objs = obj_logits.shape[0]
#       rel_inds = rel_labels[:, :3].data.clone()
        num_rels = rel_inds.shape[0]  ##506
        rois = torch.cat((im_inds[:, None].float(), boxes_priors), 1)  #[23, 5]
##################################################################
        obj_fmaps = self.obj_feature_map(fmaps, rois)  ##[23, 512, 7, 7]

##############################################################################################################################
###############################################################################################################################
        reduce_obj_fmaps_0 = torch.nn.functional.interpolate(self.reduce_obj_fmaps(obj_fmaps[0]), size=(4,4), scale_factor=None, mode='nearest'$
##        reduce_obj_fmaps_1 = self.reduce_obj_fmaps(obj_fmaps[2])  ##[23, 256, 5, 5]
        reduce_obj_fmaps_1 = torch.nn.functional.interpolate(self.reduce_obj_fmaps(obj_fmaps[1]), size=(4,4), scale_factor=None, mode='nearest'$
        reduce_obj_fmaps_2 = self.reduce_obj_fmaps(obj_fmaps[2])  ##[23, 256, 1, 1]

        S_fmaps_0 = reduce_obj_fmaps_0[rel_inds[:, 1]]  ##[506, 256, 4, 4]
        O_fmaps_0 = reduce_obj_fmaps_0[rel_inds[:, 2]]

        S_fmaps_1 = reduce_obj_fmaps_1[rel_inds[:, 1]]  ##[506, 256, 4, 4]
        O_fmaps_1 = reduce_obj_fmaps_1[rel_inds[:, 2]]

        S_fmaps_2 = reduce_obj_fmaps_2[rel_inds[:, 1]]  ##[506, 256, 1, 1]
        O_fmaps_2 = reduce_obj_fmaps_2[rel_inds[:, 2]]

        pooling_size_sq = 16

        S_fmaps_trans_0 = S_fmaps_0.view(num_rels, self.reduce_dim, pooling_size_sq).transpose(2, 1)
        O_fmaps_trans_0 = O_fmaps_0.view(num_rels, self.reduce_dim, pooling_size_sq).transpose(2, 1)

        S_fmaps_trans_1 = S_fmaps_1.view(num_rels, self.reduce_dim, pooling_size_sq).transpose(2, 1)
        O_fmaps_trans_1 = O_fmaps_1.view(num_rels, self.reduce_dim, pooling_size_sq).transpose(2, 1)

        S_fmaps_trans_2 = S_fmaps_2.view(num_rels, self.reduce_dim, 1).transpose(2, 1)
        O_fmaps_trans_2 = O_fmaps_2.view(num_rels, self.reduce_dim, 1).transpose(2, 1)

                SO_fmaps_extend_0 = (S_fmaps_trans_0.unsqueeze(1).expand(-1, pooling_size_sq, -1, -1) + O_fmaps_trans_0.unsqueeze(2).expand(-1, -1, pooling_size_sq, -1)).view(num_rels, pooling_size_sq*pooling_size_sq, self.reduce_dim)
        SO_fmaps_extend_1 = (S_fmaps_trans_1.unsqueeze(1).expand(-1, pooling_size_sq, -1, -1) + O_fmaps_trans_1.unsqueeze(2).expand(-1, -1, pooling_size_sq, -1)).view(num_rels, pooling_size_sq*pooling_size_sq, self.reduce_dim)
        SO_fmaps_extend_2 = (S_fmaps_trans_2.unsqueeze(1).expand(-1, 1, -1, -1) + O_fmaps_trans_2.unsqueeze(2).expand(-1, -1, 1, -1)).view(num_rels, 1, self.reduce_dim)

        presence_vector_0 =  torch.sigmoid(self.similar_fun(SO_fmaps_extend_0)).view(num_rels, pooling_size_sq, pooling_size_sq)
        presence_vector_1 =  torch.sigmoid(self.similar_fun(SO_fmaps_extend_1)).view(num_rels, pooling_size_sq, pooling_size_sq)
        presence_vector_2 =  torch.sigmoid(self.similar_fun(SO_fmaps_extend_2))

        m_vector_0 = torch.mean(SO_fmaps_extend_0.transpose(2,1), dim=1).view(-1)  ##129536
        m_vector_1 = torch.mean(SO_fmaps_extend_1.transpose(2,1), dim=1).view(-1)
        m_vector_2 = torch.mean(SO_fmaps_extend_2.transpose(2,1), dim=1).view(-1) ##506

        latent_mask_0 = torch.where(F.softmax(presence_vector_0.view(-1), dim=0) >= (1/len(m_vector_0)), torch.ones_like(presence_vector_0.view(-1)), torch.zeros_like(presence_vector_0.view(-1)))
        latent_mask_1 = torch.where(F.softmax(presence_vector_1.view(-1), dim=0) >= (1/len(m_vector_1)), torch.ones_like(presence_vector_1.view(-1)), torch.zeros_like(presence_vector_1.view(-1)))
        latent_mask_2 = torch.where(F.softmax(presence_vector_2.view(-1), dim=0) >= (1/len(m_vector_2)), torch.ones_like(presence_vector_2.view(-1)), torch.zeros_like(presence_vector_2.view(-1)))

        attended_vector_0 =  (m_vector_0 * latent_mask_0).view(num_rels, pooling_size_sq, pooling_size_sq)
        attended_vector_1 =  (m_vector_1 * latent_mask_1).view(num_rels, pooling_size_sq, pooling_size_sq)
        attended_vector_2 =  (m_vector_2 * latent_mask_2).view(num_rels, 1, 1)

        SO_fmaps_scores_0 = F.softmax(attended_vector_0, dim=1)
        SO_fmaps_scores_1 = F.softmax(attended_vector_1, dim=1)
        SO_fmaps_scores_2 = F.softmax(attended_vector_2, dim=1)

        SO_fmaps_scores = SO_fmaps_scores_0*SO_fmaps_scores_1
        SO_fmaps_scores = torch.matmul(SO_fmaps_scores.view(num_rels,1,pooling_size_sq*pooling_size_sq).transpose(2, 1) , SO_fmaps_scores_2)  ##[506, 256, 1]

        weighted_S_fmaps = torch.matmul(SO_fmaps_scores.view(num_rels,pooling_size_sq,pooling_size_sq).transpose(2, 1), S_fmaps_trans_0)  ##[506, 16, 256]
        last_SO_fmaps = torch.cat((weighted_S_fmaps, O_fmaps_trans_0), dim=2)
        last_SO_fmaps = nn.functional.interpolate(last_SO_fmaps.transpose(2, 1).view(num_rels,self.reduce_dim*2,4,4), size=(7,7))

        obj_feats = self.roi_fmap_obj(obj_fmaps[0].view(rois.size(0), -1))

        obj_logits = self.obj_compress(obj_feats)  ##[23, 151]
        obj_dists = F.softmax(obj_logits, dim=1)  ##[23, 151]

        pred_obj_cls = obj_dists[:, 1:].max(1)[1] + 1    ##[23]
#        sub_dist = obj_distributions.view(1, obj_num, num_class).expand(obj_num, obj_num, num_class).contiguous().view(-1, num_class)
        # for relationship classification
        rel_input = self.roi_fmap(last_SO_fmaps)  ##[506, 4096]
        subobj_rep = self.post_obj(obj_feats)  ##[23, 1024]
        sub_rep = subobj_rep[:, :self.hidden_dim][rel_inds[:, 1]]  ##[506, 512]
        obj_rep = subobj_rep[:, self.hidden_dim:][rel_inds[:, 2]]

        last_rel_input = self.reduce_rel_input(rel_input)  ##[506, 1536]
        last_obj_input = self.mapping_x(torch.cat((sub_rep, obj_rep), 1))
        triple_rep = nn.ReLU(inplace=True)(last_obj_input + last_rel_input) - (last_obj_input - last_rel_input).pow(2) ##[506, 1536]
        rel_logits = self.rel_compress(triple_rep)

        # follow neural-motifs paper
        if self.use_bias:
            if self.mode in ['sgcls', 'sgdet']:
                rel_logits = rel_logits + self.freq_bias.index_with_labels(
                    torch.stack((
                        pred_obj_cls[rel_inds[:, 1]],  ##result.obj_preds[rel_inds[:, 1]],
                        pred_obj_cls[rel_inds[:, 2]],
                        ), 1))
            elif self.mode == 'predcls':
                rel_logits = rel_logits + self.freq_bias.index_with_labels(
                torch.stack((
                        obj_labels[rel_inds[:, 1]],
                        obj_labels[rel_inds[:, 2]],
                        ), 1))
            else:
                raise NotImplementedError

        return pred_obj_cls, obj_logits, rel_logits

class RelModelAlign(nn.Module):

    def __init__(self, classes, rel_classes, mode='sgdet', num_gpus=1, use_vision=True, require_overlap_det=True,
                 embed_dim=200, hidden_dim=512, pooling_dim=2048, use_resnet=False, thresh=0.01,
                 use_proposals=False, rec_dropout=0.2, use_bias=True, use_tanh=True,
                 limit_vision=True, sl_pretrain=False, eval_rel_objs=False, num_iter=-1, reduce_input=False,
                 post_nms_thresh=0.5,):

        super(RelModelAlign, self).__init__()
        self.classes = classes
        self.rel_classes = rel_classes
        self.num_gpus = num_gpus
        assert mode in MODES
        self.mode = mode

        self.pooling_size = 7
        ##self.pooling_size = conf.pooling_size
        self.embed_dim = embed_dim
        self.hidden_dim = hidden_dim
        self.obj_dim = 2048 if use_resnet else 4096
        self.pooling_dim = pooling_dim

        self.use_bias = use_bias
        self.use_vision = use_vision
        self.use_tanh = use_tanh
        self.limit_vision = limit_vision
        self.num_iter = num_iter
        self.require_overlap = require_overlap_det and self.mode == 'sgdet'
        self.sl_pretrain = sl_pretrain
        ################################### Added  ################################################
        self.co_occour = np.load(CO_OCCOUR_PATH)
        self.co_occour = self.co_occour / self.co_occour.sum()
        ###########################################################################################
        self.detector = ObjectDetector(
            classes=classes,
            mode=('proposals' if use_proposals else 'refinerels') if mode == 'sgdet' else 'gtbox',
            use_resnet=use_resnet,
            thresh=thresh,
            max_per_img=64,
        )

        self.context = DynamicFilterContext(self.classes, self.rel_classes, mode=self.mode,
                                            use_vision=self.use_vision, embed_dim=self.embed_dim,
                                            hidden_dim=self.hidden_dim, obj_dim=self.obj_dim,
                                            pooling_dim=self.pooling_dim, pooling_size=self.pooling_size,
                                            dropout_rate=rec_dropout,
                                            use_bias=self.use_bias, use_tanh=self.use_tanh,
                                            limit_vision=self.limit_vision,
                                            sl_pretrain = self.sl_pretrain,
                                            num_iter=self.num_iter,
                                            use_resnet=use_resnet,
                                            reduce_input=reduce_input,
                                            post_nms_thresh=post_nms_thresh,)
        @property
        def num_classes(self):
            return len(self.classes)
       @property
        def num_rels(self):
            return len(self.rel_classes)

    def get_reward_rel_inds(self, im_inds, box_priors, box_score):    ##get_rel_inds(self, rel_labels, im_inds, box_priores)

        rel_cands = im_inds.data[:, None] == im_inds.data[None]
        rel_cands.view(-1)[diagonal_inds(rel_cands)] = 0

        if self.require_overlap:
            rel_cands = rel_cands & (bbox_overlaps(box_priors.data,
                                                   box_priors.data) > 0)
        rel_cands = rel_cands.nonzero()
        if rel_cands.dim() == 0:
            rel_cands = im_inds.data.new(1, 2).fill_(0)
        rel_cands = rel_cands.nonzero()
        if rel_cands.dim() == 0:
            rel_cands = im_inds.data.new(1, 2).fill_(0)

        rel_inds = torch.cat((im_inds.data[rel_cands[:, 0]][:, None], rel_cands), 1)
        return rel_inds

    def get_msg_rel_inds(self, im_inds, box_priors, box_score):

        rel_cands = im_inds.data[:, None] == im_inds.data[None]
        rel_cands.view(-1)[diagonal_inds(rel_cands)] = 0

        if self.require_overlap:
            rel_cands = rel_cands & (bbox_overlaps(box_priors.data,
                                                   box_priors.data) > conf.overlap_thresh)
        rel_cands = rel_cands.nonzero()
        if rel_cands.dim() == 0:
            rel_cands = im_inds.data.new(1, 2).fill_(0)

        rel_inds = torch.cat((im_inds.data[rel_cands[:, 0]][:, None], rel_cands), 1)
        return rel_inds

    def get_rel_inds(self, rel_labels, im_inds, box_priors, box_score):

        if self.training:
            rel_inds = rel_labels[:, :3].data.clone()
        else:
            rel_cands = im_inds.data[:, None] == im_inds.data[None]
            rel_cands.view(-1)[diagonal_inds(rel_cands)] = 0

            # Require overlap for detection
            # Require overlap in the test stage
            if self.require_overlap:
                rel_cands = rel_cands & (bbox_overlaps(box_priors.data,
                                                       box_priors.data) > 0)
            rel_cands = rel_cands.nonzero()
            if rel_cands.dim() == 0:
                rel_cands = im_inds.data.new(1, 2).fill_(0)
            rel_inds = torch.cat((im_inds.data[rel_cands[:, 0]][:, None], rel_cands), 1)
        return rel_inds

    def forward(self, x, im_sizes, image_offset,
                gt_boxes=None, gt_classes=None, gt_rels=None, proposals=None, train_anchor_inds=None,
                return_fmap=False):

        result = self.detector(x, im_sizes, image_offset, gt_boxes, gt_classes, gt_rels, proposals,
                               train_anchor_inds, return_fmap=True)

        if result.is_none():
            return ValueError("heck")

        im_inds = result.im_inds - image_offset
        boxes = result.rm_box_priors
        # boxes = result.boxes_assigned
        boxes_deltas = result.rm_box_deltas # sgcls is None
        boxes_all = result.boxes_all # sgcls is None
        if (self.training) and (result.rel_labels is None):
            import pdb; pdb.set_trace()
            print('debug')
            assert self.mode == 'sgdet'
            result.rel_labels = rel_assignments(im_inds.data, boxes.data, result.rm_obj_labels.data, result.rm_obj_dists.data,
                                                gt_boxes.data, gt_classes.data, gt_rels.data,
                                                image_offset, filter_non_overlap=True,
                                                num_sample_per_gt=1)

        rel_inds = self.get_rel_inds(result.rel_labels, im_inds, boxes, result.rm_obj_dists.data)
        ############################## change to ##########################################
##        rel_inds = self.get_reward_rel_inds(result.rel_labels, im_inds, boxes, result.rm_obj_dists.data)

        reward_rel_inds = None
        #########################################
        if self.mode == 'sgdet':
            msg_rel_inds = self.get_msg_rel_inds(im_inds, boxes, result.rm_obj_dists.data)
            reward_rel_inds = self.get_reward_rel_inds(im_inds, boxes, result.rm_obj_dists.data)

        if self.mode == 'sgdet':
            result.rm_obj_dists_list, result.obj_preds_list, result.rel_dists_list, result.bbox_list, result.offset_list, \
                result.rel_dists, result.obj_preds, result.boxes_all, result.all_rel_logits = self.context(
                                            result.fmap.detach(), result.rm_obj_dists.detach(), im_inds, rel_inds, msg_rel_inds,
                                            reward_rel_inds, im_sizes, boxes.detach(), boxes_deltas.detach(), boxes_all.detach(),
                                            result.rm_obj_labels if self.training or self.mode == 'predcls' else None)

        elif self.mode in ['sgcls', 'predcls']:
           # import pdb; pdb.set_trace()
            result.obj_preds, result.rm_obj_logits, result.rel_logits = self.context(
                                            result.fmap.detach(), result.rm_obj_dists.detach(),
                                            im_inds, rel_inds, None, None, im_sizes, boxes.detach(), None, None, self.co_occour,
                                            result.rm_obj_labels if self.training or self.mode == 'predcls' else None)

        else:
            raise NotImplementedError

        # result.rm_obj_dists = result.rm_obj_dists_list[-1]

        if self.training:
            return result

        if self.mode == 'predcls':
            import pdb; pdb.set_trace()
            print('debug..')
            result.obj_preds = result.rm_obj_labels
            result.obj_scores = Variable(torch.from_numpy(np.ones(result.obj_preds.shape[0],)).float().cuda())
        else:
            twod_inds = arange(result.obj_preds.data) * self.num_classes + result.obj_preds.data
            result.obj_scores = F.softmax(result.rm_obj_logits, dim=1).view(-1)[twod_inds]

        # # Bbox regression
        if self.mode == 'sgdet':
            if conf.use_postprocess:
                bboxes = result.boxes_all.view(-1, 4)[twod_inds].view(result.boxes_all.size(0), 4)
            else:
                bboxes = result.rm_box_priors
        else:
            bboxes = result.rm_box_priors

        rel_scores = F.softmax(result.rel_logits, dim=1)
##        return filter_dets(bboxes, result.obj_scores,
##                           result.obj_preds, rel_inds[:, 1:], rel_scores)
       #################################################################################################
        return filter_dets(bboxes, result.obj_scores,
                           result.obj_preds, rel_inds[:, 1:], rel_scores, gt_boxes, gt_classes, gt_rels)
       ##################################################################################################
    def __getitem__(self, batch):
        """ Hack to do multi-GPU training"""

        batch.scatter()
        if self.num_gpus == 1:
            return self(*batch[0])

        replicas = nn.parallel.replicate(self, devices=list(range(self.num_gpus)))
        outputs = nn.parallel.parallel_apply(replicas, [batch[i] for i in range(self.num_gpus)])

        if self.training:
            return gather_res(outputs, 0, dim=0)
        return outputs


text:
import numpy as np
from scipy.special import erf
import math
from scipy.integrate import quad
import resfuncRead as rfr
import scipy.optimize as so
import pandas as pds
import lmfit as lmf
import nrfano_stats as nfs
import copy

#bootstrapping
import bootstrapped.bootstrap as bs
import bootstrapped.stats_functions as bs_stats

def QEr_Ebin(Q, Ernr, bins=[5, 10, 20, 30, 40, 50, 70,150],silent=False):

    #create a dataframe
    nr_df = pds.DataFrame(data={'yield':Q, 'energy':Ernr})

    #bin the data
    nr_df['binned'] = pds.cut(nr_df['energy'],bins)

    #print stats in each bin
    s = nr_df.groupby(pds.cut(nr_df['energy'], bins=bins)).size()
    if not silent:
      print (s)

    #create list of vectors for histogrammin'
    bindf = nr_df.groupby(pds.cut(nr_df['energy'], bins=bins))['yield'].apply(list)
    bindfE = nr_df.groupby(pds.cut(nr_df['energy'], bins=bins))['energy'].apply(list)
    #print(bindf)

    return bindf,bindfE

def QEr_Qboot(bindf,bins=[5, 10, 20, 30, 40, 50, 70,150],silent=False):

    qbootsigs = np.zeros((np.shape(bins)[0]-1,))
    qbootsigerrsu = np.zeros((np.shape(bins)[0]-1,))
    qbootsigerrsl = np.zeros((np.shape(bins)[0]-1,))

    for i,Qv in enumerate(bindf):
      if not silent:
        print(np.shape(Qv))
      Qv = np.asarray(Qv)
      #print(Qv[0:10])
      try:
        bsr = bs.bootstrap(Qv, stat_func=bs_stats.std,iteration_batch_size=100)
      except MemoryError as e:
        print('There was a memory error - too much memory to be allocated')

      if not silent:
        print(bsr)
      qbootsigs[i] = np.std(Qv)
      qbootsigerrsu[i] = bsr.upper_bound
      qbootsigerrsl[i] = bsr.lower_bound

    #change over to size of error bars, not confidence interval
    qbootsigerrsu = qbootsigerrsu - qbootsigs
    qbootsigerrsl = -qbootsigerrsl + qbootsigs

    return qbootsigs,qbootsigerrsl,qbootsigerrsu

def QEr_QbootBC(bindf,qbootsigs,qbootEs,n=10,bins=[5, 10, 20, 30, 40, 50, 70,150],silent=False):

    qbootcorrs = np.ones((np.shape(bins)[0]-1,))

    for i,Ev in enumerate(bindf):
        if((i>0)&(i<(np.shape(qbootsigs)[0]-1))):
          m = (qbootsigs[i+1] -qbootsigs[i-1])/(qbootEs[i+1]-qbootEs[i-1])
        elif (i>0):
          m = (qbootsigs[i] -qbootsigs[i-1])/(qbootEs[i]-qbootEs[i-1])
        elif (i<(np.shape(qbootsigs)[0]-1)):
          m = (qbootsigs[i+1] -qbootsigs[i])/(qbootEs[i+1]-qbootEs[i])
        intercept = qbootsigs[i]
        #print(qbootEs[i])
        #print(qbootsigs[i])
        fsig = lambda E: m*(E-qbootEs[i]) + intercept
        #print(fsig(qbootEs[i]))
        sigcorr = bc_corr(Ev,fsig,n)
        #print(qbootsigs[i]/sigcorr)
        qbootcorrs[i] = (qbootsigs[i]/sigcorr)

    #first two are absurd because of negative projected sigma (FIXME)
    #qbootcorrs[0] = 1
    #qbootcorrs[1] = 1

    if not silent:
      print(qbootcorrs)
    return qbootcorrs

def QEr_QbootBC_iterative(bindf,qbootsigs,qbootEs,n=10,bins=[5, 10, 20, 30, 40, 50, 70,150],silent=False):

    qbootsigs_cp = copy.deepcopy(qbootsigs)
    qbootEs_cp = copy.deepcopy(qbootEs)

    for i in np.arange(5):
      qbootcorrs = QEr_QbootBC(bindf,qbootsigs,qbootEs,n,bins,silent)
      qbootsigs = qbootsigs*qbootcorrs

    print(qbootsigs/qbootsigs_cp)
    return qbootsigs/qbootsigs_cp

def QEr_Qhist(bindf, qbins=np.linspace(0,0.6,40)):

    xcq = (qbins[:-1] + qbins[1:]) / 2

    qhistos = np.zeros((np.shape(qbins)[0]-1,0))
    qerrs = np.zeros((np.shape(qbins)[0]-1,0))

    #get errors for 0-20 here
    fcerrs = nfs.largestErr_fast()

    for i,Qv in enumerate(bindf):
      n,nx = np.histogram(Qv,bins=qbins)
      n = np.reshape(n,(np.shape(n)[0],1))
      qhistos = np.append(qhistos,n,axis=1)
      qerrs0 = np.sqrt(n)
      qerrs0[n<=20] = fcerrs[n[n<=20]]
      qerrs = np.append(qerrs,qerrs0,axis=1)
      #qerrs[n<=20] = fcerrs[n[n<=20]]
      #use gaussian errors
      #qerrs = np.append(qerrs,np.sqrt(n),axis=1)
      #qerrs[qerrs==0]=1
      #use FC poissonian errors
      #qerrs = np.append(qerrs,nfs.largestErr(n),axis=1)

    return qhistos,qerrs

def QEr_Qfit(qhistos,qerrs, qbins=np.linspace(0,0.6,40),damps=0.1,dmu=1.0,dsig=0.1,silent=False):

    xcq = (qbins[:-1] + qbins[1:]) / 2

    qamps = np.zeros((np.shape(qhistos)[1],))
    qampserrs = np.zeros((np.shape(qhistos)[1],))
    qmus = np.zeros((np.shape(qhistos)[1],))
    qmuerrs = np.zeros((np.shape(qhistos)[1],))
    qsigs = np.zeros((np.shape(qhistos)[1],))
    qsigerrs = np.zeros((np.shape(qhistos)[1],))

    startamps = damps*np.ones((np.shape(qhistos)[1],))
    startmus = dmu*np.ones((np.shape(qhistos)[1],))
    startsigs = dsig*np.ones((np.shape(qhistos)[1],))

    for i,h in enumerate(qhistos[0,:]):
      if not silent:
        print('fitting {}'.format(i))

      qsum = np.sum(qhistos[:,i])
      #do it with lmfit
      params = lmf.Parameters()
      params.add('amp', value=startamps[i])
      params.add('mean', value=startmus[i])
      params.add('sig', value=startsigs[i])
      lmfout = lmf.minimize(gauss_residual, params, args=(xcq, qhistos[:,i]/qsum, qerrs[:,i]/qsum))
      #print(lmf.fit_report(lmfout))
      if not silent:
        print('lmfit results')
        print(lmf.report_fit(lmfout.params))
      qamps[i] = lmfout.params['amp'].value
      qmus[i] = lmfout.params['mean'].value
      qsigs[i] = lmfout.params['sig'].value

      #somtimes covariance doesn't exist (if bad fit)
      try:
        qampserrs[i] = np.sqrt(lmfout.covar[0,0])
        qmuerrs[i] = np.sqrt(lmfout.covar[1,1])
        qsigerrs[i] = np.sqrt(lmfout.covar[2,2])
      except:
        print('bad fit')
        qampserrs[i] = -1
        qmuerrs[i] = -1
        qsigerrs[i] = -1

    if not silent:
      print(qsigs)
      print(qsigerrs)

    return qamps,qampserrs,qmus,qmuerrs,qsigs,qsigerrs

def gauss_residual(params, x, data, eps_data):
    amp = params['amp']
    mean = params['mean']
    sig = params['sig']

    model = amp * np.exp(-(x-mean)**2/(2*sig**2))

    return (data-model) / eps_data

#a bin-centering correction
def bc_corr(E,sig,n=10,Qbar=lambda E: 0.16*E**(0.18)):

    Ev = np.linspace(np.amin(E),np.amax(E),n)
    Ec = (Ev[:-1] + Ev[1:]) / 2

    #create a dataframe
    df = pds.DataFrame(data={'energy':E})

    s = df.groupby(pds.cut(df['energy'], bins=Ev)).size()[:]
    s = np.asarray(s)
    s = s/np.sum(s)
    #print(s)

    #make a function
    func = {}
    for i,Er in enumerate(Ec):
      idx = 'f{}'.format(i)
      idx = str(idx)
      prefac = copy.copy(s[i])
      #func[idx] = lambda Q,prefac=prefac,Er=Er: prefac*(1/np.sqrt(2*np.pi*sig(Er)**2))*np.exp(-(Qbar(Er)-Q)**2/(2*sig(Er)**2))
      func[idx] = (lambda prefac,Er: (lambda Q: prefac*(1/np.sqrt(2*np.pi*sig(Er)**2))*np.exp(-(Qbar(Er)-Q)**2/(2*sig(Er)**2))))(prefac,Er)

    fmean = lambda Q: Q*fsum(Q,func)
    f2mom = lambda Q: Q**2*fsum(Q,func)

    #print(quad(fsum,-10,10,args=(func,)))
    mean = quad(fmean,-2,2)[0]
    tmom = quad(f2mom,-2,2)[0]
    #print(mean)
    #print(tmom)

    #print(np.sqrt(tmom-mean**2))

    return (np.sqrt(tmom-mean**2))

def fsum(x,f):

    fsum=0
    for fi in f:
      fsum+=f[fi](x)

    return fsum

reason: import numpy,import scipy,from scipy