#A Sage script  to fit the Hill-Keller model v'(t) = P - k*v(t) to the data from
#Usain Bolt's 2008 Beijing gold medal Olympic race.  This notebook fits both k and P
#to fit the data, in a least-squares sense.

#The Data: Here is the data for Usain Bolt's 2008 Beijing race. The data consists
#of pairs [t, d], where d is 0, 10, 20, ..., 100 meters and t is the corresponding
#split time.

bolttimes = [0.165, 1.85, 2.87, 3.78, 4.65, 5.50, 6.32, 7.14, 7.96, 8.79, 9.69];
dists = [10*i for i in range(11)];
data = list(zip(bolttimes,dists));

#A plot
plt1 = scatter_plot(data,rgbcolor=[1,0,0]);
show(plt1)

#Solving The ODE: Solve the Hill-Keller ODE, using P = 11 and treating "k" as
#unknown. Also use initial condition v(0.165) = 0.
var('k P t');
v = function('v')(t)
t0 = 0.165;
de = diff(v,t) == P - k*v
vsol(t) = desolve(de,v,[t0,0],t);

#Integrate to obtain the position in terms of t, k, and P:
var('tau');
X(t,k,P) = integral(vsol(tau),tau,t0,t);

#Let's guess a value of k and plot X(t) with the data.
plt2 = plot(X(t,1,11), t, 0, 9.69);
pp = plt1 + plt2;
show(plt1+plt2)

#The Optimal Choice for k and P: Not bad, but we can do better by forming a sum of squares SS and minimizing with respect to k and P.
SS = function('SS')(k,P);
SS(k,P) = add((X(bolttimes[i])-dists[i])^2 for i in range(11));

#A quick plot of log(SS(k,P))
plot3d(log(SS), (k,0.7,1.1), (P,10,11)).show()

#Somewhere around k = 0.9, P=10.5 it appears there is a minimum.
#We could minimize using Calculus 3 techniques, but Sage has a
#built-in command "find_fit" to find the least-squares model fit to the data:
model(t) = X(t,k,P); #Specify the model to be fit, with "t" as the independent variable
sol = find_fit(data,model,parameters=[k,P],initial_guess=[0.9,10.5]) #Fit the model by adjusting k and b
kbest = sol[0].rhs();
Pbest = sol[1].rhs();
bestX(t) = model(k=kbest,P=Pbest) #Define the best fit function of the form in "model"

#Plot the function X(t) with this best choice for k and P.
plt2 = plot(bestX(t), t, 0, 9.69);
pp = plt1 + plt2;
show(plt1+plt2)