Skip to content

Commit

Permalink
Merge branch 'stable'
Browse files Browse the repository at this point in the history
  • Loading branch information
fabian.froehlich committed Sep 5, 2016
2 parents ecb4c4d + 264a42a commit 9f4aa72
Show file tree
Hide file tree
Showing 50 changed files with 2,529 additions and 1,168 deletions.
32 changes: 6 additions & 26 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,32 +49,6 @@ index.html

AMICI.mlappinstall

examples/example_1/simulate_model_example_1.m

examples/example_2/simulate_model_example_2.m

examples/example_3/simulate_model_example_3.m

examples/example_4/simulate_model_example_4.m

examples/example_5/simulate_model_example_5.m

examples/example_6/simulate_model_example_6.m

examples/geneExpression/simulate_geneExpression.m

examples/example_1/html/*

examples/example_2/html/*

examples/example_3/html/*

examples/example_4/html/*

examples/example_5/html/*

examples/example_6/html/*

*.mat

mtoc/makeExampleDoc.m
Expand Down Expand Up @@ -170,3 +144,9 @@ examples/example_dirac_adjoint/dxdotdp_model_dirac_adjoint.m
examples/example_adjoint/J_model_adjoint.m

examples/example_adjoint/dxdotdp_model_adjoint.m

examples/example_jakstat_adjoint_hvp/simulate_model_jakstat_adjoint_hvp.m

examples/example_dirac_adjoint_hessVecProd/simulate_model_dirac_adjoint_hessVecProd.m

examples/example_adjoint_hessian/simulate_model_adjoint_hessian.m
51 changes: 40 additions & 11 deletions @amifun/gccode.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,15 +22,35 @@
this.sym = subs(this.sym,sym('D([2], am_min)'),sym('D2am_min'));
end

if(model.splineflag)
for nodes = [3,4,5,10]
for ideriv = 1:nodes
this.sym = subs(this.sym,sym(['D([' num2str(ideriv*2+1) '], spline_pos' num2str(nodes) ')']),sym(['D' num2str(ideriv*2+1) 'spline_pos' num2str(nodes)]));
this.sym = subs(this.sym,sym(['D([' num2str(ideriv*2+1) '], spline' num2str(nodes) ')']),sym(['D' num2str(ideriv*2+1) 'spline' num2str(nodes)]));
% If we have spline, we need to parse them to get derivatives
if (model.splineflag)
symstr = char(this.sym);
if (strfind(symstr, 'spline'))
tokens = regexp(symstr, 't\,\s(\w+\.\w+)\,', 'tokens');
nNodes = round(str2double(tokens{1}));
end
if (regexp(symstr, 'D\(\[(\w+|\w+\,\w+)\]\,.am_spline'))
isDSpline = true;
else
isDSpline = false;
end

if (isDSpline)
[~, nCol] = size(this.sym);
for iCol = 1 : nCol
for iNode = 1 : nNodes
if (model.o2flag)
for jNode = 1:nNodes
this.sym(:,iCol) = subs(this.sym(:,iCol),sym(['D([' num2str(iNode*2+2) ', ' num2str(jNode*2+2) '], am_spline_pos)']),sym(['D' num2str(iNode*2+2) 'D' num2str(jNode*2+2) 'am_spline_pos']));
this.sym(:,iCol) = subs(this.sym(:,iCol),sym(['D([' num2str(iNode*2+2) ', ' num2str(jNode*2+2) '], am_spline)']),sym(['D' num2str(iNode*2+2) 'D' num2str(jNode*2+2) 'am_spline']));
end
end
this.sym(:,iCol) = subs(this.sym(:,iCol),sym(['D([' num2str(iNode*2+2) '], am_spline_pos)']),sym(['D' num2str(iNode*2+2) 'am_spline_pos']));
this.sym(:,iCol) = subs(this.sym(:,iCol),sym(['D([' num2str(iNode*2+2) '], am_spline)']),sym(['D' num2str(iNode*2+2) 'am_spline']));
end
end
end
end

end

cstr = ccode(this.sym);
if(~strcmp(cstr(3:4),'t0'))
Expand All @@ -46,8 +66,17 @@
cstr = strrep(cstr,'log','amilog');
% fix derivatives again (we cant do this before as this would yield
% incorrect symbolic expressions
cstr = regexprep(cstr,'D([0-9]*)([\w]*)\(','D$2\($1,');
cstr = regexprep(regexprep(cstr,'D([0-9]*)([\w]*)\(','D$2\($1,'),'DD([0-9]*)([\w]*)\(','DD$2\($1,');
cstr = strrep(strrep(cstr, 'DDam_spline', 'am_DDspline'), 'Dam_spline', 'am_Dspline');

if (model.splineflag)
if (strfind(symstr, 'spline'))
% The floating numbers after 't' must be converted to integers
cstr = regexprep(cstr, '(spline|spline_pos)\(t\,\w+\.\w+\,', ['$1\(t\,', num2str(nNodes), '\,']);
cstr = regexprep(cstr, '(spline|spline_pos)\((\w+)\,t\,\w+\.\w+\,', ['$1\($2\,t\,', num2str(nNodes), '\,']);
cstr = regexprep(cstr, '(spline|spline_pos)\((\w+)\,(\w+)\,t\,\w+\.\w+\,', ['$1\($2\,$3\,t\,', num2str(nNodes), '\,']);
end
end

if(numel(cstr)>1)

Expand Down Expand Up @@ -83,7 +112,7 @@
cstr = regexprep(cstr,'xdotdp([0-9]+)',['xdotdp\[$1 + ip*' num2str(model.nx) '\]']);

if(strcmp(this.cvar,'qBdot'))
cstr = regexprep(cstr,'qBdot\[([0-9]*)\]','qBdot\[ip]');
cstr = regexprep(cstr,'qBdot\[([0-9]*)\]','qBdot\[ip+np*$1]');
elseif(strcmp(this.cvar,'stau'))
cstr = regexprep(cstr,'stau\[([0-9]*)\]','stau\[ip]');
elseif(strcmp(this.cvar,'y'))
Expand Down Expand Up @@ -124,7 +153,7 @@
cstr = regexprep(cstr,'dydx_([0-9]+)','dydx\[$1]');
cstr = regexprep(cstr,'dydp_([0-9]+)',['dydp\[$1+ip*' num2str(model.ny) ']']);
cstr = regexprep(cstr,'my_([0-9]+)','my\[it+nt*$1]');
cstr = regexprep(cstr,'sdy_([0-9]+)','sd_y\[$1\]');
cstr = regexprep(cstr,'sigma_y_([0-9]+)','sigma_y\[$1\]');
cstr = regexprep(cstr,'dsdydp\[([0-9]*)\]','dsigma_ydp\[$1\]');
if(strcmp(this.cvar,'sJy'))
cstr = regexprep(cstr,'sy_([0-9]+)','sy\[$1\]');
Expand All @@ -141,7 +170,7 @@
cstr = regexprep(cstr,'dzdp_([0-9]+)',['dzdp\[$1+ip*' num2str(model.nz) ']']);
cstr = regexprep(cstr,'mz_([0-9]+)','mz\[nroots[ie]+nmaxevent*$1]');
cstr = regexprep(cstr,'sz_([0-9]+)',['sz\[nroots[ie]+nmaxevent*\($1+ip*' num2str(model.nz) '\)\]']);
cstr = regexprep(cstr,'sdz_([0-9]+)','sd_z\[$1\]');
cstr = regexprep(cstr,'sigma_z_([0-9]+)','sigma_z\[$1\]');
cstr = regexprep(cstr,'dsdzdp\[([0-9]*)\]','dsigma_zdp\[$1\]');
cstr = regexprep(cstr,'z_([0-9]+)','z\[nroots[ie]+nmaxevent*$1\]');
cstr = strrep(cstr,'=','+=');
Expand Down
14 changes: 7 additions & 7 deletions @amifun/getArgs.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,21 +120,21 @@
case 's2root'
this.argstr = '(realtype t, int ie, int *nroots, realtype *s2root, N_Vector x, N_Vector *sx, void *user_data)';
case 'Jy'
this.argstr = '(realtype t, int it, realtype *Jy, realtype *y, N_Vector x, realtype *my, realtype *sd_y, void *user_data)';
this.argstr = '(realtype t, int it, realtype *Jy, realtype *y, N_Vector x, realtype *my, realtype *sigma_y, void *user_data)';
case 'dJydx'
this.argstr = '(realtype t, int it, realtype *dJydx, realtype *y, N_Vector x, realtype *dydx, realtype *my, realtype *sd_y, void *user_data)';
this.argstr = '(realtype t, int it, realtype *dJydx, realtype *y, N_Vector x, realtype *dydx, realtype *my, realtype *sigma_y, void *user_data)';
case 'dJydp'
this.argstr = '(realtype t, int it, realtype *dJydp, realtype *y, N_Vector x, realtype *dydp, realtype *my, realtype *sd_y, realtype *dsigma_ydp, void *user_data)';
this.argstr = '(realtype t, int it, realtype *dJydp, realtype *y, N_Vector x, realtype *dydp, realtype *my, realtype *sigma_y, realtype *dsigma_ydp, void *user_data)';
case 'sJy'
this.argstr = '(realtype t, int it, realtype *sJy, realtype *dJydy, realtype *dJydp, realtype *sy, void *user_data)';
case 'Jz'
this.argstr = '(realtype t, int ie, realtype *Jz, realtype *z, N_Vector x, realtype *mz, realtype *sd_z, void *user_data, void *temp_data)';
this.argstr = '(realtype t, int ie, realtype *Jz, realtype *z, N_Vector x, realtype *mz, realtype *sigma_z, void *user_data, void *temp_data)';
case 'dJzdx'
this.argstr = '(realtype t, int ie, realtype *dJzdx, realtype *z, N_Vector x, realtype *dzdx, realtype *mz, realtype *sd_z, void *user_data, void *temp_data)';
this.argstr = '(realtype t, int ie, realtype *dJzdx, realtype *z, N_Vector x, realtype *dzdx, realtype *mz, realtype *sigma_z, void *user_data, void *temp_data)';
case 'dJzdp'
this.argstr = '(realtype t, int ie, realtype *dJzdp, realtype *z, N_Vector x, realtype *dzdp, realtype *mz, realtype *sd_z, realtype *dsigma_zdp, void *user_data, void *temp_data)';
this.argstr = '(realtype t, int ie, realtype *dJzdp, realtype *z, N_Vector x, realtype *dzdp, realtype *mz, realtype *sigma_z, realtype *dsigma_zdp, void *user_data, void *temp_data)';
case 'sJz'
this.argstr = '(realtype t, int ie, realtype *sJz, realtype *z, N_Vector x, realtype *dzdp, realtype *sz, realtype *mz, realtype *sd_z, realtype *dsigma_zdp, void *user_data, void *temp_data)';
this.argstr = '(realtype t, int ie, realtype *sJz, realtype *z, N_Vector x, realtype *dzdp, realtype *sz, realtype *mz, realtype *sigma_z, realtype *dsigma_zdp, void *user_data, void *temp_data)';
case 'w'
this.argstr = ['(realtype t, N_Vector x, N_Vector dx, void *user_data)'];
case 'dwdp'
Expand Down
103 changes: 88 additions & 15 deletions @amifun/getSyms.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -330,7 +330,14 @@
this.strsym(idx) = Js;

case 'JB'
this.sym=-transpose(model.fun.J.sym);
this.sym = sym(zeros(model.nx, model.nx));
% Augmentation needs a transposition in the submatrices
for ix = 1 : model.ng
for jx = 1 : model.ng
this.sym((ix-1)*model.nxtrue+1:ix*model.nxtrue, (jx-1)*model.nxtrue+1:jx*model.nxtrue) = ...
-transpose(model.fun.J.sym((ix-1)*model.nxtrue+1:ix*model.nxtrue, (jx-1)*model.nxtrue+1:jx*model.nxtrue));
end
end

case 'dxdotdp'
if(~isempty(w))
Expand Down Expand Up @@ -408,19 +415,43 @@
% transform to symbolic variable
vs = sym(vs);
% multiply
this.sym = -transpose(model.fun.J.sym)*vs;
this.sym = model.fun.JB.sym*vs;

case 'xBdot'
if(strcmp(model.wtype,'iw'))
syms t
this.sym = diff(transpose(model.fun.M.sym),t)*model.fun.xB.sym + transpose(model.fun.M.sym)*model.fun.dxB.sym - transpose(model.fun.dfdx.sym)*model.fun.xB.sym;
else
this.sym = -transpose(model.fun.J.sym)*model.fun.xB.sym;
% Augmenting the system needs transposition of submatrices
% I'm sure, there is a more intelligent solution to it...
this.sym = sym(zeros(nx,1));
if(model.o2flag)
for ix = 1 : model.ng
for jx = 1 : model.ng
this.sym((ix-1)*model.nxtrue+1 : ix*model.nxtrue) = this.sym((ix-1)*model.nxtrue+1 : ix*model.nxtrue) - ...
transpose(model.fun.J.sym((ix-1)*model.nxtrue+1 : ix*model.nxtrue , (jx-1)*model.nxtrue+1 : jx*model.nxtrue)) ...
* model.fun.xB.sym((jx-1)*model.nxtrue+1 : jx*model.nxtrue);
end
end
else
this.sym = -transpose(model.fun.J.sym) * model.fun.xB.sym;
end
end

case 'qBdot'
this.sym = -transpose(model.fun.xB.sym)*model.fun.dxdotdp.sym;

% If we do second order adjoints, we have to augment
if (model.nxtrue < nx)
this.sym = sym(zeros(model.ng, model.np));
this.sym(1,:) = -transpose(model.fun.xB.sym(1:model.nxtrue)) * model.fun.dxdotdp.sym(1:model.nxtrue, :);
for ig = 2 : model.ng
this.sym(ig,:) = ...
-transpose(model.fun.xB.sym(1:model.nxtrue)) * model.fun.dxdotdp.sym((ig-1)*model.nxtrue+1 : ig*model.nxtrue, :) ...
-transpose(model.fun.xB.sym((ig-1)*model.nxtrue+1 : ig*model.nxtrue)) * model.fun.dxdotdp.sym(1:model.nxtrue, :);
end
else
this.sym = -transpose(model.fun.xB.sym)*model.fun.dxdotdp.sym;
end

case 'dsigma_ydp'
this.sym = jacobian(model.fun.sigma_y.sym,p);
this = makeStrSyms(this);
Expand Down Expand Up @@ -490,7 +521,7 @@
staus = sym(staus);
% multiply
this.strsym = staus;

case 'deltax'
if(nevent>0)
this.sym = [model.event.bolus];
Expand All @@ -513,7 +544,7 @@
for ievent = 1:nevent
this.sym(:,ievent,:) = jacobian(model.fun.deltax.sym(:,ievent),x);
end

case 'ddeltaxdt'
this.sym = diff(model.fun.deltax.sym,'t');

Expand Down Expand Up @@ -547,9 +578,18 @@
end

case 'deltaqB'
this.sym = sym(zeros(np,nevent));
if (model.nxtrue < nx)
ng_tmp = round(nx / model.nxtrue);
this.sym = sym(zeros(np*ng_tmp,nevent));
else
this.sym = sym(zeros(np,nevent));
end

for ievent = 1:nevent
this.sym(:,ievent) = transpose(model.fun.xB.sym)*squeeze(model.fun.ddeltaxdp.sym(:,ievent,:));
this.sym(1:np,ievent) = transpose(model.fun.xB.sym)*squeeze(model.fun.ddeltaxdp.sym(:,ievent,:));
% This is just a very quick fix. Events in adjoint systems
% have to be implemented in a way more rigorous way later
% on... Some day...
end

case 'deltaxB'
Expand Down Expand Up @@ -625,19 +665,52 @@

case 'Jy'
this.sym = model.sym.Jy;
% replace unify symbolic expression
this.sym = mysubs(this.sym,model.sym.y,model.fun.y.strsym);
case 'dJydy'
this.sym = jacobian(model.fun.Jy.sym,model.fun.y.strsym);
this.sym = sym(zeros(model.nytrue, model.ng, model.ny));
for iy = 1 : model.nytrue
this.sym(iy,:,:) = jacobian(model.fun.Jy.sym(iy,:),model.fun.y.strsym);
end
this = makeStrSyms(this);
case 'dJydx'
this.sym = model.fun.dJydy.sym*model.fun.dydx.strsym;
this.sym = sym(zeros(model.nytrue, model.nxtrue, model.ng)); % Maybe nxtrue is sufficient...
dJydy_tmp = sym(zeros(model.ng, model.ny));
for iy = 1 : model.nytrue
dJydy_tmp(:,:) = model.fun.dJydy.sym(iy,:,:);
this.sym(iy,:,:) = transpose(dJydy_tmp * model.fun.dydx.strsym(:,1:model.nxtrue));
% Transposition is necessary to have things sorted
% correctly in gccode.m
end
disp('');
case 'dJydsigma'
this.sym = jacobian(model.fun.Jy.sym,model.fun.sigma_y.strsym);
this.sym = sym(zeros(model.nytrue, model.ng, model.nytrue));
for iy = 1 : model.nytrue
this.sym(iy,:,:) = jacobian(model.fun.Jy.sym(iy,:),model.fun.sigma_y.strsym(1:model.nytrue));
end
case 'dJydp'
this.sym = model.fun.dJydy.sym*model.fun.dydp.strsym + model.fun.dJydsigma.sym*model.fun.dsigma_ydp.strsym;
this.sym = sym(zeros(model.nytrue, model.np, model.ng));
dJydy_tmp = sym(zeros(model.ng, model.ny));
dJydsigma_tmp = sym(zeros(model.ng, model.nytrue));
for iy = 1 : model.nytrue
dJydy_tmp(:,:) = model.fun.dJydy.sym(iy,:,:);
dJydsigma_tmp(:,:) = model.fun.dJydsigma.sym(iy,:,:);
this.sym(iy,:,:) = transpose(dJydy_tmp * model.fun.dydp.strsym ...
+ dJydsigma_tmp * model.fun.dsigma_ydp.strsym(1:model.nytrue,:));
% Transposition is necessary to have things sorted
% correctly in gccode.m
end
this = makeStrSyms(this);
case 'sJy'
this.sym = model.fun.dJydy.strsym*model.fun.sy.strsym + model.fun.dJydp.strsym;

this.sym = sym(zeros(model.nytrue, model.np, model.ng));
dJydy_tmp = sym(zeros(model.ng, model.ny));
for iy = 1 : model.nytrue
dJydy_tmp(:,:) = model.fun.dJydy.strsym(iy,:,:);
this.sym(iy,:,:) = transpose(dJydy_tmp*model.fun.sy.strsym);
% Transposition is necessary to have things sorted
% correctly in gccode.m
end
this.sym = this.sym + model.fun.dJydp.strsym;
case 'Jz'
this.sym = model.sym.Jz(:);
case 'dJzdz'
Expand Down
4 changes: 2 additions & 2 deletions @amifun/printLocalVars.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ function printLocalVars(this,model,fid)
fprintf(fid,'int ix;\n');
fprintf(fid,['memset(xBdot_tmp,0,sizeof(realtype)*' num2str(nx) ');\n']);
case 'qBdot'
fprintf(fid,'memset(qBdot_tmp,0,sizeof(realtype)*np);\n');
fprintf(fid,'memset(qBdot_tmp,0,sizeof(realtype)*np*ng);\n');
case 'x0'
fprintf(fid,['memset(x0_tmp,0,sizeof(realtype)*' num2str(nx) ');\n']);
case 'dx0'
Expand Down Expand Up @@ -112,7 +112,7 @@ function printLocalVars(this,model,fid)
case 'deltaxB'
fprintf(fid,['memset(deltaxB,0,sizeof(realtype)*' num2str(nx) ');\n']);
case 'deltaqB'
fprintf(fid,['memset(deltaqB,0,sizeof(realtype)*np);\n']);
fprintf(fid,['memset(deltaqB,0,sizeof(realtype)*np*ng);\n']);
case 'deltasx'
fprintf(fid,['memset(deltasx,0,sizeof(realtype)*' num2str(nx) '*np);\n']);
case 'stau'
Expand Down
2 changes: 1 addition & 1 deletion @amifun/writeCcode.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ function writeCcode(this,model,fid)
if(any(any(nonzero)))
for iy = 1:model.nytrue
fprintf(fid,['if(!mxIsNaN(my[' num2str(iy-1) '*nt+it])){\n']);
tmpfun.sym = this.sym(iy,:);
tmpfun.sym = permute(this.sym(iy,:,:),[2,3,1]);
tmpfun.gccode(model,fid);
fprintf(fid,'}\n');
end
Expand Down
17 changes: 16 additions & 1 deletion @amifun/writeCcode_sensi.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ function writeCcode_sensi(this,model,fid)
% void

np = model.np;
ng = model.ng;

nonzero_idx = find(this.sym);
nonzero = zeros(size(this.sym));
Expand All @@ -18,7 +19,7 @@ function writeCcode_sensi(this,model,fid)
if(strcmp(this.funstr,'deltaqB'))
if(any(nonzero))
tmpfun = this;
for ip=1:np
for ip=1:np*ng
if(nonzero(ip))
fprintf(fid,[' case ' num2str(ip-1) ': {\n']);
tmpfun.sym = this.sym(ip,:);
Expand Down Expand Up @@ -54,6 +55,20 @@ function writeCcode_sensi(this,model,fid)
end
end
end
elseif(strcmp(this.funstr,'qBdot'))
nonzero = this.sym ~=0;
if(any(any(nonzero)))
tmpfun = this;
for ip=1:np
if(any(nonzero(:,ip)))
fprintf(fid,[' case ' num2str(ip-1) ': {\n']);
tmpfun.sym = this.sym(:,ip);
tmpfun.writeCcode(model,fid);
fprintf(fid,'\n');
fprintf(fid,' } break;\n\n');
end
end
end
else
nonzero = this.sym ~=0;
if(any(any(nonzero)))
Expand Down
Loading

0 comments on commit 9f4aa72

Please sign in to comment.