def tryFits1(fName):
    distances, heights = getTrajectoryData(fName)
    distances = pylab.array(distances)*36
    totHeights = pylab.array([0]*len(distances))
    for h in heights:
        totHeights = totHeights + pylab.array(h)
    pylab.title('Trajectory of Projectile (Mean of 4 Trials)')
    pylab.xlabel('Inches from Launch Point')
    pylab.ylabel('Inches Above Launch Point')
    meanHeights = totHeights/float(len(heights))
    pylab.plot(distances, meanHeights, 'bo')
    a,b = pylab.polyfit(distances, meanHeights, 1)
    altitudes = a*distances + b
    pylab.plot(distances, altitudes, 'r',
               label = 'Linear Fit' + ', R2 = '
               + str(round(rSquare(meanHeights, altitudes), 4)))
    a,b,c = pylab.polyfit(distances, meanHeights, 2)
    altitudes = a*(distances**2) + b*distances + c
    pylab.plot(distances, altitudes, 'g',
               label = 'Quadratic Fit' + ', R2 = '
               + str(round(rSquare(meanHeights, altitudes), 4)))
    pylab.legend()

##tryFits1('launcherData.txt')
##pylab.show()

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 processTrajectories(filename):
    dist, heights = getTrajectory(filename)
    trials = len(heights)
    dist = pylab.array(dist)

    # get array combiming mean height at each distances
    totHeights = pylab.array([0] * len(dist))
    for h in heights:
        totHeights = totHeights + pylab.array(h)
    meanHeights = totHeights/len(heights)

    # start plotting :)
    pylab.title("Trajectile of Projectile " + "(Mean of " + str(trials)
                + " trails)")
    pylab.xlabel("Inches from launch point")
    pylab.ylabel("Inches above launch point")
    pylab.plot(dist, meanHeights, 'bo')

    a, b = pylab.polyfit(dist, meanHeights, 1)
    altitudes = a * dist + b
    pylab.plot(dist, altitudes, 'b', label='linear fit')
    print 'RSquared of linear fit: ', rSquared(meanHeights, altitudes)

    a, b, c = pylab.polyfit(dist, meanHeights, 2)
    altitudes = a*(dist**2) + b*(dist) + c
    pylab.plot(dist, altitudes, 'b:', label="Quadratic Fit")
    print 'RSquared of quadratic fit: ', rSquared(meanHeights, altitudes)
    pylab.legend()

    pylab.show()

def tryFits1(fName):
    distances, heights = getTrajectoryData(fName)
    distances = pylab.array(distances)*36
    totHeights = pylab.array([0]*len(distances))
    for h in heights:
        totHeights = totHeights + pylab.array(h)
    pylab.title('Trajectory of Projectile (Mean of 4 Trials)')
    pylab.xlabel('Inches from Launch Point')
    pylab.ylabel('Inches Above Launch Point')
    meanHeights = totHeights/float(len(heights))
    pylab.plot(distances, meanHeights, 'bo')
    a,b = pylab.polyfit(distances, meanHeights, 1)
    altitudes = a*distances + b
    pylab.plot(distances, altitudes, 'r',
               label = 'Linear Fit' + ', R2 = '
               + str(round(rSquare(meanHeights, altitudes), 4)))
    a,b,c = pylab.polyfit(distances, meanHeights, 2)
    altitudes = a*(distances**2) + b*distances + c
    pylab.plot(distances, altitudes, 'g',
               label = 'Quadratic Fit' + ', R2 = '
               + str(round(rSquare(meanHeights, altitudes), 4)))
    pylab.legend()

##tryFits1('launcherData.txt')
##pylab.show()

# Calculating the location of the peak of height:
# −(b/2a)=547.1

