text:
import numpy as np
import glob
import sys
import getopt
import os
from astropy.io import fits
from scipy.interpolate import interp1d, splrep, splev, sproot
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
from matplotlib.widgets import Button, SpanSelector, CheckButtons
from matplotlib import gridspec
from mpl_toolkits.axes_grid.inset_locator import inset_axes
#from stsci.convolve import boxcar
from astropy.convolution import convolve, Box1DKernel

# Usage: run CrossCor.py imagename ('objid')
# Usage: run CrossCor.py imagename ('unsure')
# Usage: run CrossCor.py imagename ('inter')

# Requires the .big files in the same location as the cor_.fits files
'''
Inputs:
imagename - string -  set to be your big.fits file, including the path.
objid - optional string -  set to the 6-letter object id of a single object to open only that one
usure - optional string -  set to the string 'unsure' to open only those objects which the user flagged as unsure
inter - optional string -  set to the string 'inter' to run interactively; this is automatically on if looking at 1 objid
emiss - optional string -  set to the string 'emiss' to run through a list of emission-line only galaxies

Change to your dataset:
wave1,wave2 - (int,int) - the wavelength range over which to correlate the data, can be changed by the user in the GUI
zrange - (int,int) - the range of redshifts over which to search
emclip - boolean - set True to clip emission lines by default, or False to keep them in.
                   This setting can always be toggled within the GUI, but this value is what displays first.
specplot - boolean - set to 1 if you want a second plot showing common spectral lines overlaid on the galaxy at the currect redshift.
                     Can be kind of clunky, so default is off, but could be useful to enable for 'unsure' galaxies.
                     Can be toggled in the GUI. Requires galaxlylines.dat in the same folder as your image

Other:
outloc - string - where your files will be output, defaults to the same as input
temploc - string - where you store your templates. This defaults to you imagename location
outname - string - name of your output files. Defaults to 'cc_' + (your big.fits name). e.g. 'a' creates a.txt and verb_a.txt.
imloc - string - this will be set to the path to the big.fits file automatically
imname - string - this will be set to just the big.fits file automatically
tcorr - boolean - will be automatically set to 0 if the correction code has not been run.

NOTE: The outfile must correspond to the image that it was generated from.
      The file is generated once with all objects in the mask if it does not exist, and then is modified from there.
'''

wave1, wave2 = (3900, 8200)
zrange = (0.01, 0.4)
emclip = True
markbads = False  # Read in a separate file and auto-mark bads
haselist = False  # Read in list of objects to assume emission
specplot = 0

outloc = '/Users/galaxies-air/COSMOS/COSMOSData/corFitsFileOut/'
outname = 0
temploc = sys.path[0]+'/Templates/'
readunsure = 0
interactive = 1
objid = 0
tcorr = 1

fullCmdArguments = sys.argv
argumentList = fullCmdArguments[2:]
imagename = fullCmdArguments[1]
unixOptions = "iueo:b:l:t:"
gnuOptions = ["inter", "unsure", "emiss",
              "objid=", "bads=", "elist=", "temploc="]

try:
    arguments, values = getopt.getopt(argumentList, unixOptions, gnuOptions)
except getopt.error as err:
    # output error, and return with an error code
    print(str(err))
    sys.exit(2)

for currentArgument, currentValue in arguments:
    if currentArgument in ("-i", "--inter"):
        interactive = 1
    elif currentArgument in ("-u", "--unsure"):
        readunsure = 1
    elif currentArgument in ("-o", "--objid"):
        objid = currentValue
        interactive = 1
    elif currentArgument in ("-e", "--emiss"):
        emclip = False
    elif currentArgument in ("-b", "--bads"):
        badfile = currentValue
        markbads = True
    elif currentArgument in ("-l", "--elist"):
        emissfile = currentValue
        haselist = True
    elif currentArgument in ("-t", "--temploc"):
        temploc = currentValue

imname = imagename
while imname.find('/') != -1:
    imname = imname[imname.find('/')+1:]
imloc = imagename.replace(imname, '')
if not temploc:
    temploc = imloc
if not outloc:
    outloc = imloc
if not outname:
    outname = 'cc_' + imname.replace('.fits', '')

outfile = outloc + outname + '.txt'
outfilev = outloc + 'verb_' + outname + '.txt'

outfile = outloc + outname + '.txt'
# Trying to fix chanigng to a new computer, set to the mask you are doing
#outfilev = outfile.replace('cc_feb16_abig','verb_a6')
#outfilev = outfile.replace('cc_feb16_bbig','verb_b6')
#outfilev = outfile.replace('cc_feb17_bbig','verb_b7')
#outfilev = outfile.replace('cc_feb16_dbig','verb_d6')
#outfilev = outfile.replace('cc_feb16_ebig','verb_e6')
#outfilev = outfile.replace('cc_feb17_ebig','verb_e7')
#outfilev = outfile.replace('cc_feb16_fbig','verb_f6')
#outfilev = outfile.replace('cc_feb16_gbig','verb_g6')
#outfilev = outfile.replace('cc_feb16_hbig','verb_h6')
#outfilev = outfile.replace('cc_feb16_ibig','verb_i6')
outfilev = outfile.replace('cc_feb17_jbig', 'verb_j7')

# Check if telluric correction is done
if not glob.glob(imloc+'cor_??????_'+imname):
    tcorr = 0

# Check if the text files is set up. If not, create one
if not os.path.exists(outfile):
    print('No output text files exists. Creating ' + outfile + ' and ' + outfilev)
    if tcorr:
        imarr = glob.glob(imloc + 'cor_??????_' + imname)
    else:
        imarr = glob.glob(imloc + '??????_' + imname)
    f = open(outfile, 'w+')
    f2 = open(outfilev, 'w+')
    f.write('#OBJID  temp    z           dzhi        dzlo        ccmax     chi2    rchi2  eclip')
    f.write(
        '     S/N      Star Bad  Unsure    ImageName            Revisit  Note  Unusable\n')
    f2.write(
        '#OBJID z23 dzhi23 dzlo23 ccmax23 chi223 rchi223 z24 dzhi24 dzlo24 ccmax24 chi224 ')
    f2.write(
        'rchi224 z25 dzhi25 dzlo25 ccmax25 chi225 rchi225 z26 dzhi26 dzlo26 ccmax26 chi226 ')
    f2.write(
        'rchi226 z27 dzhi27 dzlo27 ccmax27 chi227 rchi227 temp eclip S/N Star Bad Unsure ')
    f2.write('ImageName Revisit Note Unusable\n')
    for i in range(len(imarr)):
        objname = imarr[i].replace(imloc, '')
        f.write(('%06d   %d    %.6f    %.6f   ' % (0, 0, 0, 0) +
                 '%.6f    %7.2f   %7.2f   %2.2f     ' % (0, 0, 0, 0) +
                 '%d      %2.2f      %d    %d     ' % (0, 0, 0, 0) +
                 '%d ' % (0)) + objname +
                ('      %d       %d       %d' % (0, 0, 0)) + '\n')
        f2.write(('%06d %.6f %.6f %.6f ' % (0, 0, 0, 0) +
                  '%7.3f %7.3f %2.2f %.6f ' % (0, 0, 0, 0) +
                  '%.6f %.6f %7.3f %7.3f ' % (0, 0, 0, 0) +
                  '%2.2f %.6f %.6f %.6f ' % (0, 0, 0, 0) +
                  '%7.3f %7.3f %2.2f %.6f ' % (0, 0, 0, 0) +
                  '%.6f %.6f %7.3f %7.3f ' % (0, 0, 0, 0) +
                  '%2.2f %.6f %.6f %.6f ' % (0, 0, 0, 0) +
                  '%7.3f %7.3f %2.2f %d ' % (0, 0, 0, 0) +
                  '%d %3.2f %d %d ' % (0, 0, 0, 0) +
                  '%d' % (0)) + ' ' + objname +
                 (' %d %d %d' % (0, 0, 0)) + '\n')
    f.close()
    f2.close()

skylines = [5577.0, 5890.0, 6300, 6364]
telluric = [(7584, 7650)]

def divz(X, Y):
    return X/np.where(Y, Y, Y+1)*np.not_equal(Y, 0)

def svdfit(b, y):
    decomp = np.linalg.svd(b, full_matrices=False)
    sol1 = np.transpose(decomp[2])
    sol2 = divz(1.0, decomp[1])
    sol3 = np.dot(np.transpose(decomp[0]), y)
    if np.sometrue(sol3):
        solr = (sol2*sol3)
        soll = np.dot(sol1, solr)
    else:
        soll = np.zeros(sol3.shape)
    return soll

def legendre(x, nl):
    x = np.asarray(x, dtype='float64')
    l = (nl+1,)+x.shape
    p = np.zeros(l, dtype='float64')
    p[0] = np.ones(x.shape, dtype='float64')
    if nl > 0:
        p[1] = x
    if nl > 1:
        for j in range(1, nl):
            p[j+1] = ((2*j+1)*x*p[j] - j*p[j-1])/(j+1)
    return p

def outGraph(xdata, ydata, objid, kind, lim=0):
    # kind is a short string describing the graph
    filepath = outloc
    fig, ax1 = plt.subplots(figsize=(10, 8))
    ax1.plot(xdata, ydata, color='red')
    ax1.set_title(('%06d' % objid) + ', ' + kind)
    if lim != 0:
        ax1.set_ylim(lim)
    fig.savefig(filepath + ('%06d' % objid) + '_' + kind)
    plt.close(fig)

def FitPoly(s, n, g, w, o, w1, w2):
    u = g*np.greater_equal(w, w1)*np.less_equal(w, w2)
    x = np.arange(len(s))
    xu, su, nu, gu, wu = np.compress(u, [x, s, n, g, w], 1)
    xmin, xmax = xu.min(), xu.max()
    xc = 0.5*(xmax+xmin)
    xs = 0.5*(xmax-xmin)
    xl = (xu-xc)/xs
    basis = legendre(xl, o)
    wt = divz(gu, nu)
    co = svdfit(np.transpose(basis)*wt[::, np.newaxis], su*wt)
    xl = (x-xc)/xs
    basis = legendre(xl, o)
    pl = np.dot(co, basis) * np.greater_equal(x, xmin) * np.less_equal(x,
                                                                       xmax) * np.greater_equal(w, w1) * np.less_equal(w, w2)
    return pl

def Smooth(y, good, p=50, h=25):
    m = np.zeros(y.shape, y.dtype)
    for j in range(len(y)):
        a, b = np.clip([j-h, j+h+1], 0, len(y)-1)
        u = np.compress(good[a:b], y[a:b])
        if len(u):
            if p == 50:
                m[j] = np.median(u)
            else:
                m[j] = np.sort(u)[len(u)*p/100]
    return m

