#!/usr/bin/python

import os
import numpy as np
import matplotlib
matplotlib.use( 'Agg' )
import matplotlib.pyplot as plt
import scipy.optimize as optimize
from util import print_h_m_s , print_m_s

def lin_func( x , *p ) :
    return p[0] + p[1] * x + p[2] * x**2

def pow_func( x , *p ) :
    return p[0] * x**p[1]

def do_fit( data , func , p0 ) :
    p , c = optimize.curve_fit( func , data[:,0] , data[:,1] , p0 = p0 )
    l , v = np.linalg.eig( c )
    sig = v.dot( np.sqrt( np.diag( l ) ) ).dot( v.T )
    dp = np.sqrt( np.sum( sig.dot( v )**2 , axis=0 ) )
    
    return p , dp

def read_result_file( fname ) :
    running_paces = []
    f = open( fname , 'r' )
    for line in f :
        e = line.split()
        dist_meters = float(e[0])
        pace_minute = int(e[1])
        pace_second = int(e[2])
        race_year = int(e[3])
        include_race = int(e[4])
        race_name = ' '.join( e[5:] )
        t = ( pace_minute + pace_second/60. )
        # running_paces.append( running_pace( d = dist_meters , t = ( pace_minute*60. + pace_second ) , y = race_year , inc = include_race , name = race_name ) )
        if include_race :
            running_paces.append( [ dist_meters , t ] )
    
    rp = np.array( running_paces )
    p , dp = do_fit( rp , lin_func , p0 = [ 1 , 1 , 1 ] )
    pp , pm = p+dp , p-dp
    plt.scatter( rp[:,0] , rp[:,1] )
    plt.xlim( [0,60] )
    plt.ylim( [0,16] )

    x = np.linspace( 0 , 60 , 100 )
    plt.plot( x , lin_func( x , *p ) , 'r' , linewidth=2.5)
    plt.plot( x , lin_func( x , *pp ) , 'r--' )
    plt.plot( x , lin_func( x , *pm ) , 'r--' )

    p , dp = do_fit( rp , pow_func , p0 = [ 0.5 , 0.5 ] )
    pp , pm = p+dp , p-dp
    plt.plot( x , pow_func( x , *p ) , 'b' , linewidth=2.5)
    plt.plot( x , pow_func( x , *pp ) , 'b--' )
    plt.plot( x , pow_func( x , *pm ) , 'b--' )

    plt.savefig('running_pace.png')

if __name__ == '__main__' :
    read_result_file( 'running_paces.txt' )

@@ -339,24 +339,68 @@
 nc , p0 ) :%0A
+    datax = data%5B:,0%5D%0A    datay = data%5B:,1%5D%0A
     p , c = 
@@ -434,25 +434,17 @@
 data
-%5B:,0%5D
+x
  , data
-%5B:,1%5D
+y
  , p
@@ -605,20 +605,1184 @@
     
-return p , d
+errfunc = lambda p, x, y: func(x,*p) - y%0A    %0A    residuals = errfunc( p, datax, datay)%0A    s_res = np.std(residuals)%0A    ps = %5B%5D%0A    datayerrors = None%0A    # 100 random data sets are generated and fitted%0A    for _ in range(100):%0A      if datayerrors is None:%0A          randomDelta = np.random.normal(0., s_res, len(datay))%0A          randomdataY = datay + randomDelta%0A      else:%0A          randomDelta =  np.array( %5B %5C%0A                             np.random.normal(0., derr,1)%5B0%5D %5C%0A                             for derr in datayerrors %5D )%0A          randomdataY = datay + randomDelta%0A      randomfit, randomcov = %5C%0A          optimize.leastsq( errfunc, p0, args=(datax, randomdataY),%5C%0A                            full_output=0)%0A      ps.append( randomfit )%0A%0A    ps = np.array(ps)%0A    mean_pfit = np.mean(ps,0)%0A    Nsigma = 1. # 1sigma gets approximately the same as methods above%0A                # 1sigma corresponds to 68.3%25 confidence interval%0A                # 2sigma corresponds to 95.44%25 confidence interval%0A    err_pfit = Nsigma * np.std(ps,0)%0A%0A    pfit_bootstrap = mean_pfit%0A    perr_bootstrap = err_pfit%0A    %0A    # return p , dp%0A    return pfit_bootstrap , perr_bootstra
 p%0A%0Ad
@@ -2103,24 +2103,24 @@
 in( e%5B5:%5D )%0A
+
         t = 
@@ -2157,170 +2157,8 @@
 . )%0A
-        # running_paces.append( running_pace( d = dist_meters , t = ( pace_minute*60. + pace_second ) , y = race_year , inc = include_race , name = race_name ) )%0A
     

add code for bootstrap confidence intervals