# Once we have the curve fitting value, we can derrive from it other measurements like speed, strenght etc, using known physics.

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 tryFits(fName):
    '''
    Linear fit does not fit the data. Not a logical assumption that the arrow
    flies in a straight line to the target.

    Quadratic fit mirrors a parabolic pathway.
    '''
    distances, heights = getTrajectoryData(fName)
    distances = pylab.array(distances)*36 # Convert yards to feet
    totHeights = pylab.array([0]*len(distances))
    for h in heights:
        totHeights = totHeights + pylab.array(h) # Get one avg measurement of height
    pylab.title('Trajectory of Projectile (Mean of 4 Trials)')
    pylab.xlabel('Inches from Launch Point')
    pylab.ylabel('Inches Above Launch Point')
    meanHeights = totHeights/float(len(heights))
    pylab.plot(distances, meanHeights, 'bo')
    a,b = pylab.polyfit(distances, meanHeights, 1)
    altitudes = a*distances + b
    pylab.plot(distances, altitudes, 'r',
               label = 'Linear Fit')
    a,b,c = pylab.polyfit(distances, meanHeights, 2)
    altitudes = a*(distances**2) + b*distances + c 
    pylab.plot(distances, altitudes, 'g',
               label = 'Quadratic Fit')
    pylab.legend()

def visualise(fileName, title="", linearFit=False, polyFit=False):
    """ Draw graph representing the result. A bit verbose as that seem to be the only way to make borders gray.
    fileName = path and name of the pickled results
    """
    with open(fileName, "rb") as f:
        data = pickle.load(f)
        fig = plt.figure()
        p = fig.add_subplot(111)
        if not linearFit and not polyfit:
            p.plot(data[0], data[1], 'bo-', label="sentiment")
            p.plot([data[0][0], data[0][-1]], [data[1][0], data[1][-1]],
                   'g', label="straight line through first and last point")
        elif linearFit:
            fit = polyfit(data[0], data[1], 1)
            fitFunc = poly1d(fit)
            p.plot(data[0], data[1], '-ro', label='sentiment')
            p.plot(data[0], fitFunc(data[0]), "--k", label="linear fit")
        elif polyFit:
            fit = polyfit(data[0], data[1], 2)
            f = [d*d*fit[0] + d*fit[1] + fit[2] for d in data[0]]
            p.plot(data[0], data[1], '-ro', label='sentiment')
            p.plot(data[0], f, "--k", label="polynomial fit")
        p.legend(prop={'size': 10}, frameon=False)
        plt.ylabel("Average happiness")
        plt.xlabel("Rating")
        for e in ['bottom', 'top', 'left', 'right']:
            p.spines[e].set_color('gray')
        if title:
            plt.title(title)
        plt.show()

def tryFits1(fName):
    '''
    Calculate the coefficient of determination (R**2) to determine how
    well the model fits the data and ergo could make predictions.
    '''
    distances, heights = getTrajectoryData(fName)
    distances = pylab.array(distances)*36
    totHeights = pylab.array([0]*len(distances))
    for h in heights:
        totHeights = totHeights + pylab.array(h)
    pylab.title('Trajectory of Projectile (Mean of 4 Trials)')
    pylab.xlabel('Inches from Launch Point')
    pylab.ylabel('Inches Above Launch Point')
    meanHeights = totHeights/float(len(heights))
    pylab.plot(distances, meanHeights, 'bo')
    a,b = pylab.polyfit(distances, meanHeights, 1)
    altitudes = a*distances + b
    pylab.plot(distances, altitudes, 'r',
               label = 'Linear Fit' + ', R2 = '
               + str(round(rSquare(meanHeights, altitudes), 4)))
    a,b,c = pylab.polyfit(distances, meanHeights, 2)
    altitudes = a*(distances**2) + b*distances + c
    pylab.plot(distances, altitudes, 'g',
               label = 'Quadratic Fit' + ', R2 = '
               + str(round(rSquare(meanHeights, altitudes), 4)))
    pylab.legend()

