def main():

    is_transparent = False 
    
    f = open("pi_data.txt","r")
    
    # this is a little different than normal becase of the complex data for the floquet stability
    # multipliers. When we use the "dtype" option we get a single array of tuples so slicing is a
    # little more awkward has to look like data[#][#] to get a single value NOT data[#,#].
    data = pl.genfromtxt(f,comments="e",dtype="complex,complex,float")
   
    eigs1 = pl.array([])
    eigs2 = pl.array([])
    A = pl.array([])
    
    for i,j in enumerate(data):
        eigs1 = pl.append(eigs1,j[0])
        eigs2 = pl.append(eigs2,j[1])
        A = pl.append(A,j[2])

    fig1, ax1 = pl.subplots(2,2,sharex=True)
    ax1[0,0].plot(A,[k.real for k in eigs1],color = "Black")
    ax1[1,0].plot(A,[k.imag for k in eigs1],color = "Black")
    ax1[0,1].plot(A,[k.real for k in eigs2],color = "Black")
    ax1[1,1].plot(A,[k.imag for k in eigs2],color = "Black")

    ax1[0,0].set_ylabel("Re[$\lambda_1$]",fontsize=25)
    ax1[1,0].set_ylabel("Im[$\lambda_1$]",fontsize=25)
    ax1[0,1].set_ylabel("Re[$\lambda_2$]",fontsize=25)
    ax1[1,1].set_ylabel("Im[$\lambda_2$]",fontsize=25)
    ax1[1,0].set_xlabel("$A$",fontsize=25)
    ax1[1,1].set_xlabel("$A$",fontsize=25)
    fig1.tight_layout()
    fig1.savefig("paper_A_vs_eigs.png",dpi=300,transparent=is_transparent)
    os.system("open paper_A_vs_eigs.png")

    def f(self,x,t):
        # This one line is different than in ECclass.py
        N = 2
        xdot = pl.array([])

        # modulus the x for periodicity.
        x[N:2*N]= x[N:2*N]%self.d
        # HERE ---->> 1Dify
        for i in range(N):
            temp = 0.0
            for j in range(N):
                if i == j:
                    continue
                #repulsive x interparticle force of j on i
                temp += self.qq*(x[N+i]-x[N+j])/(pl.sqrt((x[N+i]-x[N+j])**2)**3)
                # All of the forces coming from the 'same' paricle but from other 'cells' due to the
                # periodic contrains can be wraped up in a sum that converges to an aswer that can
                # be expressed in terms of polygamma functions (se pg 92 of notebook).
                # Note on the sign (xi-xj or xj-xi). Changing the sign of the xi-xj term (i.e. which
                # particle are we considering forces on) changes the direction of the force
                # apropriately.
                temp += self.qq*(polygamma(1,(self.d+x[N+i]-x[N+j])/self.d)-polygamma(1,1.0-((x[N+i]-x[N+j])/self.d)))/(self.d**2)
            # periodic force on particle i
            temp += self.A*pl.sin(x[N+i])*pl.cos(t)
            temp -= self.beta*x[i]
            xdot = pl.append(xdot,temp)
        for i in range(N):
            xdot = pl.append(xdot,x[i])
        return xdot

def euler_1(top):
    '''
    inputs:
        top = maximum value considered for suming
    outputs:
        sum = sum of all numbers below top that are mulitples of 3
        or five
    '''
    
    #define varibles
    mfive = []
    mthree = []
    mboth = []
    
    #deine win condition
    win = 0
    
    for i in range(0, top):
        
        if(i%3==0 and i%5==0):
            mboth = pylab.append(i, mboth)
            
        elif(i%3==0):
            mthree = pylab.append(i, mthree)
            
        elif(i%5==0):
            mfive = pylab.append(i, mfive)
            
    
    win = sum(mboth) + sum(mfive) + sum(mthree)
    
                            
    print win

