Trying to Fit a Function to a Specific Region of Data

x = s.t #time
    data = s.firstharmonic #harmonics (y values)
    
    c = (diff(sign(diff(data))) < 0).nonzero()[0] + 1 # local max - finding the peaks of the Harmonics
	
    timemax = [] #list for time values (x values)
    maxh = [] #list for peak Harmonics (y values)
    
    for i in c:
        timemax.append(x[i])
        maxh.append(data[i])
    
    maxdata = max(maxh) #This is the maximum value/peak in the list
 
    def cutoff(maxh):
        for i in range(len(maxh)-1):
            if not maxh[i] > maxh[i+1]:
		           startindex = i
                    startvalue =  maxh[i]
		    
		    return startindex

  polyfunc = polyfit(s.t,s.firstharmonic,1) #1st Order polynomial
    yfit = polyval(polyfunc,s.t)

for i in range(len(maxh)-1):  
		if not maxh[i] > maxh[i+1]:

x = s.t #time
    data = s.firstharmonic #harmonics y values
    
    c = (diff(sign(diff(data))) < 0).nonzero()[0] + 1 # local max peak finder
	
    timemax = []
    maxh = []
    rmaxh = [] #empty list to fill with reverse maxh values
    
    for i in c:
        timemax.append(x[i])
        maxh.append(data[i])
    
    #print maxh
    
    for i in reversed(maxh):
        print i
        rmaxh.append(i)
    #print rmaxh
    
    def crap(rmaxh):             
	    for j in range(len(rmaxh)-1):  
		if not rmaxh[j] > rmaxh[j+1]:
		    sliceindex = j
                    slicevalue = rmaxh[j]
		    
		    return sliceindex
    
    sliceindex = crap(rmaxh)   
    
    def twostream(timemax, rmaxh):
        timestream = timemax[: sliceindex]
        timenoisy = timemax[sliceindex :]
        rmaxhstream = rmaxh[: sliceindex]
        rmaxhnoisy= rmaxh[sliceindex:]  
        
        print timestream
        
        print timenoisy
  
            
        #coefs = polyfit(timestream, maxhstream,1)
        #straightfit = polyval(polyfunc,timestream)
        #return timestream, straightfit
    #timestream, straightfit, = twostream(timemax, log(rmaxh))

#!/usr/bin/env python
#
# Electrostatic PIC code in a 1D cyclic domain

from numpy import *
from scipy import *

import matplotlib.pyplot as plt # Matplotlib plotting library

try:
    import matplotlib.gridspec as gridspec  # For plot layout grid
    got_gridspec = True
except:
    got_gridspec = False

# Need an FFT routine, either from SciPy or NumPy
try:
    from scipy.fftpack import fft, ifft
except:
    # No SciPy FFT routine. Import NumPy routine instead
    from fft import fft, ifft

def rk4step(f, y0, dt, args=()):
    """ Takes a single step using RK4 method """
    k1 = f(y0, *args)
    k2 = f(y0 + 0.5*dt*k1, *args)
    k3 = f(y0 + 0.5*dt*k2, *args)
    k4 = f(y0 + dt*k3, *args)

    return y0 + (k1 + 2.*k2 + 2.*k3 + k4)*dt / 6.

def calc_density(position, ncells, L):
    """ Calculate charge density given particle positions
    
    Input
      position  - Array of positions, one for each particle
                  assumed to be between 0 and L
      ncells    - Number of cells
      L         - Length of the domain

    Output
      density   - contains 1 if evenly distributed
    """
    # This is a crude method and could be made more efficient
    
    density = zeros([ncells])
    nparticles = len(position)
    
    dx = L / ncells       # Uniform cell spacing
    for p in position / dx:    # Loop over all the particles, converting position into a cell number
        plower = int(p)        # Cell to the left (rounding down)
        offset = p - plower    # Offset from the left
        density[plower] += 1. - offset
        density[(plower + 1) % ncells] += offset
    # nparticles now distributed amongst ncells
    density *= float(ncells) / float(nparticles)  # Make average density equal to 1
    return density

def periodic_interp(y, x):
    """
    Linear interpolation of a periodic array y at index x
    
    Input

    y - Array of values to be interpolated
    x - Index where result required. Can be an array of values
    
    Output
    
    y[x] with non-integer x
    """
    ny = len(y)
    if len(x) > 1:
        y = array(y) # Make sure it's a NumPy array for array indexing
    xl = floor(x).astype(int) # Left index
    dx = x - xl
    xl = ((xl % ny) + ny) % ny  # Ensures between 0 and ny-1 inclusive
    return y[xl]*(1. - dx) + y[(xl+1)%ny]*dx