class Spectrum():
    def __init__(self, file, clip=0, verb=0, good=0, changeclip=0):
        self.file = file
        self.bigfile = imloc + imname
        self.verb = verb
        self.clipemission = clip
        self.Read()
        self.bighdu = fits.open(self.bigfile)[0]
        self.image_data = self.bighdu.data[self.head['csecta']:(
            self.head['csectb']+1)]
        self.lobar = self.head['lowbound']-self.head['csecta']-1
        self.hibar = self.head['upbound']-self.head['csecta']+1
        self.bighead = self.bighdu.header
        if tcorr:
            self.spec = self.data[0]
            self.orig = self.data[1]
            self.noise = self.data[2]
            self.flat = self.data[3]
            self.sky = self.data[4]
            self.tc = self.data[5]
            self.mask = self.data[6]
            self.fit = self.data[7]
        else:
            self.spec = self.data[0]
            self.noise = self.data[1]
            self.flat = self.data[2]
            self.sky = self.data[3]
        self.objid = self.head['OBJID']
        self.findSN()
        self.GetBadPixels(good=good, changeclip=changeclip)

    def Read(self):
        self.hdu = fits.open(self.file)[0]
        self.head = self.hdu.header
        self.data = self.hdu.data.astype('float32')
        self.crval1 = self.head["crval1"]
        self.crpix1 = self.head["crpix1"]
        self.cdelt1 = self.head["cdelt1"]
        self.naxis1 = self.head["naxis1"]
        self.dcflag = self.head["dc-flag"]
        self.wavelength = (1.0+np.arange(self.naxis1) -
                           self.crpix1)*self.cdelt1 + self.crval1
        if self.dcflag:
            self.wavelength = np.power(10.0, self.wavelength)

    def findSN(self):
        sig = np.median(self.spec[450:650])
        noi = np.median(self.noise[450:650])
        if noi:
            self.signoise = sig/noi
        else:
            self.signoise = 0

    def GetBadPixels(self, good=0, changeclip=0):
        if hasattr(self, "noise"):
            #            bad = 1.*np.logical_or(np.less_equal(self.spec,1e-4),np.less_equal(self.noise,1e-4))
            bad = 1.*np.logical_or(np.less_equal(self.spec, -
                                                 3*self.noise), np.less_equal(self.noise, 1e-4))
            for skyline in skylines:
                bad = 1.*np.logical_or(bad, np.greater(self.wavelength, skyline -
                                                       3*self.cdelt1)*np.less(self.wavelength, skyline+3*self.cdelt1))
            for absorb in telluric:
                bad = 1.*np.logical_or(bad, np.greater(self.wavelength,
                                                       absorb[0])*np.less(self.wavelength, absorb[-1]))
            # bad = 1.*np.greater(boxcar(bad,(8,)),0)       #stsci boxcar
            # astropy boxcar
            bad = 1.*np.greater(convolve(bad, Box1DKernel(8)), 0)
        else:
            bad = np.zeros(self.spec.shape, dtype='float')
        sm = Smooth(self.spec, np.logical_not(bad))
        if self.verb:
            lim = (np.median(self.spec)*0.25, np.median(self.spec)*1.75)
            outGraph(self.wavelength, self.spec, self.objid, 'spec', lim=lim)
            outGraph(self.wavelength, sm, self.objid, 'sm50', lim=lim)
        sr = self.spec-sm
        ss = 1.49*Smooth(abs(sr), np.logical_not(bad), h=50)
        if self.verb:
            outGraph(self.wavelength, ss, self.objid, 'ss50')
        # KEEPS EMISSION LINES IN
        self.good2 = np.ones(sr.shape, dtype='float32')
        if changeclip:
            if self.clipemission:
                self.good = good*np.less(sr, 2.5*ss)
            else:
                self.good = good * \
                    np.ones(sr.shape, dtype='float32') * \
                    np.less(sr, 2.5*ss)+np.greater(sr, 2.5*ss)
            # NOTE THAT POLYNOMIAL FIT TO CONTINUUM MASKS OUT EMISSION/STRONG ABSORPTION
            self.continuum = 1.*np.less(abs(sr), 2.5*ss)
        else:
            self.good = np.logical_not(bad)
            if self.clipemission:
                # CLIPS OUT EMISSION LINES AND BAD POSITIVE SKY LINE RESIDUALS
                self.good = self.good*np.less(sr, 2.5*ss)
            else:
                # KEEPS EMISSION LINES IN
                self.good = self.good*np.ones(sr.shape, dtype='float32')
            # NOTE THAT POLYNOMIAL FIT TO CONTINUUM MASKS OUT EMISSION/STRONG ABSORPTION
            self.continuum = 1.*np.less(abs(sr), 2.5*ss)
            if self.verb:
                outGraph(self.wavelength, self.spec*self.good,
                         self.objid, 'sg50', lim=lim)
        self.sm = sm
        self.ss = ss

class Template(Spectrum):
    def __init__(self, n, clip=0, verb=0):
        self.file = temploc + "spDR2-%03d.fit" % (n)
        self.verb = verb
        self.clipemission = clip
        self.Read()
        self.spec = self.data[0]
        self.verb = 0
        self.GetBadPixels()
        wu, su = np.compress(self.good, [self.wavelength, self.spec], 1)
        self.tck = splrep(wu, su, s=0, k=3)
        self.interp = lambda w: splev(
            w, self.tck)*np.greater_equal(w, wu.min())*np.less_equal(w, wu.max())
        #self.interp = interp1d(wu,su,fill_value=0,bounds_error=0)

    def Redshift(self, z, w):
        return self.interp(w)/(1.0+z)

