def plotPopulations(numSteps):
    """ Plots populations of Foxes & Rabbits for given timesteps. """
    rab_pop, fox_pop = runSimulation(numSteps)
    # for i in range(len(rab_pop)):
    #     print(rab_pop[i], fox_pop[i])

    r_style = 'bo'    # blue - continuous line
    f_style = 'ro'    # red  - continuous line

    pylab.figure('Fox / Rabit Populations')

    pylab.plot(rab_pop, r_style, label='Rabbit Pop')
    pylab.plot(fox_pop, f_style, label='Fox Pop')

    pylab.title('Fox / Rabit Populations: {} timesteps'.format(numSteps))
    pylab.xlabel('Time Steps')
    pylab.ylabel('Population')
    pylab.legend(loc='best')

    degree = 2
    rab_coeff = pylab.polyfit(range(len(rab_pop)), rab_pop, degree)
    pylab.plot(pylab.polyval(rab_coeff, range(len(rab_pop))), 'b-')
    fox_coeff = pylab.polyfit(range(len(fox_pop)), fox_pop, degree)
    pylab.plot(pylab.polyval(fox_coeff, range(len(fox_pop))), 'r-')

    pylab.show()

def getApices(y):
    """ 
    returns the time (in frames) and position of initial and final apex
    height from a given trajectory y, which are obtained by fitting a cubic
    spline 
    
    ==========
    parameter:
    ==========
    y : *array* (1D)
        the trajectory. Should ideally start ~1 frame before an apex and end ~1
        frame behind an apex

    ========
    returns:
    ========

    [x0, xF], [y0, yF] : location (in frames) and value of the first and final
    apices

    """
    # the math behind here is: fitting a 2nd order polynomial and finding 
    # the root of its derivative. Here, only the results are applied, that's
    # why it appears like "magic" numbers
    c = dot(array([[.5, -1, .5], [-1.5, 2., -.5], [1., 0., 0.]]), y[:3])
    x0 = -1. * c[1] / (2. * c[0])
    y0 = polyval(c, x0)

    c = dot(array([[.5, -1, .5], [-1.5, 2., -.5], [1., 0., 0.]]), y[-3:])
    xF = -1. * c[1] / (2. * c[0])
    yF = polyval(c, xF)
    xF += len(y) - 3

    return [x0, xF], [y0, yF]

def runSimulation2(numSteps):
    """
    Runs the simulation for `numSteps` time steps.

    Returns a tuple of two lists: (rabbit_populations, fox_populations)
      where rabbit_populations is a record of the rabbit population at the
      END of each time step, and fox_populations is a record of the fox population
      at the END of each time step.

    Both lists should be `numSteps` items long.
    """

    rabbitPopulationOverTime = []
    foxPopulationOverTime = []
    for step in range(numSteps):
        rabbitGrowth()
        rabbitPopulationOverTime.append(CURRENTRABBITPOP)
        foxGrowth()
        foxPopulationOverTime.append(CURRENTFOXPOP)
    print "CURRENTRABBITPOP", CURRENTRABBITPOP, rabbitPopulationOverTime
    print "CURRENTFOXPOP", CURRENTFOXPOP, foxPopulationOverTime
    pylab.plot(range(numSteps), rabbitPopulationOverTime, '-g', label='Rabbit population')
    pylab.plot(range(numSteps), foxPopulationOverTime, '-o', label='Fox population')
    rabbit_coeff = pylab.polyfit(range(len(rabbitPopulationOverTime)), rabbitPopulationOverTime, 2)
    pylab.plot(pylab.polyval(rabbit_coeff, range(len(rabbitPopulationOverTime))))
    fox_coeff = pylab.polyfit(range(len(foxPopulationOverTime)), foxPopulationOverTime, 2)
    pylab.plot(pylab.polyval(fox_coeff, range(len(rabbitPopulationOverTime))))
    pylab.title('Fox and rabbit population in the wood')
    xlabel = "Plot for simulation of {} steps".format(numSteps)
    pylab.xlabel(xlabel)
    pylab.ylabel('Current fox and rabbit population')
    pylab.legend(loc='upper right')
    pylab.tight_layout()
    pylab.show()
    pylab.clf()

def plotFittingCurve(rabbits, foxes):
    r_coeff = pylab.polyfit(range(len(rabbits)), rabbits, 2) 
    f_coeff = pylab.polyfit(range(len(foxes)), foxes, 2)
    pylab.plot(pylab.polyval(r_coeff, range(len(rabbits))), label = "Rabbits Curve")
    pylab.plot(pylab.polyval(f_coeff, range(len(foxes))), label = "Foxes Curve")
    pylab.legend()
    pylab.show()