def plot_dco_values(ax, values, color="k"):
    interpol1 = {"temperature": values["temperature"][0:7], "values": values["values"][0:7]}
    fit1 = pylab.polyfit(interpol1["temperature"], interpol1["values"], 1)
    print "m={} b={}".format(fit1[0], fit1[1])
    fit_fn1 = pylab.poly1d(fit1)

    interpol2 = {"temperature": values["temperature"][6:14], "values": values["values"][6:14]}
    fit2 = pylab.polyfit(interpol2["temperature"], interpol2["values"], 1)
    print "m={} b={}".format(fit2[0], fit2[1])
    fit_fn2 = pylab.poly1d(fit2)

    plot = ax.plot(
        interpol1["temperature"],
        fit_fn1(interpol1["temperature"]),
        "k-",
        interpol2["temperature"],
        fit_fn2(interpol2["temperature"]),
        "k-",
        # values['temperature'], values['values'], '{}-'.format(color),
        values["temperature"],
        values["values"],
        "{}o".format(color),
        markersize=5,
    )
    pylab.setp(plot[0], linewidth=2)
    pylab.setp(plot[1], linewidth=2)

    return plot

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 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 fitData2(fileName):
    '''
    Models predictions using Terman's law model (cubic fit) and 
    Hooks Law (linear fit).

    Hook's Law functions up to the point where the spring reaches 
    it's elastic limit - when it stops behaving as a spring but instead 
    as a rope, etc (doesn't get longer b/c hang more weight on it)
    '''
    xVals, yVals = getData(fileName)
    extX = pylab.array(xVals + [1.05, 1.1, 1.15, 1.2, 1.25])
    xVals = pylab.array(xVals)
    yVals = pylab.array(yVals)
    xVals = xVals*9.81  # convert mass to force (F = mg)
    extX = extX*9.81    # convert mass to force (F = mg)
    pylab.plot(xVals, yVals, 'bo', label = 'Measured displacements')
    pylab.title('Measured Displacement of Spring')
    pylab.xlabel('|Force| (Newtons)')
    pylab.ylabel('Distance (meters)')
    a,b = pylab.polyfit(xVals, yVals, 1)
    estYVals = a*extX + b
    pylab.plot(extX, estYVals, label = 'Linear fit')
    a,b,c,d = pylab.polyfit(xVals, yVals, 3)
    estYVals = a*(extX**3) + b*extX**2 + c*extX + d
    pylab.plot(extX, estYVals, label = 'Cubic fit')
    pylab.legend(loc = 'best')

def findOrder(xVals, yVals, accuracy = 1.0e-1):
    # Your Code Here
    i = 0
    r = pylab.polyfit(xVals, yVals, 0, full=True)[1][0]
    while r > accuracy:
        i += 1
        r = pylab.polyfit(xVals, yVals, i, full=True)[1][0]
    s = pylab.polyfit(xVals, yVals, i, full=True)[0]
    return s

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 fitAverage(self,
                   linear_fit_min_value,
                   linear_fit_max_value,

                   degree,

                   # function taking the list of parameters
                   # and translating it into a label for writing 
                   # a text on the graph
                   # with the fit result
                   label_func,
                   ):
        """ perform a linear fit to the average values """
        import utils
        global linear_fit_min_num_vertices, linear_fit_max_num_vertices

        # filter the values for the fit
        xpos_for_fit = []
        ypos_for_fit = []
        uncert_for_fit = []

        oneSigmaIndex = 2
        assert standardQuantileHistoDefs[oneSigmaIndex]['title'] == '1 sigma'

        for x,y, sigmaUp, sigmaDown in zip(self.xpos, self.avg_values, self.quantile_values_upper[oneSigmaIndex], self.quantile_values_lower[oneSigmaIndex]):
            if x >= linear_fit_min_value and \
               x <= linear_fit_max_value:

                xpos_for_fit.append(x)

                # mean
                ypos_for_fit.append(y)
                
                # 1 sigma 
                uncert_for_fit.append(0.5 * (sigmaUp - sigmaDown))

        if len(xpos_for_fit) < 2:
            raise Exception("must have at least two points for a linear fit (linear_fit_min_value=" + str(linear_fit_min_value) + " linear_fit_max_value=" + str(linear_fit_max_value) + ")")

        # perform fit to means
        # note that pylab.polyfit returns coefficients in an order
        # where the first coefficient corresponds to the highest power of x, 
        # we reverse this
        import pylab

        fittedCoeffs = pylab.polyfit(xpos_for_fit, ypos_for_fit, degree)
        fittedCoeffs = fittedCoeffs[::-1]
        self.addFitResult('mean', linear_fit_min_value, linear_fit_max_value, self.meanFitResult, fittedCoeffs)

        self.fitResultLabel = label_func(fittedCoeffs * self.y_scale_factor, self.yaxis_unit_label)
        
        # perform fit to symmetrized 1 sigma bands
        fittedCoeffs = pylab.polyfit(xpos_for_fit, uncert_for_fit, degree)
        fittedCoeffs = fittedCoeffs[::-1]
        self.addFitResult('uncert', linear_fit_min_value, linear_fit_max_value, self.uncertFitResult, fittedCoeffs)