class Plot():
    def __init__(self, g):
        '''
        Variable initialization
        All of these should be set automatically, but for an explanation:
        image - image location + file, will be read as the spectrum
        objimage - just the image name
        verb - set to 1 to display a few more outputs
        eclip - whether to keep or clip emission lines
        wavechange - used when the user changes the wavelegnth from wave1 or wave2
        zbound - float,float - stores the lower and upper bound for the redshift range
        lowave and hiwave - the starting values of wave1 and wave2, but they get modified as the user changes them
        create - sets to 0 after creating the axes for the first time
        add - toggles between 1 (add) and 0 (subtract) for modifying the mask
        star,baddata,unsure,flag1,flag2,flag3 - 0 or 1, flags that the user can set
        failure - boolean - is set to 1 if the cross correlation fails, then the code outputs only zeros
        mkspec - boolean - toggles whether or not to display the spec plot
        '''
        self.image = g
        self.objimage = self.image
        while self.objimage.find('/') != -1:
            self.objimage = self.objimage[self.objimage.find('/')+1:]
        verb = 0
        eclip = emclip
        self.eclip = eclip
        self.wavechange = 0
        self.zbound = zrange
        self.lowave = wave1
        self.hiwave = wave2
        self.create = 1
        self.range = 100
        self.add = 0
        self.star = 0
        self.baddata = 0
        self.unsure = 0
        self.interactive = interactive
        self.failure = 0
        self.flag1 = 0
        self.flag2 = 0
        self.flag3 = 0
        self.mkspec = specplot

    def doCC(self, g, eclip, verb, newmask=0, wavechange=0, good=0, changeclip=0):
        '''
        Manages the cross correlation, sending the templates through one at a time and checking for failures.
        '''
        print('Computing Redshift for ' + self.image)
        if newmask:
            G = g
        else:
            G = Spectrum(g, clip=eclip, verb=verb,
                         good=good, changeclip=changeclip)
            self.G = G

        self.t23 = CCcalc(23, G, self.range, eclip, verb,
                          self.zbound, wavechange=wavechange)
        self.t24 = CCcalc(24, G, self.range, eclip, verb,
                          self.zbound, wavechange=wavechange)
        self.t25 = CCcalc(25, G, self.range, eclip, verb,
                          self.zbound, wavechange=wavechange)
        self.t26 = CCcalc(26, G, self.range, eclip, verb,
                          self.zbound, wavechange=wavechange)
        self.t27 = CCcalc(27, G, self.range, eclip, verb,
                          self.zbound, wavechange=wavechange)

        if self.t23.failure or self.t24.failure or self.t25.failure or self.t26.failure or self.t27.failure:
            self.failure = 1
            return

        maxarr = np.asarray((self.t23.ccmax, self.t24.ccmax,
                             self.t25.ccmax, self.t26.ccmax, self.t27.ccmax))
        chisqarr = np.asarray(
            (self.t23.rchi2, self.t24.rchi2, self.t25.rchi2, self.t26.rchi2, self.t27.rchi2))
        maxindex = np.argmax(maxarr)
        minindex = np.argmin(chisqarr)
        # If best chisq and cc_max at same point, it's easy
        if maxindex == minindex:
            self.tempid = maxindex + 23
        else:  # It just couldn't be easy, could it?
            bestzs = np.asarray(
                (self.t23.zsmax, self.t24.zsmax, self.t25.zsmax, self.t26.zsmax, self.t27.zsmax))
            z_o_a = 0.85 * (self.zbound[1] - self.zbound[0]) + self.zbound[0]
            zneighbors = [len(np.where(bestzs - bestzs[ii] <= 0.015)[0])
                          for ii in range(5)]
            # First: is one best guess at the upper-z bounds? It's probably worse
            if bestzs[maxindex] >= z_o_a and bestzs[minindex] < z_o_a:
                self.tempid = minindex + 23
            elif bestzs[maxindex] < z_o_a and bestzs[minindex] >= z_o_a:
                self.tempid = maxindex + 23
            # Second: is one a 3+ hit mode? Go with that.
            elif zneighbors[maxindex] >= 3:
                self.tempid = maxindex + 23
            elif zneighbors[minindex] >= 3:
                self.tempid = minindex + 23
            # Third: Is one a 2-hit mode and the other isn't? Try that
            elif zneighbors[maxindex] >= 2 and zneighbors[minindex] == 1:
                self.tempid = maxindex + 23
            elif zneighbors[minindex] >= 2 and zneighbors[maxindex] == 1:
                self.tempid = minindex + 23
            # Fourth: Okay, this is probably not too hot. Just go with the rchisq
            else:
                self.tempid = minindex + 23

        self.findTemp()

    def findTemp(self):
        if self.tempid == 23:
            self.temp = self.t23
        elif self.tempid == 24:
            self.temp = self.t24
        elif self.tempid == 25:
            self.temp = self.t25
        elif self.tempid == 26:
            self.temp = self.t26
        elif self.tempid == 27:
            self.temp = self.t27

    def mkspecplot(self):
        self.lfile = np.genfromtxt(
            imloc + 'galaxylines.dat', skip_header=11, dtype=None, names=('loc', 'type', 'name'))
        self.fig2 = plt.figure(figsize=(10, 6))
        self.a2x = self.fig2.add_axes([0, 0, 1, 1])
        G = self.G
        self.galplot = self.a2x.plot(G.wavelength, G.spec, color='black')
        xmin, xmax = self.a2x.get_xlim()
        self.a2x.set_xlim(xmin, xmax)
        self.a2x.set_ylim(np.median(G.spec)*-0.5, np.median(G.spec)*3.5)
        ymin, ymax = self.a2x.get_ylim()
        for i in range(len(self.lfile)):
            loc, ctype, name = self.lfile[i]
            shiftloc = loc*(1+self.temp.zsmax)
            if ctype == 2:
                color = 'indianred'
            else:
                color = 'cornflowerblue'
            if name == 'Break':
                color = 'mediumseagreen'
            elif name == 'Hbeta':
                color = 'cornflowerblue'
            self.a2x.plot((shiftloc, shiftloc), (-100000, 100000),
                          color=color, linestyle='--', alpha=0.75)

    def createPlot(self):
        '''
        Sets up the GUI for the first time
        '''
        if self.failure:
            G = self.G
            f3 = open(outfile, "r+")
            f4 = open(outfilev, 'r+')
            d = f3.readlines()
            e = f4.readlines()
            f3.seek(0)
            f4.seek(0)
            for i in d:
                if self.objimage in i:
                    f3.write(('%06d   %d    %.6f    %.6f   ' % (G.objid, 0, 0, 0) +
                              '%.6f    %7.2f   %7.2f   %2.2f     ' % (0, 0, 0, 0) +
                              '%d      %2.2f      %d    %d     ' % (0, 0, 0, 1) +
                              '%d ' % (1)) + self.objimage +
                             ('      %d      %d      %d' % (0, 0, 1)) + '\n')

                else:
                    f3.write(i)
            for i in e:
                if self.objimage in i:
                    f4.write(('%06d %.6f %.6f %.6f ' % (G.objid, 0, 0, 0) +
                              '%7.3f %7.3f %2.2f %.6f ' % (0, 0, 0, 0) +
                              '%.6f %.6f %7.3f %7.3f ' % (0, 0, 0, 0) +
                              '%2.2f %.6f %.6f %.6f ' % (0, 0, 0, 0) +
                              '%7.3f %7.3f %2.2f %.6f ' % (0, 0, 0, 0) +
                              '%.6f %.6f %7.3f %7.3f ' % (0, 0, 0, 0) +
                              '%2.2f %.6f %.6f %.6f ' % (0, 0, 0, 0) +
                              '%7.3f %7.3f %2.2f %d ' % (0, 0, 0, 0) +
                              '%d %3.2f %d %d ' % (0, 0, 0, 1) +
                              '%d' % (1)) + ' ' + self.objimage +
                             (' %d %d %d' % (0, 0, 1)) + '\n')
                else:
                    f4.write(i)
            f3.close()
            f4.close()
            return

        if self.create:
            G = self.G
            self.fig = plt.figure(figsize=(14, 8))
            self.ax1 = self.fig.add_axes([0.05, 0.7, 0.43, 0.25])
            self.ax2 = self.fig.add_axes([0.53, 0.7, 0.43, 0.20])
            self.ax2b = self.fig.add_axes([0.65, 0.40, 0.3, 0.2])
            self.ax2c = self.fig.add_axes([0.65, 0.1, 0.3, 0.2])
            self.ax3 = self.fig.add_axes([0.19, 0.1, 0.45, 0.50])
            self.ax5 = self.fig.add_axes([0.53, 0.9, 0.43, 0.08])
            self.ax1.set_title('Spectrum of ' + self.objimage +
                               ', S/N = ' + str(G.signoise))
            self.ax2.set_title('Galaxy Template Comparison')
            self.ax2b.set_title('H and K Line Zoom')
            self.ax2c.set_title('O3 and H-Beta Line Zoom')
            self.ax3.set_title('Cross Correlation                        ')
            self.refreshText()
            self.l1, = self.ax1.plot(
                G.wavelength, G.spec, color='cornflowerblue', label='Spectrum')
            self.l2, = self.ax1.plot(
                G.wavelength, G.good*G.spec, color='orange', label='Mask')
            self.l3, = self.ax1.plot(
                G.wavelength, G.sm, color='indianred', label='Smooth', visible=False)
            self.l4, = self.ax1.plot(
                G.wavelength, self.temp.Gp, color='darkgreen', label='Fit')
            self.l5, = self.ax1.plot(
                (self.lowave, self.lowave), (-100000, 100000), color='black')
            self.l6, = self.ax1.plot(
                (self.hiwave, self.hiwave), (-100000, 100000), color='black')
            self.ax1.set_ylim(np.median(G.spec)*-0.5, np.median(G.spec)*4)
            self.ax1.set_xlabel('$\lambda$')
            self.ax1.set_ylabel('Counts')
            self.ax1.legend(loc=1)
            self.k1, = self.ax2.plot(
                G.wavelength, self.temp.Gr, color='cornflowerblue')
            self.k2, = self.ax2.plot(
                G.wavelength, self.temp.Trs[self.temp.zp][self.temp.zmax], color='black', alpha=0.5)
            self.k1b, = self.ax2b.plot(
                G.wavelength, self.temp.Gr2, color='cornflowerblue')
            self.k2b, = self.ax2b.plot(
                G.wavelength, self.temp.Trs[self.temp.zp][self.temp.zmax], color='black', alpha=0.5)
            self.k1c, = self.ax2c.plot(
                G.wavelength, self.temp.Gr2, color='cornflowerblue')
            self.k2c, = self.ax2c.plot(
                G.wavelength, self.temp.Trs[self.temp.zp][self.temp.zmax], color='black', alpha=0.5)
            l1bloc = 3933.7*(1+self.temp.zsmax)
            l2bloc = 3968.5*(1+self.temp.zsmax)
            self.l1b, = self.ax2b.plot(
                (l1bloc, l1bloc), (-10, 10), color='mediumseagreen', linestyle='--')
            self.l2b, = self.ax2b.plot(
                (l2bloc, l2bloc), (-10, 10), color='mediumseagreen', linestyle='--')
            l1cloc = 4861.3*(1+self.temp.zsmax)
            l2cloc = 4959*(1+self.temp.zsmax)
            l3cloc = 5007*(1+self.temp.zsmax)
            self.l1c, = self.ax2c.plot(
                (l1cloc, l1cloc), (-10, 10), color='mediumseagreen', linestyle='--')
            self.l2c, = self.ax2c.plot(
                (l2cloc, l2cloc), (-10, 10), color='mediumseagreen', linestyle='--')
            self.l3c, = self.ax2c.plot(
                (l3cloc, l3cloc), (-10, 10), color='mediumseagreen', linestyle='--')
            self.ax2.set_xlim(self.temp.w1-100, self.temp.w2+100)
            self.ax2.set_ylim(-4, 5.5)
            self.ax2b.set_xlim(3800*(1+self.temp.zsmax),
                               4100*(1+self.temp.zsmax))
            self.ax2b.set_ylim(-4, 5.5)
            self.ax2c.set_xlim(4750*(1+self.temp.zsmax),
                               5050*(1+self.temp.zsmax))
            self.ax2c.set_ylim(-4, 5.5)
            self.ax2b.yaxis.tick_right()
            self.ax2c.yaxis.tick_right()
            self.ax2.set_xlabel('$\lambda$                                 ')
            self.ax2.set_ylabel('(g-f)/$\sigma$')
            self.ax2b.set_xlabel('$\lambda$')

            self.ax2c.set_xlabel('$\lambda$')
            self.m1, = self.ax3.plot(
                self.temp.zs[0], self.temp.cc[0], color='indianred', zorder=1)
            self.m1_23, = self.ax3.plot(self.t23.zs[0], self.t23.cc[0]*(max(self.temp.cc[0])/max(self.t23.cc[0])),
                                        color='#ff9933', alpha=0.35, zorder=0)
            self.m1_24, = self.ax3.plot(self.t24.zs[0], self.t24.cc[0]*(max(self.temp.cc[0])/max(self.t24.cc[0])),
                                        color='Yellow', alpha=0.35, zorder=0)
            self.m1_25, = self.ax3.plot(self.t25.zs[0], self.t25.cc[0]*(max(self.temp.cc[0])/max(self.t25.cc[0])),
                                        color='#33cc33', alpha=0.35, zorder=0)
            self.m1_26, = self.ax3.plot(self.t26.zs[0], self.t26.cc[0]*(max(self.temp.cc[0])/max(self.t26.cc[0])),
                                        color='#0066cc', alpha=0.35, zorder=0)
            self.m1_27, = self.ax3.plot(self.t27.zs[0], self.t27.cc[0]*(max(self.temp.cc[0])/max(self.t27.cc[0])),
                                        color='#9933ff', alpha=0.35, zorder=0)
            if self.temp.zp > 0:
                self.m2, = self.ax3.plot(
                    self.temp.zs[1], self.temp.cc[1], color='firebrick')
            else:
                self.m2, = self.ax3.plot(
                    (0, 0), (0, 0), color='firebrick', zorder=2)
            if self.temp.zp > 1:
                self.m3, = self.ax3.plot(
                    self.temp.zs[2], self.temp.cc[2], color='darkred')
            else:
                self.m3, = self.ax3.plot(
                    (0, 0), (0, 0), color='darkred', zorder=2)
            if self.temp.zp > 2:
                self.m4, = self.ax3.plot(
                    self.temp.zs[3], self.temp.cc[3], color='darkred')
            else:
                self.m4, = self.ax3.plot(
                    (0, 0), (0, 0), color='darkred', zorder=2)
            self.g1, = self.ax3.plot((self.temp.zs[self.temp.zp][self.temp.zmax],
                                      self.temp.zs[self.temp.zp][self.temp.zmax]),
                                     (-10000, 10000), color='mediumseagreen', zorder=1)
            self.ax3.set_xlim(self.zbound[0], self.zbound[1])
            self.ax3.set_ylim(
                min(self.temp.cc[0])*0.95, max(self.temp.cc[0]*1.05))
            self.ax3.set_xlabel('z')
            self.ax3.set_ylabel('Correlation Coeffecient')
            self.ax5.imshow(G.image_data, cmap='gray',
                            clim=(-240.0, 240.0), aspect='auto')
            self.ax5.plot((0, 12000), (G.lobar, G.lobar),
                          color='mediumseagreen')
            self.ax5.plot((0, 12000), (G.hibar, G.hibar),
                          color='mediumseagreen')
            self.ax5.get_xaxis().set_visible(False)
            self.ax5.get_yaxis().set_visible(False)
            ax2min, ax2max = self.ax2.get_xlim()
            self.ax5.set_xlim((ax2min-4900)/2, (ax2max-4900)/2)
            rax = plt.axes([0.01, 0.34, 0.08, 0.12])
            check = CheckButtons(
                rax, ('Spectrum', 'Mask', 'Smooth', 'Fit'), (True, True, False, True))
            self.create = 0

