#A worksheet to start the analysis of the two-compartment pharmacokinetic
#LSD model in Section 6.5.1.

#The Data: Sampling times, in hours
data = [(1/12, 11.1), (1/4, 7.4), (1/2, 6.3), (1, 6.9), (2, 5), (4, 3.1),
(8, 0.80), (1/12, 10.6), (1/4, 7.6), (1/2, 7), (1, 4.8), (2, 2.8),
(4, 2.5), (8, 2), (1/12, 8.7), (1/4, 6.7), (1/2, 5.9), (1, 4.3),
(2, 4.4), (8, 0.3), (1/12, 10.9), (1/4, 8.2), (1/2, 7.9), (1, 6.6),
(2, 5.3), (4, 3.8), (8, 1.2), (1/12, 6.4), (1/4, 6.3), (1/2, 5.1),
(1, 4.3), (2, 3.4), (4, 1.9), (8, 0.7)];

#A plot of the data:
plt1 = line([data[i] for i in range(7)],color='red',legend_label='Subject 1')
plt2 = line([data[i] for i in range(7,14)],color='blue',legend_label='Subject 2')
plt3 = line([data[i] for i in range(14,20)],color='black',legend_label='Subject 3')
plt4 = line([data[i] for i in range(20,27)],color='green',legend_label='Subject 4')
plt5 = line([data[i] for i in range(27,34)],color='purple',legend_label='Subject 5')
p = plt1+plt2+plt3+plt4+plt5;
p.set_legend_options(loc='upper center');
p.axes_labels(['Time (hours)','Concentration (mg/ml)'])
show(p)

#The ODE Model: We use cP(t) for the plasma concentration and cT(t) as the tissue
#concentration.  #We have cP(0) = 12.27 for each subject and cT(0) = 0.
#The model for cP(t) and cT(t) of Section 6.5.1 (equation (6.68)) is
var('t ka kb ke');
cP = function('cP')(t);
cT = function('cT')(t);
deP = diff(cP,t) ==  (-kb - ke)*cP + 0.706*ka*cT;
deT = diff(cT,t) == 1.407*kb*cP - ka*cT;
#where ka, kb, and ke are the constants to be determined from the concentration
#data.

#Unfortunately, Sage's desolve_system() will not solve this system, nor
#will the Laplace transform commands work (the symbolic constants ka, kb,
#and ke in the ODEs cause problems.)

#Instead, let's use eigenvalue/vector techniques to solve this
#ODE system. The system can be expressed as a homogeneous linear
#system of the form x' = A*x where
A = matrix(SR,[[-(kb+ke),0.706*ka],[1.407*kb,-ka]]);

#Eigenvalues/vectors are
eigs = A.eigenvectors_right() #Compute right eigenvectors, A*v = lambda*v

#The eigenvalues are
lambda1 = eigs[0][0];
lambda2 = eigs[1][0];

#The eigenvectors are returned in eigs[1][0] and eigs[1][1] and can be
#defined as vectors as
v1 = vector(eigs[0][1][0]);
v2 = vector(eigs[1][1][0]);

#Define a matrix P with the eigenvectors as the columns:
P = transpose(matrix([v1, v2]))
#as was done in equation (6.18) in the text.

#Let
x0 = vector(SR,[12.27,0]); #Symbolic vector with real entries
#denote the initial condition vector.
#As per the procedure in Section 6.2.2, we can construct the solution to this ODE
#system by first solving Pc = x0 for c, which in Sage takes the form
c = P\x0;

#From (6.17) in the text the solution is then given by "xsol" where
var('t');
xsol = c[0]*exp(lambda1*t)*v1 + c[1]*exp(lambda2*t)*v2;

#Assign the solution components to cPsol(t) and cTsol(t):
cPsol(t,ka,kb,ke) = xsol[0];
cTsol(t,ka,kb,ke) = xsol[1];

#We could for a sum of squarea and minimize it using Calculus 3 techniques, but in Sage it is easier
#to use Sage's built-in command "find_fit" to find the least-squares model
#fit to the data. See the worksheet "paramest_demo.sage".