Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bugfixng with the Jacobian generation #17

Closed
wants to merge 18 commits into from
Closed

Bugfixng with the Jacobian generation #17

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
ec9f899
Merge branch 'master' into stable
Jul 1, 2016
0bf1309
Merge branch 'master' into develop
Jul 1, 2016
d3986d5
fix m-function file generation
Jul 8, 2016
9d2a93a
SBML rate rules, hooray!
Jul 8, 2016
856888f
fix for piecewise functions
Jul 8, 2016
1a11ce1
fixed or, added xor
Jul 8, 2016
8790d09
fix dependencies in generated M functions
Jul 8, 2016
5e7fd72
fix for certain function definitions
Jul 8, 2016
392a5a1
rate rules for parameters, go get it gurl!
Jul 8, 2016
6dfd7f1
added my mfun
paulstapor Jul 11, 2016
8dfdd9a
removed my old mfun
paulstapor Jul 11, 2016
27e2936
tried to pull from icb_dcm stable, see what happens
paulstapor Jul 11, 2016
00c1b43
new mfun
paulstapor Jul 11, 2016
633f0f9
fixes for state dependent event ouptuts
Jul 12, 2016
b3cf36c
Merge branch 'develop' into stable
Jul 12, 2016
00711f1
fixed a bug for the Jacobian file generation
paulstapor Jul 12, 2016
c61fea3
changed little thing in mfun
paulstapor Jul 12, 2016
e32d90c
hopefully now on stable
paulstapor Jul 14, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
126 changes: 94 additions & 32 deletions @SBMLode/importSBML.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,13 +48,11 @@ function importSBML(this,modelname)
error('DAEs currently not supported!');
end

if(any(strcmp(rule_types,'SBML_RATE_RULE')))
%custom reate laws TBD
error('Sorry, custom rate laws currently not supported!');
end

rulevars = sym({model.rule.variable});
rulemath = sym({model.rule.formula});
% remove rate rules
rulevars = rulevars(not(strcmp({model.rule.typecode},'SBML_RATE_RULE')));
rulemath = rulemath(not(strcmp({model.rule.typecode},'SBML_RATE_RULE')));
repeat_idx = ismember(rulevars,symvar(rulemath));
while(any(repeat_idx))
rulemath= subs(rulemath,rulevars,rulemath);
Expand Down Expand Up @@ -120,12 +118,9 @@ function importSBML(this,modelname)
% apply rules
applyRule(this,model,'initState',rulevars,rulemath)

this.knom = double(this.initState(cond_idx));

% remove constant species
const_idx = logical([model.species.constant]) & not(cond_idx);
constant_sym = this.state(const_idx);
constants = this.initState(const_idx);
while(any(ismember(symvar(this.initState),this.state)))
this.initState = subs(this.initState,this.state,this.initState);
end

%% PARAMETERS

Expand All @@ -143,6 +138,16 @@ function importSBML(this,modelname)

np = length(this.param);

%% CONSTANTS

this.knom = double(subs(this.initState(cond_idx),parameter_sym,parameter_val));

% remove constant species
const_idx = logical([model.species.constant]) & not(cond_idx);
constant_sym = this.state(const_idx);
constants = this.initState(const_idx);



%% REACTIONS

Expand All @@ -152,13 +157,10 @@ function importSBML(this,modelname)

kLaw = cellfun(@(x) x.math,{model.reaction.kineticLaw},'UniformOutput',false);

if(any(cell2mat(strfind(kLaw,'factorial'))))
error('Sorry, factorial functions are not supported at the moment!')
end
checkIllegalFunctions(kLaw);

if(any(cell2mat(strfind(kLaw,'ceil'))))
error('Sorry, ceil functions are not supported at the moment!')
end
kLaw = replaceDiscontinuousFunctions(kLaw);
kLaw = replaceReservedFunctions(kLaw);