# ==================================================================
#Buttons and sliders
# ==================================================================
        def checkButton(label):
            if label == 'Spectrum':
                self.l1.set_visible(not self.l1.get_visible())
            elif label == 'Mask':
                self.l2.set_visible(not self.l2.get_visible())
            elif label == 'Smooth':
                self.l3.set_visible(not self.l3.get_visible())
            elif label == 'Fit':
                self.l4.set_visible(not self.l4.get_visible())
            self.ax1.legend(loc=1)
            plt.draw()
        check.on_clicked(checkButton)

        def ccslider(xmin, xmax):
            indmin, indmax = np.searchsorted(self.temp.zs[0], (xmin, xmax))
            indmax = min(len(self.temp.zs[0]) - 1, indmax)
            self.newcc = np.mean((indmin, indmax))
            self.temp.rerunCC(self.newcc, self.range, self.zbound)
            self.updatePlot()
        span3 = SpanSelector(self.ax3, ccslider, 'horizontal',
                             useblit=True, rectprops=dict(alpha=0.5, facecolor='red'))

        def waveChangeHiF(event):
            self.hiwave += 50
            updateWave()
        wavechangehif = plt.axes([0.298, 0.643, 0.015, 0.02])
        bwavechangehif = Button(wavechangehif, '>')
        bwavechangehif.on_clicked(waveChangeHiF)

        def waveChangeHiB(event):
            self.hiwave -= 50
            updateWave()
        wavechangehib = plt.axes([0.282, 0.643, 0.015, 0.02])
        bwavechangehib = Button(wavechangehib, '<')
        bwavechangehib.on_clicked(waveChangeHiB)

        def waveChangeLoF(event):
            self.lowave += 50
            updateWave()
        wavechangelof = plt.axes([0.232, 0.643, 0.015, 0.02])
        bwavechangelof = Button(wavechangelof, '>')
        bwavechangelof.on_clicked(waveChangeLoF)

        def waveChangeLoB(event):
            self.lowave -= 50
            updateWave()
        wavechangelob = plt.axes([0.216, 0.643, 0.015, 0.02])
        bwavechangelob = Button(wavechangelob, '<')
        bwavechangelob.on_clicked(waveChangeLoB)

        def waveChangeReset(event):
            self.lowave = wave1
            self.hiwave = wave2
            updateWave()
        wavechangereset = plt.axes([0.176, 0.643, 0.036, 0.02])
        bwavechangereset = Button(wavechangereset, 'Reset')
        bwavechangereset.on_clicked(waveChangeReset)

        def updateWave():
            self.wavechange = (self.lowave, self.hiwave)
            self.l5.set_xdata((self.lowave, self.lowave))
            self.l6.set_xdata((self.hiwave, self.hiwave))
            plt.draw()

        def adjRange(event):
            self.range -= 40
            if self.range == -20:
                self.range = 100
            badjrange.label.set_text('Range: ' + ('%3d' % self.range))
            plt.draw()
        adjrange = plt.axes([0.01, 0.11, 0.12, 0.035])
        badjrange = Button(adjrange, 'Range: ' + ('%3d' % self.range))
        badjrange.on_clicked(adjRange)

        def switchTempF(event):
            self.tempid += 1
            if self.tempid == 28:
                self.tempid = 23
            bswitchtemp.label.set_text('Template = ' + ('%2d' % self.tempid))
            self.findTemp()
            self.updatePlot()

        def switchTempB(event):
            self.tempid -= 1
            if self.tempid == 22:
                self.tempid = 27
            bswitchtemp.label.set_text('Template = ' + ('%2d' % self.tempid))
            self.findTemp()
            self.updatePlot()
        switchtemp = plt.axes([0.025, 0.15, 0.09, 0.055])
        bswitchtemp = Button(switchtemp, 'Template = ' + ('%2d' % self.tempid))
        switchtempF = plt.axes([0.118, 0.15, 0.012, 0.055])
        bswitchtempF = Button(switchtempF, '>')
        bswitchtempF.on_clicked(switchTempF)
        switchtempB = plt.axes([0.010, 0.15, 0.012, 0.055])
        bswitchtempB = Button(switchtempB, '<')
        bswitchtempB.on_clicked(switchTempB)

        def toggleClip(event):
            self.eclip = not self.eclip
            btoggleclip.label.set_text('Clip = ' + str(self.eclip))
            self.worktext = self.fig.text(0.5, 0.01, 'Working...', fontsize=24)
            plt.pause(.1)
            self.doCC(self.image, self.eclip, 0,
                      wavechange=self.wavechange, good=self.G.good, changeclip=1)
            # self.doCC(self.image,self.eclip,0)
            plt.gcf().texts.remove(self.worktext)
            bswitchtemp.label.set_text('Template = ' + ('%2d' % self.tempid))
            G = self.G
            self.l2.set_ydata(G.good*G.spec)
            self.l3.set_ydata(G.sm)
            self.l4.set_ydata(self.temp.Gp)
            self.k1.set_ydata(self.temp.Gr)
            self.k1b.set_ydata(self.temp.Gr2)
            self.k1c.set_ydata(self.temp.Gr2)
            self.updatePlot()
        toggleclip = plt.axes([0.010, 0.21, 0.12, 0.035])
        btoggleclip = Button(toggleclip, 'Clip = ' + str(self.eclip))
        btoggleclip.on_clicked(toggleClip)

        def maskMode(event):
            self.add = not self.add
            if self.add:
                bmode.label.set_text('Mode: Add')
            else:
                bmode.label.set_text('Mode: Sub')
            plt.draw()
        mode = plt.axes([0.010, 0.31, 0.06, 0.035])
        bmode = Button(mode, 'Mode: Sub')
        bmode.on_clicked(maskMode)

        def maskslider(xmin, xmax):
            indmin, indmax = np.searchsorted(G.wavelength, (xmin, xmax))
            indmax = min(len(G.wavelength) - 1, indmax)
            if indmin == indmax:
                hidiff = abs(self.hiwave - G.wavelength[indmin])
                lodiff = abs(self.lowave - G.wavelength[indmin])
                if hidiff > lodiff:
                    self.lowave = G.wavelength[indmin]
                else:
                    self.hiwave = G.wavelength[indmin]
                updateWave()
            else:
                if self.add == 1:
                    self.G.good = np.concatenate(
                        (self.G.good[:indmin], np.ones(indmax-indmin), self.G.good[indmax:]))
                    self.G.continuum = np.concatenate(
                        (self.G.continuum[:indmin], np.ones(indmax-indmin), self.G.continuum[indmax:]))
                if self.add == 0:
                    self.G.good = np.concatenate(
                        (self.G.good[:indmin], np.zeros(indmax-indmin), self.G.good[indmax:]))
                    self.G.continuum = np.concatenate(
                        (self.G.continuum[:indmin], np.zeros(indmax-indmin), self.G.continuum[indmax:]))
                self.l2.set_ydata(self.G.good*self.G.spec)
                plt.draw()
        span = SpanSelector(self.ax1, maskslider, 'horizontal',
                            useblit=True, rectprops=dict(alpha=0.5, facecolor='red'))

        def buttonGo(event):
            self.worktext = self.fig.text(0.5, 0.01, 'Working...', fontsize=24)
            plt.pause(.1)
            self.doCC(self.G, self.eclip, 0, newmask=1,
                      wavechange=self.wavechange)
            plt.gcf().texts.remove(self.worktext)
            bswitchtemp.label.set_text('Template = ' + ('%2d' % self.tempid))
            self.k1.set_ydata(self.temp.Gr)
            self.k1b.set_ydata(self.temp.Gr2)
            self.k1c.set_ydata(self.temp.Gr2)
            self.l4.set_ydata(self.temp.Gp)
            self.updatePlot()
        go = plt.axes([0.070, 0.31, 0.02, 0.035])
        bgo = Button(go, 'Go')
        bgo.on_clicked(buttonGo)

        def Bsave(event):
            self.savetext = self.fig.text(0.5, 0.01, 'Saving...', fontsize=24)
            plt.pause(.1)
            plt.gcf().texts.remove(self.savetext)
            G = self.G
            f3 = open(outfile, "r+")
            f4 = open(outfilev, 'r+')
            d = f3.readlines()
            e = f4.readlines()
            f3.seek(0)
            f4.seek(0)
            for i in d:
                if self.objimage in i:
                    f3.write(('%06d   %d    %.6f    %.6f   ' % (G.objid, self.tempid, self.temp.zsmax, self.temp.dzhi) +
                              '%.6f    %7.2f   %7.2f   %2.2f     ' % (self.temp.dzlo, self.temp.ccmax, self.temp.chi2, self.temp.rchi2) +
                              '%d      %2.2f      %d    %d     ' % (self.eclip, G.signoise, self.star, self.baddata) +
                              '%d ' % (self.unsure)) + self.objimage +
                             ('      %d      %d      %d' % (self.flag1, self.flag2, self.flag3)) + '\n')
                else:
                    f3.write(i)
            for i in e:
                if self.objimage in i:
                    f4.write(('%06d %.6f %.6f %.6f ' % (G.objid, self.t23.zsmax, self.t23.dzhi, self.t23.dzlo) +
                              '%7.3f %7.3f %2.2f %.6f ' % (self.t23.ccmax, self.t23.chi2, self.t23.rchi2, self.t24.zsmax) +
                              '%.6f %.6f %7.3f %7.3f ' % (self.t24.dzhi, self.t24.dzlo, self.t24.ccmax, self.t24.chi2) +
                              '%2.2f %.6f %.6f %.6f ' % (self.t24.rchi2, self.t25.zsmax, self.t25.dzhi, self.t25.dzlo) +
                              '%7.3f %7.3f %2.2f %.6f ' % (self.t25.ccmax, self.t25.chi2, self.t25.rchi2, self.t26.zsmax) +
                              '%.6f %.6f %7.3f %7.3f ' % (self.t26.dzhi, self.t26.dzlo, self.t26.ccmax, self.t26.chi2) +
                              '%2.2f %.6f %.6f %.6f ' % (self.t26.rchi2, self.t27.zsmax, self.t27.dzhi, self.t27.dzlo) +
                              '%7.3f %7.3f %2.2f %d ' % (self.t27.ccmax, self.t27.chi2, self.t27.rchi2, self.tempid) +
                              '%d %3.2f %d %d ' % (self.eclip, G.signoise, self.star, self.baddata) +
                              '%d' % (self.unsure)) + ' ' + self.objimage +
                             (' %d %d %d' % (self.flag1, self.flag2, self.flag3)) + '\n')

                else:
                    f4.write(i)
            f3.close()
            f4.close()
            filelocation2 = outloc + 'cc_' + \
                self.objimage.replace('.fits', '.png')
            filelocation3 = outloc + 'cc_' + \
                self.objimage.replace('.fits', '2.png')
            self.fig.savefig(filelocation2)
            if self.mkspec:
                self.fig2.savefig(filelocation3)
            plt.clf()
            plt.close(self.fig)
            if self.mkspec:
                plt.close(self.fig2)

        def checkStarFlag():
            if self.star == 1:
                bsflag.label.set_text('Unmark Star')
                bsflag.color = 'IndianRed'
                bsflag.hovercolor = 'Red'
            else:
                bsflag.label.set_text('Mark Star')
                bsflag.color = '0.85'
                bsflag.hovercolor = '0.95'
            plt.draw()

        def starFlag(event):
            if self.star == 1:
                self.star = 0
            else:
                self.star = 1
            checkStarFlag()
        sflag = plt.axes([0.15, 0.0, 0.08, 0.05])
        if G.spec[1000:1500].mean() > 10000:
            self.star = 1
            bsflag = Button(sflag, 'Mark Star',
                            color='IndianRed', hovercolor='Red')
        else:
            self.star = 0
            bsflag = Button(sflag, 'Mark Star',
                            color='0.85', hovercolor='0.95')
        bsflag.on_clicked(starFlag)
        checkStarFlag()

        def dataFlag(event):
            if self.baddata == 1:
                self.baddata = 0
                bdflag.label.set_text('Mark Bad')
                bdflag.color = '0.85'
                bdflag.hovercolor = '0.95'
            else:
                self.baddata = 1
                bdflag.label.set_text('Unmark Bad')
                bdflag.color = 'IndianRed'
                bdflag.hovercolor = 'Red'
            plt.draw()
        dflag = plt.axes([0.25, 0.0, 0.08, 0.05])
        if self.baddata == 0:
            bdflag = Button(dflag, 'Mark Bad', color='0.85', hovercolor='0.95')
        else:
            bdflag = Button(dflag, 'Unmark Bad',
                            color='IndianRed', hovercolor='Red')
        bdflag.on_clicked(dataFlag)

        def unsureFlag(event):
            if self.unsure == 1:
                self.unsure = 0
                buflag.label.set_text('Mark Unsure')
            else:
                self.unsure = 1
                buflag.label.set_text('Unmark Unsure')
            plt.draw()
        uflag = plt.axes([0.35, 0.0, 0.08, 0.05])
        buflag = Button(uflag, 'Mark Unsure')
        buflag.on_clicked(unsureFlag)

        def eFlag1(event):
            if self.flag1 == 1:
                self.flag1 = 0
                bflag1.label.set_text('Mark Revisit')
            else:
                self.flag1 = 1
                bflag1.label.set_text('Unmark Revisit')
            plt.draw()
        flag1 = plt.axes([0.65, 0.0, 0.08, 0.05])
        bflag1 = Button(flag1, 'Mark Revisit')
        bflag1.on_clicked(eFlag1)

        def eFlag2(event):
            if self.flag2 == 1:
                self.flag2 = 0
                bflag2.label.set_text('Mark Note')
            else:
                self.flag2 = 1
                bflag2.label.set_text('Unmark Range')
            plt.draw()
        flag2 = plt.axes([0.75, 0.0, 0.08, 0.05])
        bflag2 = Button(flag2, 'Mark Range')
        bflag2.on_clicked(eFlag2)

        def eFlag3(event):
            if self.flag3 == 1:
                self.flag3 = 0
                bflag3.label.set_text('Mark Unusable')
            else:
                self.flag3 = 1
                bflag3.label.set_text('Unmark Unusable')
            plt.draw()
        flag3 = plt.axes([0.85, 0.0, 0.08, 0.05])
        bflag3 = Button(flag3, 'Mark Unusable')
        bflag3.on_clicked(eFlag3)

        def changeZ(event):
            self.inputtext = self.fig.text(
                0.45, 0.01, 'Input Lower Bound (float)', fontsize=14)
            plt.pause(.1)
            loz = input('Input lower bound (float): ')
            plt.gcf().texts.remove(self.inputtext)
            self.inputtext = self.fig.text(
                0.45, 0.01, 'Input Upper Bound (float)', fontsize=14)
            plt.pause(.1)
            hiz = input('Input upper bound (float): ')
            plt.gcf().texts.remove(self.inputtext)
            self.zbound = loz, hiz
            self.worktext = self.fig.text(0.5, 0.01, 'Working...', fontsize=24)
            plt.pause(.1)
            self.doCC(self.G, self.eclip, 0, newmask=1,
                      wavechange=self.wavechange)
            plt.gcf().texts.remove(self.worktext)
            bswitchtemp.label.set_text('Template = ' + ('%2d' % self.tempid))
            self.k1.set_ydata(self.temp.Gr)
            self.k1b.set_ydata(self.temp.Gr2)
            self.k1c.set_ydata(self.temp.Gr2)
            self.l4.set_ydata(self.temp.Gp)
            self.updatePlot()
        changez = plt.axes([0.0, 0.06, 0.1, 0.03])
        bchangez = Button(changez, 'Change z Range')
        bchangez.on_clicked(changeZ)

        def checkSpecFlag():
            if self.star == 1:
                bsflag.label.set_text('Unmark Star')
                bsflag.color = 'Red'
            else:
                bsflag.label.set_text('Mark Star')
                bsflag.color = 0.85
            plt.draw()

        def mkSpec(event):
            if self.mkspec == 1:
                self.mkspec = 0
                plt.close(self.fig2)
                bmkspec.label.set_text('Show Spec Plot')
            else:
                self.mkspec = 1
                self.mkspecplot()
                plt.show(block=False)
                plt.figure(1)
                bmkspec.label.set_text('Hide Spec Plot')
            plt.draw()
        mkspec = plt.axes([0.01, 0.63, 0.08, 0.03])
        if specplot:
            bmkspec = Button(mkspec, 'Hide Spec Plot')
        else:
            bmkspec = Button(mkspec, 'Show Spec Plot')
        bmkspec.on_clicked(mkSpec)

        save = plt.axes([0.0, 0.0, 0.1, 0.05])
        if self.interactive:
            bsave = Button(save, 'Save')
            bsave.on_clicked(Bsave)
            if self.mkspec:
                self.mkspecplot()
            plt.show()
        else:
            Bsave(self)