def runSimulation(numSteps):
    """
    Runs the simulation for `numSteps` time steps.

    Returns a tuple of two lists: (rabbit_populations, fox_populations)
      where rabbit_populations is a record of the rabbit population at the 
      END of each time step, and fox_populations is a record of the fox population
      at the END of each time step.

    Both lists should be `numSteps` items long.
    """
    
    rabbitPop = []
    foxPop = []
    for i in range(numSteps):
        rabbitGrowth()
        foxGrowth()
        rabbitPop.append(CURRENTRABBITPOP)
        foxPop.append(CURRENTFOXPOP)
    #return (rabbitPop,foxPop)
    pylab.plot(rabbitPop)
    pylab.plot(foxPop)
    pylab.show()
    rabbitPopulationOverTime = rabbitPop[:]
    coeff = pylab.polyfit(range(len(rabbitPopulationOverTime)), rabbitPopulationOverTime, 2)
    pylab.plot(pylab.polyval(coeff, range(len(rabbitPopulationOverTime))))
    pylab.show()
    rabbitPopulationOverTime = foxPop[:]
    coeff = pylab.polyfit(range(len(rabbitPopulationOverTime)), rabbitPopulationOverTime, 2)
    pylab.plot(pylab.polyval(coeff, range(len(rabbitPopulationOverTime))))
    pylab.show()

def equil():


    eq_m = -0.0001
    rise_rate = 0.0125
    mag = 20
    ts = linspace(0, 400, 500)
    f = F(ts, mag=mag)
#     g = eq_m * f
    g = G(ts, f, rate=eq_m)
    b = rise_rate * ts

    fg = f + g + b
    plot(ts, f, label='F Equilibration')
    plot(ts, g, label='G Consumption')
    plot(ts, fg, label='F-G (Sniff Evo)')

    rf, rollover_F, vi_F = calc_optimal_eqtime(ts, f)
    rfg, rollover_FG, vi_FG = calc_optimal_eqtime(ts, fg)
    m = 2 * eq_m * mag

    axvline(x=rollover_F, ls='--', color='blue')
    axvline(x=rollover_FG, ls='--', color='red')

    idx = list(ts).index(rollover_F)

    b = fg[idx] - m * rollover_F
    evo = polyval((m, b), ts)
    plot(ts, evo, ls='-.', color='blue', label='Static Evo. A')
    # ee = where(evo > mag)[0]
    # to_F = ts[max(ee)]
    # print 'F', rollover_F, b, to_F


    b = vi_FG - m * rollover_FG
    evo = polyval((m, b), ts)
    plot(ts, evo, ls='-.', color='red', label='Static Evo. B')

    print polyval((m, b), 200)
    ee = where(evo > mag)[0]
#     to_FG = ts[max(ee)]
#     print 'FG', rollover_FG, b, to_FG

#     axvline(x=to_FG, ls='-', color='red')
#     axvline(x=to_F, ls='-', color='blue')

    axhline(y=mag, ls='-', color='black')
#    plot([ti], [mag], 'bo')

    legend(loc=0)
#     ylim(2980, 3020)
    ylim(18, 21)
    xlim(0, 20)
    ylabel('Intensity')
    xlabel('t (s)')
#    fig = gcf()
#    fig.text(0.1, 0.01, 'asdfasfasfsadfsdaf')


    show()

def make_joining_whisker(px,py,dist,lthick,lscore,rthick,rscore):
  w = Whisker_Seg()
  tt = linspace(0,1,round(dist))
  w.x = polyval(px,tt).astype(float32)
  w.y = polyval(py,tt).astype(float32) 
  w.thick  = polyval( [rthick-lthick,lthick], tt ).astype(float32) 
  w.scores = polyval( [rscore-lscore,lscore], tt ).astype(float32)  
  return w

def plotFitSimulation(numSteps):
    ans = runSimulation(numSteps)
    rabbitCoeff = pylab.polyfit(range(numSteps), ans[0], 2)
    foxCoeff = pylab.polyfit(range(numSteps), ans[1], 2)
    print rabbitCoeff, foxCoeff
    pylab.plot(pylab.polyval(rabbitCoeff, range(numSteps)), 'r')
    pylab.plot(pylab.polyval(foxCoeff, range(numSteps)), 'g')
    pylab.title("polyfit result")
    pylab.show()