this.flux = sym(kLaw);
this.flux = this.flux(:);
Expand Down Expand Up @@ -231,6 +233,27 @@ function importSBML(this,modelname)
this.xdot = sym(zeros(size(this.state)));
end

%% RATE RULES

for irule = 1:length(model.rule)
if(strcmp(model.rule(irule).typecode,'SBML_RATE_RULE'))
state_rate_idx = find(this.state == sym(model.rule(irule).variable));
param_rate_idx = find(parameter_sym == sym(model.rule(irule).variable));
if(~isempty(state_rate_idx))
this.xdot(fstate_rate_idx) = sym(model.rule(irule).formula);
elseif(~isempty(param_rate_idx))
this.state = [this.state; this.param(param_rate_idx)];
this.xdot = [this.xdot; sym(model.rule(irule).formula)];
this.initState = [this.initState; parameter_val(param_rate_idx)];
this.volume = [this.volume; 1];
nx = nx + 1;
parameter_val(param_rate_idx) = [];
parameter_sym(param_rate_idx) = [];
this.param(param_rate_idx) = [];
np = np - 1;
end
end
end
%% CONVERSION FACTORS/VOLUMES

fprintf('converting to concentrations ...\n')
Expand Down Expand Up @@ -279,7 +302,7 @@ function importSBML(this,modelname)
assignments = sym(cat(2,tmp{:}));

tmp = cellfun(@(x) {x.math},{model.event(ievent).eventAssignment},'UniformOutput',false);
assignments_math = sym(cat(2,tmp{:}));
assignments_math = sym(replaceReservedFunctions(cat(2,tmp{:})));

for iassign = 1:length(assignments)
state_assign_idx = find(assignments(iassign)==species_sym);
Expand All @@ -288,15 +311,15 @@ function importSBML(this,modelname)
bound_assign_idx = find(assignments(iassign)==boundary_sym);

if(np>0 && ~isempty(param_assign_idx))
this.param(assignments_param_idx) = this.param(assignments_param_idx)*heaviside(this.trigger(ievent)) + assignments_math(iassign)*heaviside(-this.trigger(ievent));
this.param(param_assign_idx) = this.param(param_assign_idx)*heaviside(this.trigger(ievent)) + assignments_math(iassign)*heaviside(-this.trigger(ievent));
end

if(nk>0 && ~isempty(assignments_cond_tidx))
conditions(assignments_cond_idx) = conditions(assignments_cond_idx)*heaviside(this.trigger(ievent)) + assignments_math(iassign)*heaviside(-this.trigger(ievent));
if(nk>0 && ~isempty(cond_assign_idx))
conditions(cond_assign_idx) = conditions(cond_assign_idx)*heaviside(this.trigger(ievent)) + assignments_math(iassign)*heaviside(-this.trigger(ievent));
end

if(length(boundaries)>0 && ~isempty(bound_assign_idx))
boundaries(assignments_bound_idx) = conditions(assignments_bound_idx)*heaviside(this.trigger(ievent)) + assignments_math(iassign)*heaviside(-this.trigger(ievent));
boundaries(bound_assign_idx) = conditions(bound_assign_idx)*heaviside(this.trigger(ievent)) + assignments_math(iassign)*heaviside(-this.trigger(ievent));
end

if(length(this.state)>0 && ~isempty(state_assign_idx))
Expand Down Expand Up @@ -326,23 +349,19 @@ function importSBML(this,modelname)
tmpfun = cellfun(@(x) ['fun_' num2str(x)],num2cell(1:length(model.functionDefinition)),'UniformOutput',false);
this.funmath = strrep(this.funmath,{model.functionDefinition.id},tmpfun);
% replace helper functions
this.funmath = strrep(this.funmath,'ge(','am_ge(');
this.funmath = strrep(this.funmath,'gt(','am_gt(');
this.funmath = strrep(this.funmath,'le(','am_le(');
this.funmath = strrep(this.funmath,'lt(','am_lt(');
this.funmath = strrep(this.funmath,'piecewise(','am_piecewise(');
this.funmath = strrep(this.funmath,'max(','am_max(');
this.funmath = strrep(this.funmath,'min(','am_min(');