# ==============================================================

    def updatePlot(self):
        '''
        Updates the plot as changes are made (e.g. toggling the emission clip or changing the wavelength range)
        '''
        self.k2.set_ydata(self.temp.Trs[self.temp.zp][self.temp.zmax])
        self.k2b.set_ydata(self.temp.Trs[self.temp.zp][self.temp.zmax])
        self.k2c.set_ydata(self.temp.Trs[self.temp.zp][self.temp.zmax])
        self.m1.set_ydata(self.temp.cc[0])
        self.m1.set_xdata(self.temp.zs[0])
        self.m1_23.set_ydata(
            self.t23.cc[0]*(max(self.temp.cc[0])/max(self.t23.cc[0])))
        self.m1_23.set_xdata(self.t23.zs[0])
        self.m1_24.set_ydata(
            self.t24.cc[0]*(max(self.temp.cc[0])/max(self.t24.cc[0])))
        self.m1_24.set_xdata(self.t24.zs[0])
        self.m1_25.set_ydata(
            self.t25.cc[0]*(max(self.temp.cc[0])/max(self.t25.cc[0])))
        self.m1_25.set_xdata(self.t25.zs[0])
        self.m1_26.set_ydata(
            self.t26.cc[0]*(max(self.temp.cc[0])/max(self.t26.cc[0])))
        self.m1_26.set_xdata(self.t26.zs[0])
        self.m1_27.set_ydata(
            self.t27.cc[0]*(max(self.temp.cc[0])/max(self.t27.cc[0])))
        self.m1_27.set_xdata(self.t27.zs[0])
        if self.temp.zp > 0:
            self.m2.set_ydata(self.temp.cc[1])
            self.m2.set_xdata(self.temp.zs[1])
        else:
            self.m2.set_ydata((0, 0))
            self.m2.set_xdata((0, 0))
        if self.temp.zp > 1:
            self.m3.set_ydata(self.temp.cc[2])
            self.m3.set_xdata(self.temp.zs[2])
        else:
            self.m3.set_ydata((0, 0))
            self.m3.set_xdata((0, 0))
        if self.temp.zp > 2:
            self.m4.set_ydata(self.temp.cc[3])
            self.m4.set_xdata(self.temp.zs[3])
        else:
            self.m4.set_ydata((0, 0))
            self.m4.set_xdata((0, 0))
        self.g1.set_xdata(
            (self.temp.zs[self.temp.zp][self.temp.zmax], self.temp.zs[self.temp.zp][self.temp.zmax]))
        self.ax3.set_xlim(self.zbound[0], self.zbound[1])
        self.ax3.set_ylim(min(self.temp.cc[0])*0.95, max(self.temp.cc[0]*1.05))
        self.ax2.set_xlim(self.temp.w1-100, self.temp.w2+100)
        l1bloc = 3933.7*(1+self.temp.zsmax)
        l2bloc = 3968.5*(1+self.temp.zsmax)
        self.l1b.set_xdata((l1bloc, l1bloc))
        self.l2b.set_xdata((l2bloc, l2bloc))
        l1cloc = 4861.3*(1+self.temp.zsmax)
        l2cloc = 4959*(1+self.temp.zsmax)
        l3cloc = 5007*(1+self.temp.zsmax)
        self.l1c.set_xdata((l1cloc, l1cloc))
        self.l2c.set_xdata((l2cloc, l2cloc))
        self.l3c.set_xdata((l3cloc, l3cloc))
        self.ax2b.set_xlim(3800*(1+self.temp.zsmax), 4100*(1+self.temp.zsmax))
        self.ax2c.set_xlim(4750*(1+self.temp.zsmax), 5050*(1+self.temp.zsmax))
        ax2min, ax2max = self.ax2.get_xlim()
        self.ax5.set_xlim((ax2min-4900)/2, (ax2max-4900)/2)
        self.refreshText(remove=1)
        if self.mkspec:
            plt.close(self.fig2)
            self.mkspecplot()
            plt.show(block=False)
            plt.figure(1)
        plt.draw()

    def refreshText(self, remove=0):
        if remove:
            plt.gcf().texts.remove(self.ztext)
            plt.gcf().texts.remove(self.zhitext)
            plt.gcf().texts.remove(self.zlotext)
            plt.gcf().texts.remove(self.temptext)
            plt.gcf().texts.remove(self.linetext)
            plt.gcf().texts.remove(self.t23text)
            plt.gcf().texts.remove(self.t24text)
            plt.gcf().texts.remove(self.t25text)
            plt.gcf().texts.remove(self.t26text)
            plt.gcf().texts.remove(self.t27text)
        self.ztext = self.fig.text(
            0.4285, 0.610, "z = %.6f" % (self.temp.zsmax), fontsize=12)
        self.zhitext = self.fig.text(
            0.5095, 0.624, "+%.6f" % (self.temp.dzhi), fontsize=8)
        self.zlotext = self.fig.text(
            0.513, 0.604, "-%.6f" % (self.temp.dzlo), fontsize=8)
        self.temptext = self.fig.text(
            0.005, 0.58, '$T$     $ccmax$         $z$                 $\chi^2$', fontsize=9)
        self.linetext = self.fig.text(
            0.005, 0.58, '__   _______    _________     ______', fontsize=9)
        self.t23text = self.fig.text(0.005, 0.56, "23: %7.2f   %.6f   %7.2f" % (
            self.t23.ccmax, self.t23.zsmax, self.t23.chi2), fontsize=9)
        self.t24text = self.fig.text(0.005, 0.545, "24: %7.2f   %.6f   %7.2f" % (
            self.t24.ccmax, self.t24.zsmax, self.t24.chi2), fontsize=9)
        self.t25text = self.fig.text(0.005, 0.53, "25: %7.2f   %.6f   %7.2f" % (
            self.t25.ccmax, self.t25.zsmax, self.t25.chi2), fontsize=9)
        self.t26text = self.fig.text(0.005, 0.515, "26: %7.2f   %.6f   %7.2f" % (
            self.t26.ccmax, self.t26.zsmax, self.t26.chi2), fontsize=9)
        self.t27text = self.fig.text(0.005, 0.50, "27: %7.2f   %.6f   %7.2f" % (
            self.t27.ccmax, self.t27.zsmax, self.t27.chi2), fontsize=9)