def fft_integrate(y):
    """ Integrate a periodic function using FFTs
    """
    n = len(y) # Get the length of y
    
    f = fft(y) # Take FFT
    # Result is in standard layout with positive frequencies first then negative
    # n even: [ f(0), f(1), ... f(n/2), f(1-n/2) ... f(-1) ] Only want the first and the last
    # n odd:  [ f(0), f(1), ... f((n-1)/2), f(-(n-1)/2) ... f(-1) ]
    
    """Just interested in the first harmonic for the fastest growing mode"""
    
   
    if n % 2 == 0: # If an even number of points (if not then n is odd)
        k = concatenate( (arange(0, n/2+1), arange(1-n/2, 0)) )
    else:
        k = concatenate( (arange(0, (n-1)/2+1), arange( -(n-1)/2, 0)) )
    k = 2.*pi*k/n
    
    # Modify frequencies by dividing by ik
    f[1] /= (1j * k[1]) 
    f[n-1] /= (1j * k[n-1]) #Have cut out alot of harmonics here
    f[2:n] = 0.0
    f[0] = 0. # Set the arbitrary zero-frequency term to zero
    
    return ifft(f).real # Reverse Fourier Transform
   

def pic(f, ncells, L):
    """ f contains the position and velocity of all particles
    """
    nparticles = len(f)/2     # Two values for each particle
    pos = f[0:nparticles] # Position of each particle
    vel = f[nparticles:]      # Velocity of each particle

    dx = L / float(ncells)    # Cell spacing

    # Ensure that pos is between 0 and L
    pos = ((pos % L) + L) % L
    
    # Calculate number density, normalised so 1 when uniform
    density = calc_density(pos, ncells, L)
    
    # Subtract ion density to get total charge density
    rho = density - 1.
    
    # Calculate electric field
    E = -fft_integrate(rho)*dx
    
    # Interpolate E field at particle locations
    accel = -periodic_interp(E, pos/dx)

    # Put back into a single array
    return concatenate( (vel, accel) )

####################################################################

def run(pos, vel, L, ncells=None, out=[], output_times=linspace(0,20,100), cfl=0.5):
    
    if ncells == None:
        ncells = int(sqrt(len(pos))) # A sensible default

    dx = L / float(ncells)
    
    f = concatenate( (pos, vel) )   # Starting state
    nparticles = len(pos)
    
    time = 0.0
    for tnext in output_times:
        # Advance to tnext
        stepping = True
        while stepping:
            # Maximum distance a particle can move is one cell
            dt = cfl * dx / max(abs(vel))
            if time + dt >= tnext:
                # Next time will hit or exceed required output time
                stepping = False
                dt = tnext - time
            #print "Time: ", time, dt
            f = rk4step(pic, f, dt, args=(ncells, L))
            time += dt
            
        # Extract position and velocities
        pos = ((f[0:nparticles] % L) + L) % L
        vel = f[nparticles:]
        
        # Send to output functions
        for func in out:
            func(pos, vel, ncells, L, time)
        
    return pos, vel

####################################################################
# 
# Output functions and classes
#

class Plot:
    """
    Displays three plots: phase space, charge density, and velocity distribution
    """
    def __init__(self, pos, vel, ncells, L):
        
        d = calc_density(pos, ncells, L)
        vhist, bins  = histogram(vel, int(sqrt(len(vel))))
        vbins = 0.5*(bins[1:]+bins[:-1])
        
        # Plot initial positions
        if got_gridspec:
            self.fig = plt.figure()
            self.gs = gridspec.GridSpec(4, 4)
            ax = self.fig.add_subplot(self.gs[0:3,0:3])
            self.phase_plot = ax.plot(pos, vel, '.')[0]
            ax.set_title("Phase space")
            
            ax = self.fig.add_subplot(self.gs[3,0:3])
            self.density_plot = ax.plot(linspace(0, L, ncells), d)[0]
            
            ax = self.fig.add_subplot(self.gs[0:3,3])
            self.vel_plot = ax.plot(vhist, vbins)[0]
        else:
            self.fig = plt.figure()
            self.phase_plot = plt.plot(pos, vel, '.')[0]
            
            self.fig = plt.figure()
            self.density_plot = plt.plot(linspace(0, L, ncells), d)[0]
            
            self.fig = plt.figure()
            self.vel_plot = plt.plot(vhist, vbins)[0]
        plt.ion()
        plt.show()
        
    def __call__(self, pos, vel, ncells, L, t):
        d = calc_density(pos, ncells, L)
        vhist, bins  = histogram(vel, int(sqrt(len(vel))))
        vbins = 0.5*(bins[1:]+bins[:-1])
        
        self.phase_plot.set_data(pos, vel) # Update the plot
        self.density_plot.set_data(linspace(0, L, ncells), d)
        self.vel_plot.set_data(vhist, vbins)
        plt.draw()
        