checkIllegalFunctions(this.funmath);
this.funmath = replaceDiscontinuousFunctions(this.funmath);

this.funmath = strrep(this.funmath,tmpfun,{model.functionDefinition.id});
this.funarg = cellfun(@(x,y) [y '(' strjoin(transpose(x(1:end-1)),',') ')'],lambdas,{model.functionDefinition.id},'UniformOutput',false);
this.funarg = cellfun(@(x,y) [y '(' strjoin(transpose(x(1:end-1)),',') ')'],lambdas,replaceReservedFunctionIDs({model.functionDefinition.id}),'UniformOutput',false);

% make functions available in this file

for ifun = 1:length(this.funmath)
token = regexp(this.funarg(ifun),'\(([0-9\w\,]*)\)','tokens');
start = regexp(this.funarg(ifun),'\(([0-9\w\,]*)\)');
eval([this.funarg{ifun}(1:(start{1}-1)) ' = @(' token{1}{1}{1} ')' this.funmath{ifun} ';']);
eval([replaceReservedFunctions(this.funarg{ifun}(1:(start{1}-1))) ' = @(' token{1}{1}{1} ')' this.funmath{ifun} ';']);
end
end

Expand Down Expand Up @@ -370,7 +389,7 @@ function importSBML(this,modelname)
end

state_vars = [symvar(this.xdot),symvar(this.initState)];
event_vars = [symvar(addToBolus),symvar(this.trigger)];
event_vars = [symvar(this.bolus),symvar(this.trigger)];

% substitute with actual expressions
makeSubs(this,parameter_sym(1:np),this.param);
Expand Down Expand Up @@ -435,3 +454,46 @@ function makeSubs(this,old,new)
this.bolus = subs(this.bolus,old,new);
this.initState = subs(this.initState,old,new);
end

function checkIllegalFunctions(str)

if(any(cell2mat(strfind(str,'factorial'))))
error('Sorry, factorial functions are not supported at the moment!')
end
if(any(cell2mat(strfind(str,'ceil'))))
error('Sorry, ceil functions are not supported at the moment!')
end
if(any(cell2mat(strfind(str,'floor'))))
error('Sorry, floor functions are not supported at the moment!')
end
end


function str = replaceDiscontinuousFunctions(str)
% replace imcompatible piecewise defintion
% execute twice for directly nested calls (overlapping regexp expressions)
for logicalf = {'piecewise','and','or','lt','gt','ge','le','ge','le','xor'}
str = regexprep(str,['^' logicalf{1} '('],['am_' logicalf{1} '(']);
str = regexprep(str,['([\W]+)' logicalf{1} '('],['$1am_' logicalf{1} '(']);
str = regexprep(str,['([\W]+)' logicalf{1} '('],['$1am_' logicalf{1} '(']);
end
end

function str = replaceReservedFunctions(str)
% replace imcompatible piecewise defintion
% execute twice for directly nested calls (overlapping regexp expressions)
for logicalf = {'divide'}
str = regexprep(str,['^' logicalf{1} '('],['am_' logicalf{1} '(']);
str = regexprep(str,['([\W]+)' logicalf{1} '('],['$1am_' logicalf{1} '(']);
str = regexprep(str,['([\W]+)' logicalf{1} '('],['$1am_' logicalf{1} '(']);
end
end

function str = replaceReservedFunctionIDs(str)
% replace imcompatible piecewise defintion
% execute twice for directly nested calls (overlapping regexp expressions)
for logicalf = {'divide'}
str = regexprep(str,['^' logicalf{1} '$'],['am_' logicalf{1} '']);
end
end

12 changes: 9 additions & 3 deletions @amifun/gccode.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@