def get_first_init(x0,epsilon,N):
    x_new = pl.copy(x0)
    print('getting the first initial condition')
    print('fiducial initial: '+str(x0))
    # multi particle array layout [nth particle v, (n-1)th particle v , ..., 0th v, nth particle x, x, ... , 0th particle x]
    
    # we will use a change of coordinates to get the location of the particle relative to x0. First
    # we just find some random point a distace epsilon from the origin.
    # need 2DN random angles 
    angle_arr = pl.array([])
    purturbs = pl.array([])

    # This is just an n-sphere 
    for i in range(2*N):
        angle_arr = pl.append(angle_arr,random.random()*2.0*pl.pi)
        cur_purt = epsilon
        for a,b in enumerate(angle_arr[:-1]):
            cur_purt *= pl.sin(b)
        if i == (2*N-1):
            cur_purt = pl.sin(angle_arr[i])
        else:
            cur_purt = pl.cos(angle_arr[i])

        purturbs = pl.append(purturbs,cur_purt)

    print('sqrt of sum of squars should be epsilon -> is it? --> ' +str(pl.sqrt(pl.dot(purturbs,purturbs))))
    print('len(purturbs) == 2N ? ' +str(len(purturbs)==(2*N)))

    return x_new+purturbs

def f(x, t):
    # for now masses just = 1.0

    # the 4.0 only works for 2D
    N = len(x) / 4
    xdot = pl.array([])

    for i in range(N):
        temp = 0.0
        for j in range(N):
            if i == j:
                continue
            temp += -(x[2 * N + i] - x[2 * N + j]) / (
                pl.sqrt((x[2 * N + i] - x[2 * N + j]) ** 2 + (x[3 * N + i] - x[3 * N + j]) ** 2) ** 3
            )
        xdot = pl.append(xdot, temp)
    for i in range(N):
        temp = 0.0
        for j in range(N):
            if i == j:
                continue
            temp += -(x[3 * N + i] - x[3 * N + j]) / (
                pl.sqrt((x[2 * N + i] - x[2 * N + j]) ** 2 + (x[3 * N + i] - x[3 * N + j]) ** 2) ** 3
            )
        xdot = pl.append(xdot, temp)
    for i in range(N):
        xdot = pl.append(xdot, x[i])
    for i in range(N):
        xdot = pl.append(xdot, x[N + i])

    print("len xdot is: " + str(len(xdot)))
    return xdot

def read_file(datafile,obsv):

    x=[];y=[];t=[];pt=[];z=[];px=[];py=[] #l is the t coordinate, z is the obs location
    for o in range(obsv):
        x.append([]);y.append([]);t.append([])
        px.append([]);py.append([]);pt.append([])
    with open(datafile) as f:
        data = f.read()
    data=data.split('\n')
    print "data length:", (len(data))
    n=0;i=0
    nparts= (len(data)-8)/obsv

    for row in data:
        if row.startswith('@') or row.startswith('*') \
            or row.startswith('$') or row.startswith('#')\
            or len(row)==0: continue
        else:
            i+=1
            w=" ".join(row.split())
            s=w.split()
            x[n].append(float(s[2]))
            px[n].append(float(s[3]))
            y[n].append(float(s[4]))
            py[n].append(float(s[5]))
            t[n].append(float(s[6]))
            pt[n].append(float(s[7]))
            if i%nparts==0:
                z.append(float(s[8]))
                i=0;n+=1

    return x,px,y,py,t,pt,z