class CCcalc:
    '''
    Performs the cross correlation
    G - galaxy spectrum after going the class Spectrum()
    range - int - number of points per iteration.
            Default 100, user can set to 60 or 20 in the GUI if they want to zoom in on a small peak that is beside a large one
    t - int - template number, either 23,24,25,26,27
    failure - boolean - set to 1 if the correlation fails
    z1, z2 - float,float - bounds over which to fit redshift
    dz0 - float - starting stepsize for the first iteration
    nz - int - number of steps
    zp - int - iteration number (first iteration is 0)
    keepgoing - boolean - true if another iteration should be done,
                automatically set to false if enough iterations have been done to accurately compute errors
    w1,w2 - int,int - wavelength range over which to correlate. Data outside of this range is ignored.
    wavechange - int,int - if the user has modified the wavelength in the GUI, wavechange contains the new lower and upper bound
    o - int - order of the polynomial to fit to the galaxy and templates. Scales with the length of data over which we are fitting
    cc - array of correlation coefficients
    zs - array of redshifts corresponding to those coefficients
    Trs - array of normalized templates shifted to the redshifts in zs
    '''

    def __init__(self, t, G, inrange, eclip, verb, zbound, wavechange=0):
        self.G = G
        self.range = inrange
        self.t = t
        self.failure = 0

        T = Template(t, clip=eclip, verb=verb)
        self.T = T

        z1, z2 = zbound
        dz0 = 4*G.cdelt1/6500.0
        nz = int((z2-z1)/dz0+1)
        zp = 0
        keepgoing = True
        w1, w2 = wave1, wave2
        if wavechange:
            self.w1 = wavechange[0]
            self.w2 = wavechange[1]
        else:
            self.w1 = w1
            self.w2 = w2
        o = int((w2-w1)/250)
        self.o = o

        self.cc = []
        self.zs = []
        self.Trs = []

        dz = (z2-z1)/(nz-1)

        self.computeCC(z1, z2, dz, zp, keepgoing, self.w1,
                       self.w2, o, G, T, self.t, verb, intro=1)

    def rerunCC(self, newcc, inrange, zbound):
        '''
        This is done if the user clicks on a different peak -
        instead of rerunning the whole template, we can skip the first run through and start at the second, near where the user clicked.
        See above docstring for explanation of variables
        '''
        self.range = inrange
        self.cc = [self.cc[0]]
        self.zs = [self.zs[0]]
        self.Trs = [self.Trs[0]]
        zmax = self.zs[0][int(newcc)]
        z1, z2 = zbound
        nz = 100
        dz = (z2-z1)/(nz-1)
        if self.range == 100:
            z1, z2 = zmax-5*dz, zmax+4.9*dz
        elif self.range == 60:
            z1, z2 = zmax-3*dz, zmax+2.9*dz
        elif self.range == 20:
            z1, z2 = zmax-1*dz, zmax+0.9*dz
        dz = dz/10.0
        zp = 1
        keepgoing = True
        self.zp = zp
        self.computeCC(z1, z2, dz, zp, keepgoing, self.w1,
                       self.w2, self.o, self.G, self.T, self.t, 0)

    def computeCC(self, z1, z2, dz, zp, keepgoing, w1, w2, o, G, T, t, verb, intro=0):
        '''
        Main iterative calculations are done here
        '''

        dGw = G.wavelength[1:]-G.wavelength[:-1]
        dGw = np.compress(np.greater(
            G.wavelength[1:], w1)*np.greater(G.wavelength[:-1], w2), dGw)
        Gw = np.compress(np.greater(G.wavelength[1:], w1)*np.greater(
            G.wavelength[:-1], w2), (G.wavelength[1:]+G.wavelength[:-1])/2.)
        dzfloor = dGw.mean()/Gw.mean()

        while keepgoing:
            if intro:
                # print z1,z2,dzfloor,dz
                zs = np.arange(z1, z2+dz/2, dz, dtype='float64')
            else:
                # print zsmax,z1,z2,dzfloor,dz
                zs = np.sort(np.concatenate([np.arange(
                    zsmax, z2+dz/2, dz, dtype='float64'), np.arange(zsmax-dz, z1-dz/2, -dz, dtype='float64')]))
            self.zs.append(zs)

            if verb:
                print("Using order=", o)

            try:
                Gp = FitPoly(G.spec, G.noise, G.continuum,
                             G.wavelength, o, w1, w2)
            except(ValueError):
                print('Could not fit polynomial')
                keepgoing = 0
                self.failure = 1
            if not self.failure:
                self.Gp = Gp
                Gm = 1.*G.good * \
                    np.greater_equal(G.wavelength, w1) * \
                    np.less_equal(G.wavelength, w2)
                Gw = divz(Gm, G.noise)
                Gn = divz(np.sum(G.spec*Gw), np.sum(Gw))
                if verb:
                    lim = (np.median(G.spec)*0.25, np.median(G.spec)*1.75)
                    outGraph(G.wavelength, Gp, G.objid, 'pl', lim=lim)

                Gp2 = FitPoly(G.spec, G.noise, G.continuum, G.wavelength, int(
                    (1e4-min([w1, 4900]))/250.), min([w1, 4900]), 10000)
                Gm2 = 1.*G.good2
                Gw2 = divz(Gm2, G.noise)
                Gr2 = (G.spec-Gp2)*Gw2
                self.Gr2 = Gr2

                Tzs = np.asarray(
                    [T.Redshift(z, G.wavelength/(1+z)) for z in zs])
                Tzm = np.add.reduce(Tzs*Gw[np.newaxis, ::], 1)/np.sum(Gw)
                Tzs = (Tzs/Tzm[::, np.newaxis]) * Gn

                Tps = np.asarray(
                    [FitPoly(Tz, G.noise, G.continuum, G.wavelength, o, w1, w2) for Tz in Tzs])

                Gr = (G.spec-Gp)*Gw
                Trs = (Tzs-Tps)*Gw
                self.Gr = Gr
                self.Trs.append(Trs)

                cc = np.add.reduce(Trs*Gr[np.newaxis, ::], 1)
                self.cc.append(cc)
                if verb:
                    outGraph(zs, cc, G.objid, "c_temp%d_%d" % (t, zp))

                zmax = np.argmax(cc)
                self.zmax = zmax
                # print "(%d,%d,%.5f)" % (t,zp,zs[zmax])
                try:
                    dz1 = cc[zmax]-cc[zmax-1]
                    dz2 = cc[zmax]-cc[zmax+1]
                except(IndexError):
                    dz1 = 1
                    dz2 = 1
                self.dz1 = dz1
                self.dz2 = dz2
                if dz1 < 0.5 and dz2 < 0.5:
                    keepgoing = False

                if intro:
                    #                    dz = 0.3/99.0
                    intro = 0

                if keepgoing:
                    cctck = splrep(zs, cc, s=0, k=3)
                    zsmax = fsolve(lambda z: splev(
                        z, cctck, 1), (zs[zmax],))[0]

                    dz = dz/10.0
                    dz = max([dz, dzfloor])
                    z1, z2 = zs[zmax]-10*dz, zs[zmax]+9.9*dz
#                    if self.range == 100:
#                        z1,z2 = zs[zmax]-5*dz,zs[zmax]+4.9*dz
#                    elif self.range == 60:
#                        z1,z2 = zs[zmax]-3*dz,zs[zmax]+2.9*dz
#                    elif self.range == 20:
#                        z1,z2 = zs[zmax]-1*dz,zs[zmax]+0.9*dz
#                    dz = dz/2.5
#                    dz = max([dz,dzfloor])
                    if zp < 10:
                        zp += 1
                    else:
                        keepgoing = False

#                    zp += 1

        if not self.failure:
            '''
            Our results of zsmax and ccmax are calculated here, after keepgoing has been set to False
            '''
            cctck = splrep(zs, cc, s=0, k=3)
            zsmax = fsolve(lambda z: splev(z, cctck, 1), (zs[zmax],))[0]
            ccmax = splev(zsmax, cctck)

            # An alternative way to compute the uncertainly, the the one below is preferred
            #zlo = fsolve(lambda z: splev(z,cctck)-(ccmax-0.5),(zs[zmax-3],))[0]
            #zhi = fsolve(lambda z: splev(z,cctck)-(ccmax-0.5),(zs[zmax+3],))[0]
            #dzlo = zs[zmax]-zlo
            #dzhi = zhi-zs[zmax]
            # print "         +%.6f" % (dzhi)
            # print "z = %.6f" % (zs[zmax])
            # print "         -%.6f" % (dzlo)

            cc2tck = splrep(zs, cc-(ccmax-0.5), s=0, k=3)
            try:
                ccroots = sproot(cc2tck, mest=2*len(cc))
                if len(ccroots) < 2:
                    shrnk_zs = zs[zmax - 10: zmax + 10]
                    shrnk_cc = cc[zmax - 10: zmax + 10] - ccmax
                    shrnk_cc2tck = splrep(shrnk_zs, shrnk_cc, s=0, k=3)
                    ccroots = sproot(shrnk_cc2tck, mest=2*len(shrnk_cc))
            except(TypeError):
                ccroots = [0]

            self.zp = zp
            # Here we estimate our uncertainly in z
            if len(ccroots) >= 2:
                zlo, zhi = np.sort(
                    np.take(ccroots, np.argsort(abs(ccroots-zsmax)))[:2])
                dzlo, dzhi = zsmax-zlo, np.abs(zhi-zsmax)
                print("         +%.6f" % (dzhi))
                print("z = %.6f" % (zsmax))
                print("         -%.6f" % (dzlo))
                self.zsmax = zsmax
                self.ccmax = ccmax
                self.dzhi = dzhi
                self.dzlo = dzlo
            else:
                zlo = fsolve(lambda z: splev(z, cctck) -
                             (ccmax-0.5), (zs[zmax-3],))[0]
                zhi = fsolve(lambda z: splev(z, cctck) -
                             (ccmax-0.5), (zs[zmax+3],))[0]
                dzlo = zs[zmax]-zlo
                dzhi = zhi-zs[zmax]
                print("Fewer than two roots: solved independently")
                print("         +%.6f" % (dzhi))
                print("z = %.6f" % (zsmax))
                print("         -%.6f" % (dzlo))
                self.zsmax = zsmax
                self.ccmax = ccmax
                self.dzhi = dzhi
                self.dzlo = dzlo

                #self.zsmax = 1
                #self.ccmax = 1
                #self.dzhi = 1
                #self.dzlo = 1
                # print "Error, fewer than 2 roots"

            # Here we compute chi2 and reduce chi2
            Tzs = T.Redshift(zsmax, G.wavelength/(1+zsmax))
            Tzm = np.add.reduce(Tzs*Gw)/np.sum(Gw)
            Tzs = Tzs/Tzm * Gn
            Tps = FitPoly(Tzs, G.noise, G.continuum, G.wavelength, o, w1, w2)
            Gm = 1.*G.good * \
                np.greater_equal(G.wavelength, w1) * \
                np.less_equal(G.wavelength, w2)
            Gr = (G.spec-Gp)*Gm
            Trs = (Tzs-Tps)*Gm
            self.chi2 = np.sum(divz(Trs-Gr, G.noise)**2*Gm)
            self.rchi2 = self.chi2/(np.sum(np.greater(Gm, 0))-1)
            print("Template %d,  Chi^2=%.2f   RChi^2=%.2f  CCmax=%.2f  zsmax=%.6f" % (
                t, self.chi2, self.rchi2, ccmax, zsmax))

'''
Shell of the code - this begins the cross correlation object by object
Plot() - sets up the object as a class that can be correlated
doCC() - performs the cross correlation
createPlot() - opens the GUI, which allows the user to modify and save the results of the correlation
'''
if tcorr:
    tcstr = 'cor_'
else:
    tcstr = ''
if markbads:
    bads = np.loadtxt(badfile, unpack=True)
if haselist:
    emitters = np.loadtxt(emissfile, unpack=True)

if readunsure == 1:
    data = np.genfromtxt(outfile, dtype=None, names=True)
    dataunsure = data[data['Unsure'] == 1]
    dataunsure = dataunsure[dataunsure['Star'] == 0]
    dataunsure = dataunsure[dataunsure['Revisit'] == 0]
    dataunsure = dataunsure[dataunsure['Note'] == 0]
    dataunsure = dataunsure[dataunsure['Unusable'] == 0]
    #dataunsure = dataunsure[dataunsure['SN'] > 5]
    tarr = [i for i in reversed(dataunsure[np.argsort(dataunsure['SN'])])]
    imarr = []
    for i in range(len(tarr)):
        imarr.append(imloc+tarr[i]['ImageName'])
    rerun = True
elif objid:
    imarr = glob.glob(imloc + tcstr + objid + '_' + imname)
    if not imarr:
        sys.exit('Cannot find ' + imloc + tcstr + objid + '_' + imname)
    rerun = True
else:
    rerun = False
    outdata = np.genfromtxt(outfile, dtype=None, names=True)
    imarr = glob.glob(imloc+tcstr+'??????_' + imname)
    if not imarr:
        sys.exit('No image files found of form ' +
                 imloc + tcstr+'??????_' + imname)

max_k = len(imarr)
for k in range(max_k):
    j = imarr[k]
    if not rerun:
        if outdata[k]['OBJID'] != 0:
            continue  # Already been worked
    im = Plot(j)
    if markbads:
        if float(im.image[:6]) in bads:
            print('Marked Bad!')
            im.baddata = 1
    if haselist:
        if float(im.image[:6]) in emitters:
            print('Assuming Emission')
            im.eclip = False
    im.doCC(im.image, im.eclip, 0)
    print('Image {0} / {1}'.format(k+1, max_k))
    print(imarr[k])
    im.createPlot()

reason: import numpy,from scipy,from astropy