cstr = ccode(this.sym);
if(~strcmp(cstr(3:4),'t0'))
if(strcmp(this.funstr,'J') || strcmp(this.funstr,'JB') || strcmp(this.funstr,'dJydp') || strcmp(this.funstr,'dJydx') || strcmp(this.funstr,'dydx') || strcmp(this.funstr,'sJy') || strcmp(this.funstr,'M') || strcmp(this.funstr,'dfdx') )
if(strcmp(this.funstr,'J') || strcmp(this.funstr,'JB') || strcmp(this.funstr,'dJydp') || strcmp(this.funstr,'dJydx') || strcmp(this.funstr,'dydx') || strcmp(this.funstr,'dzdx') || strcmp(this.funstr,'sJy') || strcmp(this.funstr,'M') || strcmp(this.funstr,'dfdx') )
cstr = regexprep(cstr,'T\[([0-9]*)\]\[([0-9]*)\]',[this.cvar '[$1+$2*' num2str(size(this.sym,1)) ']']);
else
cstr = regexprep(cstr,'T\[([0-9]*)\]\[0\]',[this.cvar '[$1]']);
Expand Down Expand Up @@ -92,7 +92,7 @@
cstr = regexprep(cstr,[this.cvar '\[([0-9+*]*)\]'],[this.cvar '\[it+nt*\(($1)+ip*' num2str(model.ny) '\)\]']);
elseif(strcmp(this.cvar,'z'))
cstr = regexprep(cstr,[this.cvar '\[([0-9+*]*)\]'],[this.cvar '\[nroots[ie] + nmaxevent*($1)\]']);
elseif(strcmp(this.cvar,'sz') || strcmp(this.cvar,'sroot'))
elseif(strcmp(this.cvar,'sz') || strcmp(this.cvar,'sroot') || strcmp(this.cvar,'sz_tf'))
cstr = regexprep(cstr,[this.cvar '\[([0-9+*]*)\]'],[this.cvar '\[nroots[ie] + nmaxevent*(($1)+ip*' num2str(model.nz) ')\]']);
elseif(strcmp(this.cvar,'s2root'))
cstr = regexprep(cstr,[this.cvar '\[([0-9+*]*)\]'],[this.cvar '\[nroots[ie] + nmaxevent*(($1)+ip*' num2str(model.nz) '*np)\]']);
Expand All @@ -117,7 +117,7 @@
elseif(strcmp(this.cvar,'dJydx'))
cstr = regexprep(cstr, [this.cvar '\[([0-9+*]*)\] = '], [ this.cvar '[it+($1)*nt] = ']);
elseif(strcmp(this.cvar,'dJzdx'))
cstr = regexprep(cstr, [this.cvar '\[([0-9+*]*)\] = '], [ this.cvar '[nroots[ie]+($1)*maxevent] = ']);
cstr = regexprep(cstr, [this.cvar '\[([0-9+*]*)\] = '], [ this.cvar '[nroots[ie]+($1)*nmaxevent] = ']);
end

if(strfind(this.cvar,'Jy'))
Expand Down Expand Up @@ -156,6 +156,9 @@
if(strcmp(this.funstr,'dydx'))
cstr = strrep(cstr,'dydx_tmp','dydx');
end
if(strcmp(this.funstr,'dzdx'))
cstr = strrep(cstr,'dzdx_tmp','dzdx');
end
if(strcmp(this.funstr,'deltax'))
cstr = strrep(cstr,'deltax_tmp','deltax');
end
Expand All @@ -171,6 +174,9 @@
if(strcmp(this.funstr,'dJydx'))
cstr = strrep(cstr,'ydx_tmp','ydx');
end
if(strcmp(this.funstr,'dJzdx'))
cstr = strrep(cstr,'zdx_tmp','zdx');
end
if(strcmp(this.funstr,'x0'))
cstr = regexprep(cstr,'x_([0-9]+)','x0_tmp\[$1\]');
end
Expand Down
2 changes: 2 additions & 0 deletions @amifun/getCVar.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,8 @@
this.cvar = 'M_tmp';
case 'dfdx'
this.cvar = 'dfdx_tmp';
case 'sz_tf'
this.cvar = 'sz';
otherwise
this.cvar = this.funstr;
end
Expand Down
3 changes: 3 additions & 0 deletions @amifun/getSensiFlag.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,9 @@
case 'sz'
this.sensiflag = true;