def main():
    f = open("final_position.txt","r")
    data = pl.genfromtxt(f,comments = "L")
    
    # need to get every other
    x = pl.array([])
    y = pl.array([])
    for i,j in enumerate(data[:-7,2]):
        if i%4 == 0:
            x = pl.append(x,data[i,4])
            y = pl.append(y,j)
    
    print(x)
    print(y)
    fit = np.polyfit(x,y,2)

    print(fit)
    
    #fited = fit[0]+fit[1]*x + fit[2]*x**2
    fited = np.poly1d(fit)
    print(fited)

    #pl.plot(pl.append(x,[.262,.264,.266]),fited(pl.append(x,[.262,.264,.266])),color="black")
    pl.scatter(x,y,color = "black")
    pl.xlabel("$A$",fontsize="30")
    pl.ylabel("$x$",fontsize="30")
    pl.savefig("fin_pts.png",transparent=True,dpi=300)
    
    os.system("open fin_pts.png")

	def poly(self,step):
		# step is the radian resolution
		cA=self.cA
		self.sort()
		st=False
		fin=False
		cdi=0
		alpha=0
		x0=self.cA[0].c+point(self.cA[0].r,0)
		xc=x0
		pl=[]
		while (alpha<(2*math.pi)):
			if self.side(xc)[-1]==0:
				pl.append(xc)
				st=True
				break
			else:
				xx=num.cos(step)*cA[cdi].r
				yy=num.sin(step)*cA[cdi].r
				xc_t=point(xx,yy)
				xc=xc_t.transform(-cA[cdi].c,-alpha)
				alpha=alpha+step
		if not st:
			return pl
		# pivoting
		pv0=self.get_x0()
		pv=point(pv0.x,pv0.y)
		v0=vec(xc,pv)
		vc=vec(xc,pv)
		spin=0.0
		pv_found=False
		while spin<(2*math.pi):
			if self.side(pv)[0]==self.n:
				pv_found=True
				break
			else:
				vc=vc/2
				if vc.mag()<10*res:
					spin=spin+10*step
					vc=v0.rot(spin)
				pv=xc+vc
		# pivoting finished
		if not pv_found:
			return [xc]
		alpha=num.linspace(0,2*math.pi,int(2*math.pi/step))
		alpha[-1]=math.pi*2
		for a in alpha:
			ray=Ray(pv,a)
			try:
				pc=self.intersect(ray)
				pl.append(pc)
			except geoError:
				print 'Unknown Error in ray-ndisc intersection'
				raise geoError, 'Unknown'
		pl[-1]=pl[0]
		return pl

def callback(msg):

    global record
    global data
    global count

    data = append(data, msg.data)
    # print len(msg.data)
    count += 1

    if count == N_CYCLES:
        signal_1 = data[::2]
        signal_2 = data[1::2]
        # signal_1 = [struct.unpack('B', struct.pack('b',d))[0] for d in data[::2]]
        # signal_2 = [struct.unpack('B', struct.pack('b',d))[0] for d in data[1::2]]

        print len(signal_1)
        print signal_1
        print "number of samples(frames): " + str(len(signal_1))
        print "SAMPLING RATE: " + str(RATE) + "(Hz)"
        print "DURATION: " + str(DURATION) + "(s)"

        helper.plot_from_rawdata(signal_1, signal_2, RATE)

        rospy.signal_shutdown("Recording finished")

def remove_discontinuity(value, xgap=10, ygap=200):
    """
    Remove discontinuity (sudden jump) in a series of values.
    Written by Denis, developed for LLC Fringe Counts data.
    value : list or numpy.array
    xgap  : "width" of index of the list/array to adjust steps
    ygap  : threshold value to detect discontinuity
    """
    difflist = pl.diff(value)
    discont_index = pl.find(abs(difflist) > ygap)

    if len(discont_index) == 0:
        return value
    else:
        discont_index = pl.append(discont_index, len(difflist))

    # find indice at discontinuities
    discont = {"start": [], "end": []}
    qstart = discont_index[0]
    for i in range(len(discont_index) - 1):
        if discont_index[i + 1] - discont_index[i] > xgap:
            qend = discont_index[i]
            discont["start"].append(qstart - xgap)
            discont["end"].append(qend + xgap)
            qstart = discont_index[i + 1]

    # add offsets at discontinuities
    result = pl.array(value)
    for i in range(len(discont["end"])):
        result[0 : discont["start"][i]] += result[discont["end"][i]] - result[discont["start"][i]]

    # remove the median
    result = result - pl.median(result)
    return result

