Program Listing for File augmento2.m
↰ Return to documentation for file (matlab/@amimodel/augmento2.m
)
function [modelo2] = augmento2(this)
% augmento2 augments the system equation to also include equations for
% sensitivity equation. This will enable us to compute second order
% sensitivities in a forward-adjoint or forward-forward apporach later on.
%
% Parameters:
%
% Return values:
% this: augmented system which contains symbolic definition of the
% original system and its sensitivities @type amimodel
syms Sx Sdot Sy S0
augmodel.nxtrue = length(this.sym.x); % number of states
augmodel.nytrue = length(this.sym.y); % number of observables
if(this.nevent>0)
augmodel.nztrue = length([this.event.z]); % number of observables
else
augmodel.nztrue = 0;
end
np = this.np;
% augment states
Sx = sym(zeros(length(this.sym.x),np));
for j = 1:length(this.sym.x)
for k = 1:np
eval(['syms S' num2str(j) '_' num2str(k)]);
eval(['Sx(j,k) = S' num2str(j) '_' num2str(k) ';']);
end
end
Sdot = jacobian(this.sym.xdot,this.sym.x)*Sx+jacobian(this.sym.xdot,this.sym.p);
% augment output
Sy = jacobian(this.sym.y,this.sym.x)*Sx+jacobian(this.sym.y,this.sym.p);
% generate deltasx
this.getFun([],'deltasx');
this.getFun([],'sz');
this.getFun([],'srz');
for ievent = 1:this.nevent;
if(numel(this.event(ievent).z)>0)
Sz = this.fun.sz.sym(this.z2event==ievent,:);
for ip = 1:np
Sz(:,ip) = subs(Sz(:,ip),this.fun.sx.sym(:,ip),Sx(:,ip));
end
Srz = this.fun.srz.sym(this.z2event==ievent,:);
for ip = 1:np
Srz(:,ip) = subs(Srz(:,ip),this.fun.sx.sym(:,ip),Sx(:,ip));
end
znew = [this.event(ievent).z,reshape(Sz,[sum(this.z2event==ievent),np])];
else
znew = this.event(ievent).z;
end
tmp=subs(this.fun.deltasx.sym(:,:,ievent),this.fun.xdot.strsym_old,this.fun.xdot.sym);
tmp=subs(tmp,this.fun.xdot.strsym,subs(this.fun.xdot.sym,this.fun.x.sym,this.fun.x.sym+this.event(ievent).bolus));
tmp=subs(tmp,this.fun.stau.strsym,this.fun.stau.sym(ievent,:));
for ip = 1:np
tmp(:,ip)=subs(tmp(:,ip),this.fun.sx.sym(:,ip), Sx(:,ip));
end
bolusnew = [this.event(ievent).bolus;tmp(:)];
% replace sx by augmented x
for ip = 1:np
bolusnew(this.nxtrue*ip+(1:this.nxtrue)) = mysubs(bolusnew(this.nxtrue*ip+(1:this.nxtrue)), this.fun.sx.sym(:,ip),Sx(:,ip));
end
augmodel.event(ievent) = amievent(this.event(ievent).trigger,bolusnew,znew);
end
% augment likelihood
this.getFun([],'dsigma_ydp');
this.getFun([],'y');
this.getFun([],'dydp');
this.getFun([],'Jy');
tmp = arrayfun(@(x) sym(['var_y_' num2str(x)]),0:(augmodel.nytrue*(1+np)-1),'UniformOutput',false);
tmp = transpose([tmp{:}]);
aug_y_strsym = reshape(tmp((augmodel.nytrue+1):end),[augmodel.nytrue,np]);% the update Jy must not contain any
% references to x (otherwise we have to deal with sx in sJy), thus
% we replace sy with the corresponding augmented y that we are about to
% create
SJy = jacobian(this.fun.Jy.sym,this.sym.p) ...
+ jacobian(this.fun.Jy.sym,this.fun.sigma_y.strsym)*this.fun.dsigma_ydp.sym ...
+ jacobian(this.fun.Jy.sym,this.fun.y.strsym)*aug_y_strsym;
this.getFun([],'dsigma_zdp');
this.getFun([],'rz');
this.getFun([],'Jz');
tmp = arrayfun(@(x) sym(['var_z_' num2str(x)]),0:(augmodel.nztrue*(1+np)-1),'UniformOutput',false);
tmp = transpose([tmp{:}]);
aug_z_strsym = reshape(tmp((augmodel.nztrue+1):end),[augmodel.nztrue,np]);% the update Jz must not contain any
% references to x (otherwise we have to deal with sx in sJy), thus
% we replace sy with the corresponding augmented y that we are about to
% create
SJz = jacobian(this.fun.Jz.sym,this.sym.p);
if(~isempty(this.fun.sigma_z.strsym))
SJz = SJz + jacobian(this.fun.Jz.sym,this.fun.sigma_z.strsym)*this.fun.dsigma_zdp.sym ...
+ jacobian(this.fun.Jz.sym,this.fun.z.strsym)*aug_z_strsym;
end
this.getFun([],'Jrz');
tmp = arrayfun(@(x) sym(['var_rz_' num2str(x)]),0:(augmodel.nztrue*(1+np)-1),'UniformOutput',false);
tmp = transpose([tmp{:}]);
aug_rz_strsym = reshape(tmp((augmodel.nztrue+1):end),[augmodel.nztrue,np]);% the update Jz must not contain any
% references to x (otherwise we have to deal with sx in sJy), thus
% we replace sy with the corresponding augmented y that we are about to
% create
SJrz = jacobian(this.fun.Jrz.sym,this.sym.p);
if(~isempty(this.fun.sigma_z.strsym))
SJrz = SJrz + jacobian(this.fun.Jrz.sym,this.fun.sigma_z.strsym)*this.fun.dsigma_zdp.sym ...
+ jacobian(this.fun.Jrz.sym,this.fun.rz.strsym)*aug_rz_strsym;
end
% augment sigmas
this.getFun([],'sigma_y');
this.getFun([],'sigma_z');
S0 = jacobian(this.sym.x0,this.sym.p);
augmodel.sym.x = [this.sym.x;Sx(:)];
augmodel.sym.xdot = [this.sym.xdot;Sdot(:)];
augmodel.sym.f = augmodel.sym.xdot;
augmodel.sym.y = [this.sym.y;Sy(:)];
augmodel.sym.x0 = [this.sym.x0;S0(:)];
augmodel.sym.Jy = [this.sym.Jy,SJy];
augmodel.sym.Jz = [this.sym.Jz,SJz];
augmodel.sym.Jrz = [this.sym.Jrz,SJrz];
if(numel([this.event.z]))
augmodel.sym.rz = [this.fun.rz.sym,Srz];
end
augmodel.sym.p = this.sym.p;
augmodel.sym.k = this.sym.k;
augmodel.sym.sigma_y = [transpose(this.sym.sigma_y(:)), reshape(transpose(this.fun.dsigma_ydp.sym), [1,numel(this.fun.dsigma_ydp.sym)])];
augmodel.sym.sigma_z = [transpose(this.sym.sigma_z(:)), reshape(transpose(this.fun.dsigma_zdp.sym), [1,numel(this.fun.dsigma_zdp.sym)])];
modelo2 = amimodel(augmodel,[this.modelname '_o2']);
modelo2.o2flag = 1;
modelo2.debug = this.debug;
modelo2.forward = this.forward;
modelo2.adjoint = this.adjoint;
modelo2.param = this.param;
end
function out = mysubs(in, old, new)
% mysubs is a wrapper for ther subs matlab function
%
% Parameters:
% in: symbolic expression in which to replace @type symbolic
% old: expression to be replaced @type symbolic
% new: replacement expression @type symbolic
%
% Return values:
% out: symbolic expression with replacement @type symbolic
if(~isnumeric(in) && ~isempty(old) && ~isempty(symvar(in)))
matVer = ver('MATLAB');
if(str2double(matVer.Version)>=8.1)
out = subs(in, old(:), new(:));
else
out = subs(in, old(:), new(:), 0);
end
else
out = in;
end
end