case 'sz_tf'
this.sensiflag = true;

case 'dtaudp'
this.sensiflag = true;

Expand Down
27 changes: 23 additions & 4 deletions @amifun/writeMcode.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,10 +21,10 @@ function writeMcode(this,model)


if(strcmp(this.funstr,'w') || strcmp(this.funstr,'dwdp') || strcmp(this.funstr,'dwdx'))
fid = fopen(fullfile(model.wrap_path,'models',model.modelname,['calc_' this.funstr '_',model.modelname,'.m']), 'w');
fprintf(fid,['%% calc_' this.funstr '_' model.modelname '.m is an additional matlab file to\n'...
fid = fopen(fullfile(model.wrap_path,'models',model.modelname,[ this.funstr '_',model.modelname,'.m']), 'w');
fprintf(fid,['%% ' this.funstr '_' model.modelname '.m is an additional matlab file to\n'...
'%% compute helping variables for the jacobian with AMICI.\n\n']);
fprintf(fid, ['function ' this.funstr ' = calc_' this.funstr '_' model.modelname '(p, x, k, t)\n\n']);
fprintf(fid, ['function ' this.funstr ' = ' this.funstr '_' model.modelname '(p, x, k, t)\n\n']);

% Assignment of inputs
for l = 1 : numel(model.sym.p)
Expand All @@ -36,6 +36,11 @@ function writeMcode(this,model)
for l = 1 : numel(model.sym.k)
fprintf(fid, ['k_' num2str(l-1) ' = k(' num2str(l) ');\n']);
end
if(~strcmp(this.funstr,'w'))
for l = 1 : numel(model.fun.w.sym)
fprintf(fid, ['w_' num2str(l-1) ' = w(' num2str(l) ');\n']);
end
end
fprintf(fid, '\n');

% Writing and giving output
Expand All @@ -57,12 +62,26 @@ function writeMcode(this,model)
fprintf(fid, 'end');
fclose(fid);
else
if(and(ismember('dwdp',model.fun.(this.funstr).deps),ismember('dwdx',model.fun.(this.funstr).deps)))
mfun(this.sym, 'file', fullfile(model.wrap_path,'models',...
model.modelname,['eval_' this.funstr '_',model.modelname,'.m']), ...
model.modelname,[ this.funstr '_',model.modelname,'.m']), ...
'vars', {model.fun.p.sym, model.fun.x.sym, model.fun.k.sym, ...
model.fun.w.strsym(find(model.fun.w.strsym)), ...
model.fun.dwdp.strsym(find(model.fun.dwdp.strsym)), ...
model.fun.dwdx.strsym(find(model.fun.dwdx.strsym)), 't'});
elseif(ismember('dwdp',model.fun.(this.funstr).deps))
mfun(this.sym, 'file', fullfile(model.wrap_path,'models',...
model.modelname,[ this.funstr '_',model.modelname,'.m']), ...
'vars', {model.fun.p.sym, model.fun.x.sym, model.fun.k.sym, ...
model.fun.w.strsym(find(model.fun.w.strsym)), ...
model.fun.dwdp.strsym(find(model.fun.dwdp.strsym)), 't'});
elseif(ismember('dwdx',model.fun.(this.funstr).deps))
mfun(this.sym, 'file', fullfile(model.wrap_path,'models',...
model.modelname,[ this.funstr '_',model.modelname,'.m']), ...
'vars', {model.fun.p.sym, model.fun.x.sym, model.fun.k.sym, ...
model.fun.w.strsym(find(model.fun.w.strsym)), ...
model.fun.dwdx.strsym(find(model.fun.dwdx.strsym)), 't'});
end
end


Expand Down
Loading