Program Listing for File augmento2vec.m¶
↰ Return to documentation for file (matlab/@amimodel/augmento2vec.m
)
function [modelo2vec] = augmento2vec(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:
% modelo2vec: 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
augmodel.nztrue = this.nz;
augmodel.coptim = this.coptim;
augmodel.debug = this.debug;
% multiplication vector (extension of kappa
vecs = cell([length(this.sym.p),1]);
for ivec = 1:length(this.sym.p)
vecs{ivec} = sprintf('k_%i', length(this.sym.k) + ivec-1);
end
vec = sym(vecs);
if(this.nevent>0)
augmodel.nztrue = length([this.event.z]); % number of observables
else
augmodel.nztrue = 0;
end
np = this.np;
% augment states
sv = sym('sv',[length(this.sym.x),1]);
Sdot = jacobian(this.sym.xdot,this.sym.x)*sv+jacobian(this.sym.xdot,this.sym.p)*vec;
% augment output
Sy = jacobian(this.sym.y,this.sym.x)*sv+jacobian(this.sym.y,this.sym.p)*vec;
% generate deltasx
this.getFun([],'deltasx');
for ievent = 1:this.nevent;
if(numel(this.event(ievent).z)>0)
Sz = jacobian(this.event(ievent).z,this.sym.x)*sv+jacobian(this.event(ievent).z,this.sym.p)*vec;
znew = [this.event(ievent).z,Sz];
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(subs(tmp,this.fun.stau.strsym,this.fun.stau.sym(ievent,:)),this.fun.sx.sym*vec,sv);
bolusnew = [this.event(ievent).bolus;tmp*vec];
% replace sx by augmented x
bolusnew(this.nxtrue+(1:this.nxtrue)) = mysubs(bolusnew(this.nxtrue+(1:this.nxtrue)), this.fun.sx.sym(:,1),sv);
hflagold = this.event(ievent).hflag;
augmodel.event(ievent) = amievent(this.event(ievent).trigger,bolusnew,znew);
augmodel.event(ievent) = augmodel.event(ievent).setHflag([hflagold;zeros([numel(sv),1])]);
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*2-1),'UniformOutput',false);
aug_y_strsym = transpose([tmp{(augmodel.nytrue+1):end}]); % 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) ...
* vec + jacobian(this.fun.Jy.sym,this.fun.y.strsym)*aug_y_strsym;
this.getFun([],'dsigma_zdp');
this.getFun([],'dzdp');
this.getFun([],'Jz');
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)*Sz;
end
% augment sigmas
this.getFun([],'sigma_y');
this.getFun([],'sigma_z');
S0 = jacobian(this.sym.x0,this.sym.p)*vec;
augmodel.sym.x = [this.sym.x;sv];
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*vec];
augmodel.sym.k = [this.sym.k;vec];
augmodel.sym.p = this.sym.p;
augmodel.sym.sigma_y = [this.sym.sigma_y, transpose(this.fun.dsigma_ydp.sym * vec)];
augmodel.sym.sigma_z = [this.sym.sigma_z, transpose(this.fun.dsigma_zdp.sym * vec)];
modelo2vec = amimodel(augmodel,[this.modelname '_o2vec']);
modelo2vec.o2flag = 2;
modelo2vec.debug = this.debug;
modelo2vec.forward = this.forward;
modelo2vec.adjoint = this.adjoint;
modelo2vec.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