class Summary:
    def __init__(self):
        self.t = []
        self.firstharmonic = []
        
    def __call__(self, pos, vel, ncells, L, t):
        # Calculate the charge density
        d = calc_density(pos, ncells, L)
        
        # Amplitude of the first harmonic
        fh = 2.*abs(fft(d)[1]) / float(ncells)
        
        print "Time:", t, "First:", fh
        
        self.t.append(t)
        self.firstharmonic.append(fh)

####################################################################
# 
# Functions to create the initial conditions
#

def landau(npart, L, alpha=0.2):
    """
    Creates the initial conditions for Landau damping
    
    """
    # Start with a uniform distribution of positions
    pos = random.uniform(0., L, npart)
    pos0 = pos.copy()
    k = 2.*pi / L
    for i in range(10): # Adjust distribution using Newton iterations
        pos -= ( pos + alpha*sin(k*pos)/k - pos0 ) / ( 1. + alpha*cos(k*pos) )
        
    # Normal velocity distribution
    vel = random.normal(0.0, 1.0, npart)
    
    return pos, vel

def twostream(npart, L, vbeam=2):
    # Start with a uniform distribution of positions
    pos = random.uniform(0., L, npart)
    # Normal velocity distribution
    vel = random.normal(0.0, 1.0, npart)
    
    np2 = int(npart / 2)
    vel[:np2] += vbeam  # Half the particles moving one way
    vel[np2:] -= vbeam  # and half the other
    
    return pos,vel

####################################################################

if __name__ == "__main__":

    #vlist = [1., 5., 10.]
        # 2-stream instability
    #for v in vlist:    
    L = 100
    ncells = 20
    pos, vel = twostream(10000, L, 3.)
   
 
    # Create some output classes
    #p = Plot(pos, vel, ncells, L) # This displays an animated figure
    s = Summary()                 # Calculates, stores and prints summary info
    
    # Run the simulation
    pos, vel = run(pos, vel, L, ncells, 
                out=[s],
                #out=[p, s],                      # These are called each output
                output_times=linspace(0.,40,100)) # The times to output
                    
                    
    x = s.t #time
    data = s.firstharmonic #harmonics y values
    
    c = (diff(sign(diff(data))) < 0).nonzero()[0] + 1 # local max
	
    timemax = []
    maxh = []
    rmaxh = []
    
    for i in c:
        timemax.append(x[i])
        maxh.append(data[i])
    
    #print maxh
    
    for i in reversed(maxh):
        print i
        rmaxh.append(i)
    #print rmaxh
    
    def crap(rmaxh):             
	    for j in range(len(rmaxh)-1):  
		if not rmaxh[j] > rmaxh[j+1]:
		    sliceindex = j
                    slicevalue = rmaxh[j]
		    
		    return sliceindex
    
    sliceindex = crap(rmaxh)   
    
    def twostream(timemax, rmaxh):
        timestream = timemax[: sliceindex]
        timenoisy = timemax[sliceindex :]
        rmaxhstream = rmaxh[: sliceindex]
        rmaxhnoisy= rmaxh[sliceindex:]  
        
        print timestream
        
        print timenoisy
  
            
        #coefs = polyfit(timestream, maxhstream,1)
        #straightfit = polyval(polyfunc,timestream)
        #return timestream, straightfit
    #timestream, straightfit, = twostream(timemax, log(rmaxh))
   

    # Summary stores an array of the first-harmonic amplitude
    # Make a semilog plot to see exponential damping
    plt.figure()
    plt.plot(s.t, s.firstharmonic)
    plt.plot(timemax, maxh, 'or')
    plt.xlabel("Time [Normalised]")
    plt.ylabel("First harmonic amplitude [Normalised]")
    #plt.yscale('log')
    #plt.legend()
    plt.ioff() # This so that the windows stay open
    plt.show()