text:
import numpy as np
import keras.backend.tensorflow_backend as backend
from keras.models import Sequential
from keras.layers import Dense, Dropout, Conv2D, MaxPooling2D, Activation, Flatten
from keras.optimizers import Adam
from keras.callbacks import TensorBoard
import tensorflow as tf
from collections import deque
import time
import random
from game import VanilaGame,Snake,Food,Board
from tqdm import tqdm
import os

import cv2
tf.executing_eagerly()
DISCOUNT = 0.99
REPLAY_MEMORY_SIZE = 50_000  # How many last steps to keep for model training
MIN_REPLAY_MEMORY_SIZE = 10_000  # Minimum number of steps in a memory to start training
MINIBATCH_SIZE = 32  # How many steps (samples) to use for training
UPDATE_TARGET_EVERY = 50  # Terminal states (end of episodes)
MODEL_NAME = 'BabyConnerPt2.0.1'
MIN_REWARD = -200  # For model save
MEMORY_FRACTION = 0.20

# Environment settings
EPISODES = 20_000

# Exploration settings
epsilon = 1  # not a constant, going to be decayed
EPSILON_DECAY = 0.99975
MIN_EPSILON = 0.001

#  Stats settings
AGGREGATE_STATS_EVERY = 10  # episodes
SHOW_PREVIEW = True

env = VanilaGame(300,300,30)

# For stats
ep_rewards = [-200]

# For more repetitive results
random.seed(1)
np.random.seed(1)
tf.compat.v1.set_random_seed(1)# tf 1.0
#tf.random.set_seed(1)

# Memory fraction, used mostly when trai8ning multiple agents
#gpu_options = tf.GPUOptions(per_process_gpu_memory_fraction=MEMORY_FRACTION)
#backend.set_session(tf.Session(config=tf.ConfigProto(gpu_options=gpu_options)))

# Create models folder
if not os.path.isdir('models'):
    os.makedirs('models')
# Own Tensorboard class
class ModifiedTensorBoard(TensorBoard):

    # Overriding init to set initial step and writer (we want one log file for all .fit() calls)
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        self.step = 1
        self.model = None
        self.TB_graph = tf.compat.v1.Graph()
        with self.TB_graph.as_default():
            self.writer = tf.summary.create_file_writer(self.log_dir, flush_millis=5000)
            self.writer.set_as_default()
            self.all_summary_ops = tf.compat.v1.summary.all_v2_summary_ops()
        self.TB_sess = tf.compat.v1.InteractiveSession(graph=self.TB_graph)
        self.TB_sess.run(self.writer.init())

    # Overriding this method to stop creating default log writer
    def set_model(self, model):
        self.model = model
        self._train_dir = self.log_dir + '\\train'

    # Overrided, saves logs with our step number
    # (otherwise every .fit() will start writing from 0th step)
    def on_epoch_end(self, epoch, logs=None):
        self.update_stats(**logs)

    # Overrided
    # We train for one batch only, no need to save anything at epoch end
    def on_batch_end(self, batch, logs=None):
        pass

    def on_train_begin(self, logs=None):
        pass

    # Overrided, so won't close writer
    def on_train_end(self, _):
        pass

    # added for performance?
    def on_train_batch_end(self, _, __):
        pass

    # Custom method for saving own metrics
    # Creates writer, writes custom metrics and closes writer
    def update_stats(self, **stats):
        self._write_logs(stats, self.step)

    def _write_logs(self, logs, index):
        for name, value in logs.items():
            self.TB_sess.run(self.all_summary_ops)
            if self.model is not None:
                name = f'{name}_{self.model.name}'
            self.TB_sess.run(tf.summary.scalar(name, value, step=index))
        self.model = None

class DQNAgent:
    def __init__(self):
        #main model gets trained every step
        self.tensorboard = ModifiedTensorBoard(log_dir = f"logs/{MODEL_NAME}-{int(time.time())}")
        tf.executing_eagerly()
        self.model = self.create_model()

        #target model this is what we predict evbery step
        self.target_model = self.create_model()
        self.target_model.set_weights(self.model.get_weights())

        self.replay_memory = deque(maxlen=REPLAY_MEMORY_SIZE)

        self.target_update_counter = 0

    def create_model(self):
        model = Sequential()
        model.add(Conv2D(256,(3,3), input_shape=env.OBSERVATION_SPACE_VALUES))
        model.add(Activation("relu"))
        model.add(MaxPooling2D(2,2))
        model.add(Dropout(0.2))

        model.add(Conv2D(256,(3,3)))
        model.add(Activation("relu"))
        model.add(MaxPooling2D(2,2))
        model.add(Dropout(0.2))

        model.add(Flatten())
        model.add(Dense(64))
        model.add(Dense(env.ACTION_SPACE_SIZE,activation="linear"))
        model.compile(loss='mse', optimizer=Adam(lr=0.001),metrics=['accuracy'])
        return model

    def update_replay_memory(self, transition):
        self.replay_memory.append(transition)

    def get_qs(self,state):
        return self.model.predict(np.array(state).reshape(-1,*state.shape)/255)[0]

    def train(self,terminal_state,step):
        if len(self.replay_memory) < MIN_REPLAY_MEMORY_SIZE:
            return

        minibatch = random.sample(self.replay_memory,MINIBATCH_SIZE)

        current_states = np.array([transition[0] for transition in minibatch])/255
        current_qs_list = self.model.predict(current_states)

        new_current_states = np.array([transition[3] for transition in minibatch ])/255
        future_qs_list = self.target_model.predict(new_current_states)

        X = []
        y = []

        for index, (current_state, action, reward, new_current_state, done) in enumerate(minibatch):
            if not done:
                max_future_q = np.max(future_qs_list[index])
                new_q = reward + DISCOUNT * max_future_q
            else:
                new_q = reward
            current_qs = current_qs_list[index]

            current_qs[action] = new_q
            X.append(current_state)
            y.append(current_qs)
        self.model.fit(np.array(X)/255,np.array(y), batch_size=MINIBATCH_SIZE,
        verbose=0,shuffle=False,callbacks=[self.tensorboard] if terminal_state else None)

        if terminal_state:
            self.target_update_counter += 1

        if self.target_update_counter > UPDATE_TARGET_EVERY:
            self.target_model.set_weights(self.model.get_weights())
            self.target_update_counter = 0

agent = DQNAgent()

for episode in tqdm(range(1,EPISODES+1),ascii=True, unit="episode"):
    agent.tensorboard.step = episode
    episode_reward = 0
    step = 1
    current_state = env.reset()

    done = False
    while not done:
        if np.random.random() >epsilon:
            action = np.argmax(agent.get_qs(current_state))
        else:
            action = np.random.randint(0,env.ACTION_SPACE_SIZE)
        new_state, reward, done = env.step(action)

        episode_reward += reward

        if SHOW_PREVIEW and not episode % AGGREGATE_STATS_EVERY:
            env.render()
        agent.update_replay_memory((current_state,action, reward, new_state, done))
        agent.train(done,step)

        current_state = new_state
        step += 1
        # Append episode reward to a list and log stats (every given number of episodes)
    ep_rewards.append(episode_reward)

    if not episode % AGGREGATE_STATS_EVERY or episode == 1:
        average_reward = sum(ep_rewards[-AGGREGATE_STATS_EVERY:])/len(ep_rewards[-AGGREGATE_STATS_EVERY:])
        min_reward = min(ep_rewards[-AGGREGATE_STATS_EVERY:])
        max_reward = max(ep_rewards[-AGGREGATE_STATS_EVERY:])
        mem_size = len(agent.replay_memory)
        agent.tensorboard.update_stats(reward_avg=average_reward, reward_min=min_reward, reward_max=max_reward, epsilon=epsilon,replay_mem_size= mem_size)

        # Save model, but only when min reward is greater or equal a set value
        if min_reward >= MIN_REWARD:
            agent.model.save(f'models/{MODEL_NAME}__{max_reward:_>7.2f}max_{average_reward:_>7.2f}avg_{min_reward:_>7.2f}min__{int(time.time())}.model')

    # Decay epsilon
    if epsilon > MIN_EPSILON:
        epsilon *= EPSILON_DECAY
        epsilon = max(MIN_EPSILON, epsilon)

ext: py

max_stars_repo_name: allepicondor/SnakePython

max_issues_repo_name: allepicondor/SnakePython

max_forks_repo_name: allepicondor/SnakePython

reason: import numpy


text:
#!/usr/bin/env python3
#  Copyright 2017 Google Inc. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS-IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from typing import Dict, Iterable

import networkx as nx

from _py2tmp.ir0 import ir

def compute_template_dependency_graph(template_defns: Iterable[ir.TemplateDefn], template_defn_by_name: Dict[str, ir.TemplateDefn]):
    template_dependency_graph = nx.DiGraph()
    for template_defn in template_defns:
        template_dependency_graph.add_node(template_defn.name)

        for identifier in template_defn.referenced_identifiers:
            if identifier in template_defn_by_name.keys():
                template_dependency_graph.add_edge(template_defn.name, identifier)
    return template_dependency_graph

max_stars_repo_name: google/tmppy

max_issues_repo_name: google/tmppy

max_forks_repo_name: google/tmppy

reason: import networkx


text:
import pycxsimulator
from pylab import *

import networkx as nx

def initialize():
    global g
    g = nx.convert_node_labels_to_integers(nx.grid_2d_graph(10, 10))
    g.pos = nx.spring_layout(g, k = 0.1)
    g.count = 0

def observe():
    global g
    cla()
    nx.draw(g, pos = g.pos)

def update():
    global g
    g.count += 1
    if g.count % 20 == 0: # rewiring once in every 20 steps
        nds = list(g.nodes)
        i = choice(nds)
        if g.degree(i) > 0:
            g.remove_edge(i, choice(list(g.neighbors(i))))
            nds.remove(i)
            for j in g.neighbors(i):
                nds.remove(j)
            g.add_edge(i, choice(nds))

    # simulation of node movement
    g.pos = nx.spring_layout(g, pos = g.pos, iterations = 3, k = 0.1)

pycxsimulator.GUI().start(func=[initialize, observe, update])

max_stars_repo_name: RachidStat/PyCX

max_issues_repo_name: RachidStat/PyCX

max_forks_repo_name: RachidStat/PyCX


text:
import catboost as cat
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import time
from config import *
from helper.utils import print_time,random_seed_cpu
import copy

def _plot_importance_(importances, features_name, fig_size=(15, 15), top=50):
    importances = [
        pd.DataFrame(
            {
                "features": importance["Feature Id"],
                "importance_" + str(i): importance["Importances"],
            }
        )
        for i, importance in enumerate(importances)
    ]

    df = pd.DataFrame()
    df["features"] = features_name
    for importance in importances:
        df = df.merge(importance, on="features", how="left")
    df["importance"] = df.drop(["features"], 1).mean(axis=1)
    df.sort_values("importance", inplace=True, ascending=False)
    plt.figure(figsize=fig_size)
    sns.barplot(data=df[:top], y="features", x="importance")

def upsampling(class_to_samples,train,columns,coff=1):
    random_seed_cpu(RANDOM_STATE)
    final_train=[train]
    for _class_ in class_to_samples :
        train_class=train[train.label==_class_]
        for col in columns :
            train_class[col]=np.random.permutation(train_class[col].values)
        train_class=train_class.sample(frac=coff)
        final_train.append(train_class)
    final_train=pd.concat(final_train)
    return final_train

