% function X = ad_demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ad_demo.m
% Using the tssolve implementation of Automatic Differentiation, generate
% the sensitivity matrices for a model (in this case, the HIV model) wrt
% all parameters. Then slim down the matrix to the parameters that we're
% interested in.
%
% Written by: Keri Rehm, klrehm@ncsu.edu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Setup
% hiv_rhs_ad is the file that contains the ODE system. Make a function handle
% and assign it to variable 'odename' to make life easier if I want to find
% the sensitivities of a different ODE.
odename = @hiv_rhs_ad;
rho = 6; % Dimension of ODE system
n = 41; % Number of time points
t0 = 0; % Initial time
tf = 200; % Final time
y0 = [0.9e+6; 4000; 1; 1; 1; 12]; % Initial conditions
%%% Make our time vector since we're not pulling in data
tvec = t0:tf/(n-1): tf;
% parameter values
lam1 = 1e+4; d1 = 0.01; epsilon = 0; k1 = 8.0e-7;
lam2 = 31.98; d2 = 0.01; f = 0.34; k2 = 1e-4;
delta = 0.7; m1 = 1.0e-5; m2 = 1.0e-5; NT = 100;
c = 13; rho1 = 1; rho2 = 1; lamE = 1;
bE = 0.3; Kb = 100; dE = 0.25; Kd = 500;
deltaE = 0.1;
%%% Tell computer which parameters we're varying in our inverse problem
% qsel is a vector (that will be used as boolean) that indicates
% which parameters we want sensitivities on. I tried to line things up with
% their spaces in 'q' to make it easy to see which ones we will get
% sensitivities out of. Keep qsel(7) = 0 when epsilon=0.
q = [lam1,d1,epsilon,k1,lam2,d2,f,k2,delta,m1,m2,NT,c,rho1,rho2,lamE,bE,Kb,dE,Kd,deltaE];
qsel = [ 1, 1, 0, 0, 0, 1,0, 1, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0];
% qsel = [ 1, 1, 1, 1, 1, 1,0, 1, 1, 1, 1, 1,1, 1, 1, 1, 1, 1, 1, 1, 1];
p=sum(qsel);
tsoptions = []; % Options for tssolve's ode solver. Ignore for now.
%%% Perform the AD calculations. See tssolve.m for structure of Dq and Dx0
[Dq Dx0] = tssolve(odename,y0,q,tvec,1,tsoptions);
%%% Put the sensitivities into the form of the problem
X = squeeze(Dq.regsens(:,1,:));
for i = 2:rho
X = [X; squeeze(Dq.regsens(:,i,:))];
end
%%% Cut out the columns of X corresponding to parameters we're not
%%% interested in.
qsel = repmat(qsel,n*rho,1);
X = reshape(X(qsel ==1),n*rho,p);
%%% For fun (and to check my work), calculate inv(X'*X)
cov = inv(X'*X);
% end