def f(filename, theClass=1):
    fs, data = wavfile.read(filename)  # load the data
    # b=[(ele/2**8.)*2-1 for ele in data] # this is 8-bit track, b is now normalized on [-1,1)
    print "Sample rates is: "
    print fs
    X = stft(data, fs, 256.0 / fs, 256.0 / fs)
    X = X[:, 0 : (X.shape[1] / 2)]
    shortTimeFFT = scipy.absolute(X.T)
    shortTimeFFT = scipy.log10(shortTimeFFT)

    # Plot the magnitude spectrogram.
    pylab.figure()
    pylab.imshow(shortTimeFFT, origin="lower", aspect="auto", interpolation="nearest")
    pylab.xlabel("Time")
    pylab.ylabel("Frequency")
    savefig(filename + "SFFT.png", bbox_inches="tight")

    features = mean(shortTimeFFT, axis=1)
    pylab.figure()
    pylab.plot(features, "r")
    savefig(filename + "AFFT.png", bbox_inches="tight")

    with open(filename + ".csv", "w") as fp:
        a = csv.writer(fp, delimiter=",")
        row = pylab.transpose(features)
        row = pylab.append(row, theClass)
        a.writerow(row)

    def __parseOD(self, ll):
        '''OD data lines parsing method'''
        ll = [float(x) for x in ll]

        # Add the current time
        self.time.append(ll[0])
        numRep = 1
        prevClone = ""
        prevCond = ""
        if self.numClones == 1:
            totalReps = len(self.conditionsNU) / self.numConditions

        for idx, od in enumerate(ll[1:]):
            clone = self.clonesNU[idx]
            source = self.sourcesNU[idx]
            condition = self.conditionsNU[idx]

            # Check which clone + replicate we are observing
            if self.numClones == 1:
                numRep = (idx % totalReps) + 1

            else:
                numRep = numRep + 1 if (clone == prevClone and
                                        condition == prevCond) else 1
            prevClone = clone
            prevCond = condition

            # Append OD reading to array
            self.dataHash[clone][numRep][source][condition]['od'] =\
                py.append(self.dataHash[clone][numRep][source]
                          [condition]['od'], od)

def readDatDirectory(key, directory):
    global stats
    #Don't read data in if it's already read
    if not key in DATA["mean"]:
        data = defaultdict(array)

        #Process the dat files
        for datfile in glob.glob(directory + "/*.dat"):
            fileHandle = open(datfile, 'rb')
            keys, dataDict = csvExtractAllCols(fileHandle)
            stats = union(stats, keys)
            for aKey in keys:
                if not aKey in data:
                    data[aKey] = reshape(array(dataDict[aKey]),
                                         (1, len(dataDict[aKey])))
                else:
                    data[aKey] = append(data[aKey],
                                        reshape(array(dataDict[aKey]),
                                                (1, len(dataDict[aKey]))),
                                        axis=0)

        #Process the div files'
        for datfile in glob.glob(directory + "/*.div"):
            fileHandle = open(datfile, 'rb')
            keys, dataDict = csvExtractAllCols(fileHandle)
            stats = union(stats, keys)
            for aKey in keys:
                if not aKey in data:
                    data[aKey] = reshape(array(dataDict[aKey]),
                                         (1, len(dataDict[aKey])))
                else:
                    data[aKey] = append(data[aKey],
                                        reshape(array(dataDict[aKey]),
                                                (1, len(dataDict[aKey]))),
                                        axis=0)

        #Iterate through the stats and calculate mean/standard deviation
        for aKey in stats:
            if aKey in data:
                DATA["mean"][key][aKey] = mean(data[aKey], axis=0)
                DATA["median"][key][aKey] = median(data[aKey], axis=0)
                DATA["std"][key][aKey] = std(data[aKey], axis=0)
                DATA["ste"][key][aKey] = std(data[aKey], axis=0)/ sqrt(len(data[aKey]))
                DATA["min"][key][aKey] = mean(data[aKey], axis=0)-amin(data[aKey], axis=0)
                DATA["max"][key][aKey] = amax(data[aKey], axis=0)-mean(data[aKey], axis=0)
                DATA["actual"][key][aKey] = data[aKey]