def compute_join_curvature( px, py ):
  from scipy.integrate import quad
  xp  = polyder( px, 1 )
  xpp = polyder( px, 2 )
  yp  = polyder( py, 1 )
  ypp = polyder( py, 2 )
  pn = polyadd( polymul( xp, ypp ), polymul( yp, xpp )) #numerator
  pd = polyadd( polymul( xp, xp ) , polymul( yp, yp ) ) #denominator
  integrand = lambda t:  fabs(polyval( pn, t )/( polyval( pd, t )**(1.5)) ) 
  return quad(integrand, 0, 1) [0]

def plotSimulation():
    rabbits, foxes = runSimulation(200)
    pylab.plot(rabbits, label='Rabbits')
    pylab.plot(foxes, label='Foxes')
    a, b, c = pylab.polyfit(range(len(rabbits)), rabbits, 2)
    pylab.plot(pylab.polyval([a, b, c], range(len(rabbits))), label='Rabbits')
    d, e, f = pylab.polyfit(range(len(foxes)), foxes, 2)
    pylab.plot(pylab.polyval([d, e, f], range(len(foxes))), label='Foxes')
    pylab.grid()
    pylab.legend()
    pylab.show()

def plotLineFit():
    rabbitPops, foxPops = runSimulation(200)
    steps = [n for n in range(1, 201)]
    coeff = pylab.polyfit(range(len(rabbitPops)), rabbitPops, 2)
    pylab.plot(pylab.polyval(coeff, range(len(rabbitPops))))
    coeff = pylab.polyfit(range(len(foxPops)), foxPops, 2)
    pylab.plot(pylab.polyval(coeff, range(len(foxPops))))
    pylab.title('Fox vs. Rabbit')
    pylab.legend(('Rabbit Pop', 'Fox Pop'))
    pylab.xlabel('Step')
    pylab.ylabel('Population')
    pylab.show()

    def createPlots(data):
        x = data[0]
        y = data[1]

        a,b,c = polyFit(x, y, 2)

        x = pylab.arange(201)
        y = a*x**2 + b*x + c

        print pylab.polyval((a,b,c), 200)
        pylab.plot(x, y)
        pylab.show()

def plot(rabbits, foxes):
    N = len(rabbits)
    pylab.plot(range(N), rabbits, 'go', label = "rabbit pop")
    pylab.plot(range(N), foxes, 'ro', label = "foxes pop")
    rab_coeff = pylab.polyfit(range(N), rabbits, 2)
    fox_coeff = pylab.polyfit(range(N), foxes, 2)
    pylab.plot(pylab.polyval(rab_coeff, range(N)), 'g-', label = "rabbit polyfit")
    pylab.plot(pylab.polyval(fox_coeff, range(N)), 'r-', label = "fox polyfit")
    pylab.xlabel("Time")
    pylab.ylabel("Population")
    pylab.title("Dynamics of Rabbit and Fox Population")
    pylab.legend()
    pylab.show()

def ddPsi(a,s):
    '''
    # second derivative of polynomial
    # a = polynomical constants (needs to be in list/array form so len() command works)
    # s = position along polynomial (needs to be in list/array form so len() command works)
    '''
    n = len(a) # number of variables in polynimal
    a1 = a[:-1]*(pl.arange(n-1,0,-1)+1)*pl.arange(n-1,0,-1)
    a2 = a*(pl.arange(n-1,-1,-1)+2)*(pl.arange(n-1,-1,-1)+1)
    
    psi = pl.polyval(a1,s) - pl.polyval(a2,s)
   
    return psi

def plotSim():
    numSim = 400
    rabbit_populations, fox_populations = runSimulation(numSim)
    pylab.figure(1)
    pylab.plot(range(numSim), rabbit_populations, label = "rabit")
    pylab.plot(range(numSim), fox_populations, label = "fox")

#    pylab.figure(2)
    coeff = pylab.polyfit(range(numSim), rabbit_populations, 2)
    pylab.plot(pylab.polyval(coeff, range(numSim)), label = "rabit fitting")
    coeff = pylab.polyfit(range(numSim), fox_populations, 2)
    pylab.plot(pylab.polyval(coeff, range(numSim)), label = "fox fitting") 
    
    pylab.legend()
    pylab.show()