def get_data_set(train, target_name, features, cat_features, weight,other_params,fold_num=None):
    if fold_num != None:
        X_train = train.loc[train.fold != fold_num]
        X_val = train.loc[train.fold == fold_num]
#         print(len(X_train))
#         X_train=upsampling(other_params["class_to_samples"],X_train,other_params["columns"],other_params["coff"])
#         print(len(X_train))
        dtrain = cat.Pool(
            data=X_train[features],
            label=X_train[target_name],
            feature_names=features,
            cat_features=cat_features,
            weight=weight,
        )
        dval = cat.Pool(
            data=X_val[features],
            label=X_val[target_name],
            cat_features=cat_features,
            feature_names=features,
            weight=weight,
        )
    else:
        X_train = train.loc[train.validation == 0]
        X_val = train.loc[train.validation == 1]

        dtrain = cat.Pool(
            data=X_train[features],
            label=X_train[target_name],
            feature_names=features,
            cat_features=cat_features,
            weight=weight,
        )
        dval = cat.Pool(
            data=X_val[features],
            label=X_val[target_name],
            feature_names=features,
            cat_features=cat_features,
            weight=weight,
        )

    return dtrain, dval, X_train, X_val

def cat_single(
    train,
    test,
    target_name,
    features,
    params,
    metric,
    other_params,
    fold_num=None,
    plot_importance=True,
    print_result=True,
    predict_train=True,
):

    dtrain, dval, X_train, X_val = get_data_set(
        train, target_name, features, other_params["cat_features"],other_params["weight"],other_params ,fold_num
    )
    dtest = cat.Pool(
        data=test[features],
        feature_names=features,
        cat_features=other_params["cat_features"],
    )

    model = cat.CatBoost(params)
    model.fit(dtrain, eval_set=dval, use_best_model=True)
    train_pred = (
        model.predict(dtrain, other_params["prediction_type"])
        if predict_train
        else None
    )
    val_pred = model.predict(dval, prediction_type=other_params["prediction_type"])
    test_pred = model.predict(dtest, prediction_type=other_params["prediction_type"])
    if ((other_params["prediction_type"]=="Probability")& (other_params["num_class"]==1)) :
        if predict_train :
            train_pred=train_pred[:,1]
        val_pred=val_pred[:,1]
        test_pred=test_pred[:,1]

    val_score = metric(X_val[target_name], val_pred)
    train_score = metric(X_train[target_name], train_pred) if predict_train else -1

    if print_result:
        print(
            "final train score : {} -validation score: {}".format(
                str(round(train_score, 5)), str(round(val_score, 5))
            )
        )
    importances = [model.get_feature_importance(prettified=True)]

    if plot_importance:
        _plot_importance_(importances, features)
    return train_pred, val_pred, test_pred, train_score, val_score, importances, model

def cat_kfold(
    train, test, target_name, features, params, metric, other_params,
):
    models = {}
    oof_train = np.zeros((len(train), other_params["num_class"]))
    oof_test = np.zeros((len(test), other_params["num_class"]))
    importances = []
    validation_scores = []
    train_scores = []
    for fold in np.sort(train.fold.unique()):

        (
            train_pred,
            val_pred,
            test_pred,
            train_score,
            val_score,
            importance,
            model,
        ) = cat_single(
            train,
            test,
            target_name,
            features,
            params,
            metric=metric,
            other_params=other_params,
            fold_num=fold,
            plot_importance=False,
            print_result=other_params["print_result"],
            predict_train=other_params["predict_train"],
        )

        models[str(fold)] = model

        oof_train[train[train.fold == fold].index, :] += np.reshape(
            val_pred, (-1, other_params["num_class"])
        )
        oof_test += np.reshape(test_pred, (-1, other_params["num_class"]))

        train_scores.append(train_score)
        validation_scores.append(val_score)
        importances.extend(importance)

        if other_params["print_kfold_eval"] :
            print(
                "Iteration : {}  - train score : {} - CV Score : {}".format(
                    str(fold + 1), str(train_score), str(val_score)
                )
            )
            print("=" * 80)
    oof_test /= train.fold.nunique()
    final_train_score = np.mean(train_scores)
    oof_score = metric(train[target_name], oof_train)

    if other_params["print_kfold_eval"] :
        print(
            "ending  training  : train score {} - oof Score {}".format(
                str(final_train_score), str(oof_score)
            )
        )

    if other_params["plot_importance_kfold"]:
        _plot_importance_(importances, features)

    return oof_train, oof_test, final_train_score, oof_score, models

def cat_train(train, test, other_params):
    start_time=time.time()
    if other_params["use_kfold"]:

        oof_train, oof_test, final_train_score, oof_score, models = cat_kfold(
            train,
            test,
            other_params["target_name"],
            other_params["features"],
            other_params["params"],
            other_params["metric"],
            other_params,
        )
        finish_time=time.time()
        if other_params["print_time"] : print_time(finish_time-start_time)
        return oof_train, oof_test, final_train_score, oof_score, models
    else:

        (
            train_pred,
            val_pred,
            test_pred,
            train_score,
            val_score,
            importances,
            model,
        ) = cat_single(
            train,
            test,
            other_params["target_name"],
            other_params["features"],
            other_params["params"],
            other_params["metric"],
            other_params,
            fold_num=None,
            plot_importance=other_params["plot_importance"],
            print_result=True,
            predict_train=other_params["predict_train"],
        )
        finish_time=time.time()
        if other_params["print_time"] : print_time(finish_time-start_time)
        return (train_pred, val_pred, test_pred, train_score, val_score, model)

def cat_train_add_features(train, test, other_params):
    start_time=time.time()
    other_params_copy=copy.deepcopy(other_params)
    best_score=other_params["best_score"]
    best_validation_prediction=[]
    best_test_prediction=[]
    to_keep_features=[]
    results=[]
    for feature in other_params["added_features"] :
        result=pd.Series()
        features_to_test=other_params["features"]+to_keep_features+[feature]
        other_params_copy["features"]=features_to_test
        train_score,validation_score,validation_prediction,test_prediction,models=train_model(train, test, other_params_copy)

        result["feature"]=feature
        result["train_score"]=train_score
        result["validation_score"]=validation_score
        result["best_iteration"]=np.mean([ model.get_best_iteration() for model in models ])

        if other_params["maximise"] :
            if validation_score >= best_score :
                to_keep_features.append(feature)
                result["improved"]=1
                best_score=validation_score
                best_validation_prediction=validation_prediction
                best_test_prediction=test_prediction
            else :
                result["improved"]=0
        else  :
            if validation_score <= best_score :
                to_keep_features.append(feature)
                result["improved"]=1
                best_score=validation_score
                best_validation_prediction=validation_prediction
                best_test_prediction=test_prediction
            else :
                result["improved"]=0
        result["model"]=models
        result["features_len"]=len(features_to_test)

        display(pd.DataFrame([result]))
        results.append(result)
    results=pd.DataFrame(results)
    finish_time=time.time()
    print_time(finish_time-start_time)
    return results ,to_keep_features ,best_validation_prediction,best_test_prediction

def cat_train_remove_features(train, test, other_params):
    start_time=time.time()
    other_params_copy=copy.deepcopy(other_params)
    best_score=other_params_copy["best_score"]
    to_remove_features=[]
    results=[]
    best_validation_prediction=[]
    best_test_prediction=[]
    for feature in other_params["removed_features"] :

        result=pd.Series()
        features_to_test=copy.deepcopy(other_params["features"])
#         features_to_test.remove(feature)
#         for rm_feature in to_remove_features :
#             features_to_test.remove(rm_feature)
#         other_params_copy["features"]=features_to_test
        features_to_test=[feat for feat in features_to_test if not feat.endswith(feature) ]
        for rm_feature in to_remove_features :
            features_to_test=[feat for feat in features_to_test if not feat.endswith(rm_feature) ]
        other_params_copy["features"]=features_to_test

        train_score,validation_score,validation_prediction,test_prediction,models=train_model(train, test, other_params_copy)

        result["feature"]=feature
        result["train_score"]=train_score
        result["validation_score"]=validation_score
        result["best_iteration"]=np.mean([ model.get_best_iteration() for model in models ])

        if other_params["maximise"] :
            if validation_score >= best_score :
                to_remove_features.append(feature)
                result["improved"]=1
                best_score=validation_score
                best_validation_prediction=validation_prediction
                best_test_prediction=test_prediction
            else :
                result["improved"]=0
        else :
            if validation_score <= best_score :
                to_remove_features.append(feature)
                result["improved"]=1
                best_score=validation_score
                best_validation_prediction=validation_prediction
                best_test_prediction=test_prediction
            else :
                result["improved"]=0
        result["features_len"]=len(features_to_test)
        result["model"]=models

        display(pd.DataFrame([result]))
        results.append(result)
    results=pd.DataFrame(results)
    finish_time=time.time()
    print_time(finish_time-start_time)
    return results ,to_remove_features,best_validation_prediction,best_test_prediction

def train_model(train, test, other_params) :
    if other_params["use_kfold"]:
        oof_train, oof_test, final_train_score, oof_score, models = cat_kfold(
            train,
            test,
            other_params["target_name"],
            other_params["features"],
            other_params["params"],
            other_params["metric"],
            other_params,
        )
        return final_train_score ,oof_score,oof_train,oof_test,list(models.values())

    else:

        (
            train_pred,
            val_pred,
            test_pred,
            train_score,
            val_score,
            importances,
            model,
        ) = cat_single(
            train,
            test,
            other_params["target_name"],
            other_params["features"],
            other_params["params"],
            other_params["metric"],
            other_params,
            fold_num=None,
            plot_importance=other_params["plot_importance"],
            print_result=other_params["print_result"],
            predict_train=other_params["predict_train"],
        )
        return train_score,val_score,val_pred,test_pred,[ model]

max_stars_repo_path: Image Classification/Radiant Earth Computer Vision for Crop Detection from Satellite Imagery/#10_Be_carefuk/ICLR_2/helper/catboost.py

max_stars_repo_name: ZindiAfrica/Computer-Vision


text:
import numpy as np
import random
"""
约束条件1 x，y，z三个方向的求和均值为0
无序       random.uniform
对于第i个元素 x**2+y**2+z**2 == 1
"""

TMP = []
for i in range(99999):
    tmp = []
    for i in range(3):
        data =  random.uniform(-6,6)     ### 标注磁矩
        tmp.append(data)
    if  0.001<tmp[0]**2 + tmp[1]**2 + tmp[2]**2 - 6**2 < 0.01 :     ### 标注磁矩
        #print(tmp)
        #print(tmp[0]**2 + tmp[1]**2 + tmp[2]**2)
        TMP.append(tmp)
        #break
    else:
        pass

#######################setp1#########################################
#######################setp2#########################################
testx = []
testy = []
testz = []
for i in range(99999999):
    set1=[]
    for i in range(0,6):    ### 几个元素
        data =  random.uniform(1,len(TMP))
        #print(i)
        set1.append(int(data))
    #print(set1)
    #print("----------------")
    testx = []
    testy = []
    testz = []
    for s in set1:
        testx.append(TMP[s][0])
        testy.append(TMP[s][1])
        testz.append(TMP[s][2])

    if  0.001<abs(np.mean(testx)) + abs(np.mean(testy))+abs(np.mean(testz)) <= 0.05:
        print(np.mean(testx))
        print(np.mean(testy))
        print(np.mean(testz))
            #print(set1)
        print("x",testx)
        print("y",testy)
        print("z",testz)

    else:
            #print("Nope")
        pass

    #break

max_stars_repo_name: suthzx/-

max_issues_repo_name: suthzx/-

max_forks_repo_name: suthzx/-