for i in range(len(maxh)-1):
        goodvalues = where(maxh[i] > min(maxh))
        print goodvalues

import numpy

################ generate some data, parabola plus sinusoid
x = numpy.arange(9999.0) / 9998.0
TAU = 2*numpy.pi                  # http://tauday.com/
xi = (8*TAU) * x
y = numpy.sin(xi) + (x-0.3)**2

################ find the positions of the local maxima
local_maxima = []                         # array of indexes
i = 0
while i < (len(y)-2):
    while (i < (len(y)-2)) and (y[i] < y[i+1]):   # increasing
        i += 1
    local_maxima.append(i)
    i += 1
    while (i < (len(y)-2)) and (y[i+1] < y[i]):   # decreasing
        i += 1
    i += 1

################ bracket the data to fit
ymaxes = numpy.array(y[local_maxima])   # values of max
right = numpy.argmax(ymaxes[1:])+1  # left because I know what the data looks like
left = numpy.argmax(-ymaxes[:-1])   # smallest maximum
R = local_maxima[right]
L = local_maxima[left]

best_fit_line = numpy.polyfit(x[L:R], y[L:R], 1)

for xy in zip(x,y):
    print('{} {}'.format(*xy))

print('')

ex = x[L:R:8]
yfit = numpy.polyval(best_fit_line,ex)
for a in zip(ex, yfit):
    print('{} {}'.format(*a))

################ bracket the data to fit
ymaxes = numpy.array(y[local_maxima])   # values of max
right = numpy.argmax(ymaxes[1:])+1  # left because I know what the data looks like
left = numpy.argmax(-ymaxes[:-1])   # smallest maximum
R = local_maxima[right]
L = local_maxima[left]

best_fit_line = numpy.polyfit(x[L:R], y[L:R], 1)

for xy in zip(x,y):
    print('{} {}'.format(*xy))

x = s.t #time
    y = s.firstharmonic #harmonics y values
    
    local_maxima = []                         # array of indexes
    i = 0
    while i < (len(y)-2):
        while (i < (len(y)-2)) and (y[i] < y[i+1]):   # increasing
            i += 1
        local_maxima.append(i)
        i += 1
        while (i < (len(y)-2)) and (y[i+1] < y[i]):   # decreasing
            i += 1
        i += 1
        
    print local_maxima
   
    ymax = array(y[local_maxima])

error: ymax = array(y[local_maxima])
list indices must be integers, not list

ymax = [y[int(j)] for j in local_maxima]

x = s.t #time
    y = s.firstharmonic #harmonics
    maxh = []
    timemax = []
    
    c = (diff(sign(diff(y))) < 0).nonzero()[0] + 1 # local maximum actual values
    for a in c:
	    timemax.append(x[a])
	    maxh.append(y[a])

################ find the positions of the local maxima
    local_maxima = []                         # array of indexes
    i = 0
    while i < (len(y)-2):
        while (i < (len(y)-2)) and (y[i] < y[i+1]):   # increasing
            i += 1
        local_maxima.append(i)
        i += 1
        while (i < (len(y)-2)) and (y[i+1] < y[i]):   # decreasing
            i += 1
        i += 1

################ bracket the data to fit

    ymax = array([y[int(j)] for j in local_maxima]) #Turning the list of ymax values into an array
    
    right = argmax(ymax[1:])+1
    left = argmax(-ymax[:-1])
    R = local_maxima[right]
    L = local_maxima[left]
    
    best_fit_line = polyfit(x[L:R], y[L:R], 1)
    yfit = polyval(best_fit_line, x)
    # Summary stores an array of the first-harmonic amplitude
    # Make a semilog plot to see exponential damping
    plt.figure()
    plt.plot(s.t, s.firstharmonic)
    plt.plot(timemax, maxh, 'or')
    plt.plot(x[L:R], yfit[L:R])
    plt.xlabel("Time [Normalised]")
    plt.ylabel("First harmonic amplitude [Normalised]")
    #plt.yscale('log')
    #plt.legend()
    plt.ioff() # This so that the windows stay open
    plt.show()