def xytToEcef (lat, long, height, bearing, radius):

    x = (radius + height) * cos (lat) * cos (long)
    y = (radius + height) * cos (lat) * sin (long)
    z = (radius + height) * sin (lat)

    Rtmp = rotecef (lat, long);

    R = coord_xfms.rotz (bearing).dot (Rtmp)

    rph = coord_xfms.rot2rph (R.transpose ());

    pose = array([])
    pose = append (pose, array ([x, y, z]))
    pose = append (pose, rph)

    return pose

def get_poin(sol,dt):
    poin = pl.array([])
    for i in range(len(sol)):
        if(((i*dt)%(2.0*pl.pi))<=dt):
            poin = pl.append(poin,sol[i,:])
    poin = poin.reshape(-1,4)
    # flip order of array
    #poin = poin[::-1,:]
    return poin

 def DiscreteSum(self,lower,upper):
     """
     Perform the rectangular integration, starting
     and ending at the values passed to it (0-->N for rect,
     1-->N for trap, etc...).
     """
     for i in range(int(lower),int(upper)):
         self.I += self.f(self.a+i*self.dx)*self.dx
         self.plot = pl.append(self.plot, self.I)
     return self.I, self.plot

def get_fit(which):
    f = open("final_position.txt","r")
    data = pl.genfromtxt(f,comments = "L")

    if which=="x":
        datnum = 2
    if which=="vx":
        datnum = 0

    x = pl.array([])
    y = pl.array([])
    for i,j in enumerate(data[:-7,datnum]):
        if i%2 == 0:
            x = pl.append(x,data[i,4])
            y = pl.append(y,j)

    fit = pl.polyfit(x,y,2)

    fitted = pl.poly1d(fit)
    return fitted 

 def parser(self, argument):
     """ Parse data from the data source and fills the values array
     """
     while True:
         buffer = self.handler.serialHandler()
         if buffer <> None:
             sample = self.unpackSample(buffer)
             if sample <> None:
                 self.values = pylab.append(self.values, pylab.array([sample]), axis=0)
         else:
             break

def main():

    # read in the forest as peak yield
    args = parseCMD()
    fileName = args.fileN
    forest = pl.loadtxt(fileName)

    # determine which sites were left unTreed at maximum yield
    sitesLeft = pl.array([])
    for i in xrange(forest.size):
        if forest[i]!=1:
            sitesLeft = pl.append(sitesLeft, i)
    sitesLeft = sitesLeft[::-1]

    # determine intervals and sort them in order of size
    intervals = pl.array([forest.size-1-sitesLeft[0]])
    for i in xrange(1,sitesLeft.size-1):
        sizeInt = sitesLeft[i]-sitesLeft[i+1]
        intervals = pl.append(intervals, sizeInt)
    intervals = pl.append(intervals, sitesLeft[-1])
    intervals = pl.sort(intervals)

    # assign each interval a number
    nums = pl.arange(1,sitesLeft.size+1)

    # main plot
    fig1 = pl.figure(1)
    ax = fig1.add_subplot(111)
    pl.ylabel('Interval Size', fontsize=20)
    pl.xlabel('(Sorted) Inverval Number', fontsize=20)

    # loop over and plot each file we find
    newt = getD(fileName)
    
    ax.plot(nums,intervals,label='(D=%s)'%(newt), marker='o', linewidth=0,
            markerfacecolor='None', markeredgecolor='Navy') 
 
    # put labels into legend
    ax.legend(loc='upper left',shadow=True)
    
    pl.show()

 def trap(self):
     """
     Trapezoidal numerical quadrature. Returns the 
     value of the definite integral as well as an array for
     plotting of the instantaneous integral values over the 
     integration variable.
     """
     self.I += (self.dx/2.0)*self.f(self.a)
     self.plot = pl.append(self.plot,self.I)
     self.I,self.plot = self.DiscreteSum(1,self.N)
     self.I += (self.dx/2.0)*self.f(self.b)
     return self.I, self.plot