def visualize(numStep):
    foxPop, rabbitPop = runSimulation(numStep)
    print foxPop
    print rabbitPop
    x = range(numStep)
    # coeff = [a, b, c] where ax^2 + bx + c
    coeff = pylab.polyfit(x, rabbitPop, 2)
    rabbitPred = pylab.polyval(coeff, x)
    coeff = pylab.polyfit(x, foxPop, 2)
    foxPred = pylab.polyval(coeff, x)
    pylab.figure()
    pylab.plot(x, foxPop, "r")
    pylab.plot(x, foxPred, "y")
    pylab.plot(x, rabbitPop, "b")
    pylab.plot(x, rabbitPred, "g")
    pylab.show()

	def PolinomialRegression(self,Data,FileOutPath,hasHeader,polinoimalDegree):
		NVars = len(Data[0])
		# Plot column name if both are the same variables (diagonal).
		# For different variables, a scatter plot with both variables and their regression.
		Header = []
		if (hasHeader):
			Header = Data[0]
			Data = Data[1:]
		else:
			for i in range(NVars):
				Header.append("column "+str(i+1))
				Data = Data

		loc = 1
		for i in range(NVars):
			for j in range(NVars):
				P.subplot(NVars, NVars, loc)
				if (i==j):
					P.text(0.2, 0.5, Header[i], size=12)
				else:
					fit = P.polyfit(Data[:,i],Data[:,j],polinoimalDegree)
					yp = P.polyval(fit,np.sort(Data[:,i]))
	
					P.plot(np.sort(Data[:,i]),yp,'g-',Data[:,i],Data[:,j],'b.')
					P.grid(True)
					
				loc = loc + 1
		
		# Store as SVG
		P.savefig(FileOutPath)

def parse():
    # Split the input up
    rows = [x.strip().split(':') for x in fileinput.input()]
    # Automatically turn numbers into the base type
    rows = [[to_float(y) for y in x] for x in rows]

    # Scan once to calculate the overhead
    r = [Record(*(x + [0, 0, 0])) for x in rows]
    bounces = pylab.array([(x.loops, x.rawtime) for x in r if x.test == 'bounce'])
    fit = pylab.polyfit(bounces[:,0], bounces[:,1], 1)

    records = []

    for row in rows:
        # Make a dummy record so we can use the names
        r1 = Record(*(row + [0, 0, 0]))

        bytes = r1.size * r1.loops
        # Calculate the bounce time
        delta = pylab.polyval(fit, [r1.loops])
        time = r1.rawtime - delta
        rate = bytes / time

        records.append(Record(*(row + [time, bytes, rate])))

    return records

def evaluate_models_on_training(x, y, models):
    """
    For each regression model, compute the R-square for this model with the
    standard error over slope of a linear regression line (only if the model is
    linear), and plot the data along with the best fit curve.

    For the plots, you should plot data points (x,y) as blue dots and your best
    fit curve (aka model) as a red solid line. You should also label the axes
    of this figure appropriately and have a title reporting the following
    information:
        degree of your regression model,
        R-square of your model evaluated on the given data points
    Args:
        x: a list of length N, representing the x-coords of N sample points
        y: a list of length N, representing the y-coords of N sample points
        models: a list containing the regression models you want to apply to
            your data. Each model is a numpy array storing the coefficients of
            a polynomial.
    Returns:
        None
    """
    for m in models:
        estYVals = pylab.polyval(m, x)
        R_v_2 = r_squared(y, estYVals)
        pylab.figure()
        pylab.plot(x, y, 'o', label = 'Actual data')
        pylab.plot(x, estYVals, 'r--', label = 'Predictive values')
        print('R-squared = ', rSquared(y, estYVals))
        pylab.legend(loc = 'best')
        pylab.title('evaluate_models_on_training')
        pylab.show()

def get_minmax(data, deg=2):
    """
    returns the interpolated extremum and position (in fractions of frames)

    :args:
        data (iterable): data to be fitted with a <deg> degree polynomial
        deg (int): degree of polynomial to fit

    :returns:
        val, pos: a tuple of floats indicating the position
    """

    x = arange(len(data))
    p = polyfit(x, data, deg)
    d = (arange(len(p))[::-1] * p)[:-1]
    r = roots(d)
    cnt = 0
    idx = None
    for nr, rx in enumerate(r):
        if isreal(r):
            idx = nr
            cnt +=1
            if cnt > 1:
                raise ValueError("Too many real roots found." + 
                        "Reduce order of polynomial!")
    x0 = r[nr].real
    return x0, polyval(p, x0)