def main(filename):
	pOut = parse( filename)
	oM, oB = polyfit(pOut[0], pOut[1], 1)
	tM, tB = polyfit(pOut[2], pOut[3], 1)
	print("rep_overlap = {} * overlap + {}".format(oM, oB))
	print("rep_tanimoto = {} * tanimoto + {}".format(tM, tB))
	oC = polyfit(pOut[0], pOut[1], 0)
	tC = polyfit(pOut[2], pOut[3], 0)
	print("rep_overlap = {} * overlap".format(oM, oB))
	print("rep_tanimoto = {} * tanimoto + {}".format(tM, tB))
	makeGraphs( pOut)

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 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 anscombe(plotPoints):
    dataFile = open('anscombe.txt', 'r')
    X1,X2,X3,X4,Y1,Y2,Y3,Y4 = [],[],[],[],[],[],[],[]
    for line in dataFile:
        x1,y1,x2,y2,x3,y3,x4,y4 = line.split()
        X1.append(float(x1))
        X2.append(float(x2))
        X3.append(float(x3))
        X4.append(float(x4))
        Y1.append(float(y1))
        Y2.append(float(y2))
        Y3.append(float(y3))
        Y4.append(float(y4))
    dataFile.close()
    xVals = pylab.array(range(21))
    if plotPoints: pylab.plot(X1, Y1, 'o')
    a, b = pylab.polyfit(X1, Y1, 1)
    yVals = a*xVals + b
    pylab.plot(xVals, yVals)
    pylab.xlim(0, 20)
    mean = sum(Y1)/float(len(Y1))
    median = Y1[len(Y1)/2 + 1]
    pylab.title('Mean = ' + str(mean) + ', Median = ' + str(median))
    pylab.figure()
    if plotPoints: pylab.plot(X2, Y2, 'o')
    a, b = pylab.polyfit(X2, Y2, 1)
    yVals = a*xVals + b
    pylab.plot(xVals, yVals)
    pylab.xlim(0, 20)
    mean = sum(Y1)/float(len(Y1))
    median = Y1[len(Y1)/2 + 1]
    pylab.title('Mean = ' + str(mean) + ', Median = ' + str(median))
    pylab.figure()
    if plotPoints: pylab.plot(X3, Y3, 'o')
    a, b = pylab.polyfit(X3, Y3, 1)
    yVals = a*xVals + b
    pylab.plot(xVals, yVals)
    pylab.xlim(0, 20)
    mean = sum(Y1)/float(len(Y1))
    median = Y1[len(Y1)/2 + 1]
    pylab.title('Mean = ' + str(mean) + ', Median = ' + str(median))
    pylab.figure()
    if plotPoints: pylab.plot(X4, Y4, 'o')
    a, b = pylab.polyfit(X4, Y4, 1)
    yVals = a*xVals + b
    pylab.plot(xVals, yVals)
    pylab.xlim(0, 20)
    mean = sum(Y1)/float(len(Y1))
    median = Y1[len(Y1)/2 + 1]
    pylab.title('Mean = ' + str(mean) + ', Median = ' + str(median))
    pylab.show()