https://github.com/Data2Dynamics/d2d
Tip revision: 4c697ee160f4d91719a65cfcfe3e3b50395b6656 authored by Helge Hass on 29 July 2015, 07:45:33 UTC
new C code version
new C code version
Tip revision: 4c697ee
arParseModel.m
% Parse Model and write function files
%
% arParseModel(forceParsing)
% forceParsing: [false]
%
% Copyright Andreas Raue 2011 (andreas.raue@fdm.uni-freiburg.de)
function arParseModel(forceParsing)
error(['The function arParseModel is deprecated. Please use arCompileAll instead. ' ...
'See https://bitbucket.org/d2d-development/d2d-software/wiki/First%20steps for description of work flow.']);
global ar
if(isempty(ar))
error('please initialize by arInit')
end
if(~exist('forceParsing','var'))
forceParsing = false;
end
usePool = exist('gcp','file')>0 && ~isempty(gcp('nocreate'));
checksum_global = addToCheckSum(ar.info.c_version_code);
for m=1:length(ar.model)
fprintf('\n');
% parse model
arParseODE(m);
% extract & parse conditions
ar.model(m).condition = [];
if(isfield(ar.model(m), 'data'))
for d=1:length(ar.model(m).data)
% conditions checksum
qdynparas = ismember(ar.model(m).data(d).p, ar.model(m).px) | ... %R2013a compatible
ismember(ar.model(m).data(d).p, ar.model(m).data(d).pu); %R2013a compatible
checksum_cond = addToCheckSum(ar.model(m).data(d).fu);
checksum_cond = addToCheckSum(ar.model(m).p, checksum_cond);
checksum_cond = addToCheckSum(ar.model(m).fv, checksum_cond);
checksum_cond = addToCheckSum(ar.model(m).N, checksum_cond);
checksum_cond = addToCheckSum(ar.model(m).cLink, checksum_cond);
checksum_cond = addToCheckSum(ar.model(m).data(d).fp(qdynparas), checksum_cond);
checkstr_cond = getCheckStr(checksum_cond);
% data checksum
checksum_data = addToCheckSum(ar.model(m).data(d).fu);
checksum_data = addToCheckSum(ar.model(m).data(d).p, checksum_data);
checksum_data = addToCheckSum(ar.model(m).data(d).fy, checksum_data);
checksum_data = addToCheckSum(ar.model(m).data(d).fystd, checksum_data);
checksum_data = addToCheckSum(ar.model(m).data(d).fp, checksum_data);
checkstr_data = getCheckStr(checksum_data);
ar.model(m).data(d).checkstr = checkstr_data;
ar.model(m).data(d).fkt = [ar.model(m).data(d).name '_' checkstr_data];
cindex = -1;
for c=1:length(ar.model(m).condition)
if(strcmp(checkstr_cond, ar.model(m).condition(c).checkstr))
cindex = c;
end
end
% global checksum
if(isempty(checksum_global))
checksum_global = addToCheckSum(ar.model(m).data(d).fkt);
else
checksum_global = addToCheckSum(ar.model(m).data(d).fkt, checksum_global);
end
if(cindex == -1) % append new condition
cindex = length(ar.model(m).condition) + 1;
ar.model(m).condition(cindex).status = 0;
ar.model(m).condition(cindex).fu = ar.model(m).data(d).fu;
ar.model(m).condition(cindex).fp = ar.model(m).data(d).fp(qdynparas);
ar.model(m).condition(cindex).p = ar.model(m).data(d).p(qdynparas);
ar.model(m).condition(cindex).checkstr = checkstr_cond;
ar.model(m).condition(cindex).fkt = [ar.model(m).name '_' checkstr_cond];
ar.model(m).condition(cindex).dLink = d;
% global checksum
checksum_global = addToCheckSum(ar.model(m).condition(cindex).fkt, checksum_global);
% link data to condition
ar.model(m).data(d).cLink = length(ar.model(m).condition);
% for multiple shooting
if(isfield(ar.model(m).data(d), 'ms_index') && ~isempty(ar.model(m).data(d).ms_index))
ar.model(m).condition(cindex).ms_index = ...
ar.model(m).data(d).ms_index;
ar.model(m).condition(cindex).ms_snip_index = ...
ar.model(m).data(d).ms_snip_index;
ar.model(m).condition(cindex).ms_snip_start = ar.model(m).data(d).tLim(1);
end
else
% link data to condition
ar.model(m).condition(cindex).dLink(end+1) = d;
ar.model(m).data(d).cLink = cindex;
% for multiple shooting
if(isfield(ar.model(m).data(d), 'ms_index') && ~isempty(ar.model(m).data(d).ms_index))
ar.model(m).condition(cindex).ms_index(end+1) = ...
ar.model(m).data(d).ms_index;
ar.model(m).condition(cindex).ms_snip_index(end+1) = ...
ar.model(m).data(d).ms_snip_index;
ar.model(m).condition(cindex).ms_snip_start(end+1) = ar.model(m).data(d).tLim(1);
end
end
end
% skip parse conditions
doskip = nan(1,length(ar.model(m).condition));
for c=1:length(ar.model(m).condition)
doskip(c) = ~forceParsing && exist(['./Compiled/' ar.info.c_version_code '/' ar.model(m).condition(c).fkt '.c'],'file');
end
% parse conditions
config = ar.config;
model.name = ar.model(m).name;
model.fv = ar.model(m).fv;
model.px0 = ar.model(m).px0;
model.sym = ar.model(m).sym;
model.t = ar.model(m).t;
model.x = ar.model(m).x;
model.u = ar.model(m).u;
model.us = ar.model(m).us;
model.xs = ar.model(m).xs;
model.vs = ar.model(m).vs;
model.N = ar.model(m).N;
condition = ar.model(m).condition;
if(usePool)
parfor c=1:length(ar.model(m).condition)
condition_new(c) = arParseCondition(config, model, condition(c), m, c, doskip(c));
end
else
for c=1:length(ar.model(m).condition)
condition_new(c) = arParseCondition(config, model, condition(c), m, c, doskip(c)); %#ok<AGROW>
end
end
ar.model(m).condition = condition_new;
% skip parse data
doskip = nan(1,length(ar.model(m).data));
for d=1:length(ar.model(m).data)
doskip(d) = ~forceParsing && exist(['./Compiled/' ar.info.c_version_code '/' ar.model(m).data(d).fkt '.c'],'file');
end
% parse data
data = ar.model(m).data;
if(usePool)
parfor d=1:length(ar.model(m).data)
c = data(d).cLink;
data_new(d) = arParseOBS(config, model, condition_new(c), data(d), m, c, d, doskip(d)); %#ok<PFBNS>
end
else
for d=1:length(ar.model(m).data)
c = data(d).cLink;
data_new(d) = arParseOBS(config, model, condition_new(c), data(d), m, c, d, doskip(d)); %#ok<AGROW>
end
end
ar.model(m).data = data_new;
else
qdynparas = ismember(ar.model(m).p, ar.model(m).px) | ... %R2013a compatible
ismember(ar.model(m).p, ar.model(m).pu); %R2013a compatible
% conditions checksum
checksum_cond = addToCheckSum(ar.model(m).fu);
checksum_cond = addToCheckSum(ar.model(m).p(qdynparas), checksum_cond);
checksum_cond = addToCheckSum(ar.model(m).fv, checksum_cond);
checksum_cond = addToCheckSum(ar.model(m).N, checksum_cond);
checksum_cond = addToCheckSum(ar.model(m).cLink, checksum_cond);
checksum_cond = addToCheckSum(ar.model(m).fp, checksum_cond);
% append condition
cindex = 1;
ar.model(m).condition(cindex).status = 0;
ar.model(m).condition(cindex).fu = ar.model(m).fu;
ar.model(m).condition(cindex).fp = ar.model(m).fp(qdynparas);
ar.model(m).condition(cindex).p = ar.model(m).p(qdynparas);
ar.model(m).condition(cindex).checkstr = getCheckStr(checksum_cond);
ar.model(m).condition(cindex).fkt = [ar.model(m).name '_' ar.model(m).condition(cindex).checkstr];
ar.model(m).condition(cindex).dLink = [];
% global checksum
if(isempty(checksum_global))
checksum_global = addToCheckSum(ar.model(m).condition(cindex).fkt);
else
checksum_global = addToCheckSum(ar.model(m).condition(cindex).fkt, checksum_global);
end
% skip parse conditions
doskip = nan(1,length(ar.model(m).condition));
for c=1:length(ar.model(m).condition)
doskip(c) = ~forceParsing && exist(['./Compiled/' ar.info.c_version_code '/' ar.model(m).condition(c).fkt '.c'],'file');
end
% parse conditions
config = ar.config;
model.name = ar.model(m).name;
model.fv = ar.model(m).fv;
model.px0 = ar.model(m).px0;
model.sym = ar.model(m).sym;
model.t = ar.model(m).t;
model.x = ar.model(m).x;
model.u = ar.model(m).u;
model.us = ar.model(m).us;
model.xs = ar.model(m).xs;
model.vs = ar.model(m).vs;
model.N = ar.model(m).N;
condition = ar.model(m).condition;
parfor c=1:length(ar.model(m).condition)
condition_new(c) = arParseCondition(config, model, condition(c), m, c, doskip(c));
end
ar.model(m).condition = condition_new;
% plot setup
if(~isfield(ar.model(m), 'plot'))
ar.model(m).plot(1).name = ar.model(m).name;
else
ar.model(m).plot(end+1).name = ar.model(m).name;
end
ar.model(m).plot(end).doseresponse = false;
ar.model(m).plot(end).dLink = 0;
ar.model(m).plot(end).ny = 0;
ar.model(m).plot(end).condition = {};
end
end
ar.checkstr = getCheckStr(checksum_global);
ar.fkt = ['arSimuCalcFun_' ar.checkstr];
% ODE
function arParseODE(m)
global ar
fprintf('parsing model m%i, %s...', m, ar.model(m).name);
% make short strings
ar.model(m).xs = {};
ar.model(m).us = {};
ar.model(m).vs = {};
for j=1:length(ar.model(m).x)
ar.model(m).xs{j} = sprintf('x[%i]',j);
end
fprintf('x=%i, ', length(ar.model(m).xs));
for j=1:length(ar.model(m).u)
ar.model(m).us{j} = sprintf('u[%i]',j);
end
fprintf('u=%i, ', length(ar.model(m).u));
for j=1:length(ar.model(m).fv)
ar.model(m).vs{j} = sprintf('v[%i]',j);
end
fprintf('v=%i, ', length(ar.model(m).fv));
% make syms
ar.model(m).sym.x = sym(ar.model(m).x);
ar.model(m).sym.xs = sym(ar.model(m).xs);
ar.model(m).sym.px0 = sym(ar.model(m).px0);
ar.model(m).sym.u = sym(ar.model(m).u);
ar.model(m).sym.us = sym(ar.model(m).us);
ar.model(m).sym.vs = sym(ar.model(m).vs);
ar.model(m).sym.fv = sym(ar.model(m).fv);
% compartment volumes
if(~isempty(ar.model(m).pc))
% make syms
ar.model(m).sym.pc = sym(ar.model(m).pc);
ar.model(m).sym.C = sym(ones(size(ar.model(m).N)));
if(~isfield(ar.model(m),'isAmountBased') || ~ar.model(m).isAmountBased)
for j=1:size(ar.model(m).N,1) % for every species j
qinfluxwitheducts = ar.model(m).N(j,:) > 0 & sum(ar.model(m).N < 0,1) > 0;
eductcompartment = zeros(size(qinfluxwitheducts));
for jj=find(qinfluxwitheducts)
eductcompartment(jj) = unique(ar.model(m).cLink(ar.model(m).N(:,jj)<0)); %R2013a compatible
end
cfaktor = sym(ones(size(qinfluxwitheducts)));
for jj=find(qinfluxwitheducts & eductcompartment~=ar.model(m).cLink(j))
cfaktor(jj) = ar.model(m).sym.pc(eductcompartment(jj)) / ...
ar.model(m).sym.pc(ar.model(m).cLink(j));
end
ar.model(m).sym.C(j,:) = transpose(cfaktor);
end
else
for j=1:size(ar.model(m).N,1) % for every species j
ar.model(m).sym.C(j,:) = ar.model(m).sym.C(j,:) / ar.model(m).sym.pc(ar.model(m).cLink(j));
end
end
else
ar.model(m).sym.C = sym(ones(size(ar.model(m).N)));
end
% derivatives
if(~isempty(ar.model(m).sym.fv))
ar.model(m).sym.dfvdx = jacobian(ar.model(m).sym.fv, ar.model(m).sym.x);
if(~isempty(ar.model(m).sym.us))
ar.model(m).sym.dfvdu = jacobian(ar.model(m).sym.fv, ar.model(m).sym.u);
else
ar.model(m).sym.dfvdu = sym(ones(length(ar.model(m).sym.fv), 0));
end
else
ar.model(m).sym.dfvdx = sym(ones(0, length(ar.model(m).sym.x)));
ar.model(m).sym.dfvdu = sym(ones(0, length(ar.model(m).sym.u)));
end
ar.model(m).qdvdx_nonzero = logical(ar.model(m).sym.dfvdx~=0);
ar.model(m).qdvdu_nonzero = logical(ar.model(m).sym.dfvdu~=0);
tmpsym = ar.model(m).sym.dfvdx;
tmpsym = mysubs(tmpsym, ar.model(m).sym.x, ones(size(ar.model(m).sym.x))/2);
tmpsym = mysubs(tmpsym, ar.model(m).sym.u, ones(size(ar.model(m).sym.u))/2);
tmpsym = mysubs(tmpsym, sym(ar.model(m).p), ones(size(ar.model(m).p))/2);
ar.model(m).qdvdx_negative = double(tmpsym) < 0;
tmpsym = ar.model(m).sym.dfvdu;
tmpsym = mysubs(tmpsym, ar.model(m).sym.x, ones(size(ar.model(m).sym.x))/2);
tmpsym = mysubs(tmpsym, ar.model(m).sym.u, ones(size(ar.model(m).sym.u))/2);
tmpsym = mysubs(tmpsym, sym(ar.model(m).p), ones(size(ar.model(m).p))/2);
ar.model(m).qdvdu_negative = double(tmpsym) < 0;
fprintf('done\n');
% Condition
function condition = arParseCondition(config, model, condition, m, c, doskip)
fprintf('parsing condition m%i c%i, %s (%s)...', m, c, model.name, condition.checkstr);
% hard code conditions
condition.sym.p = sym(condition.p);
condition.sym.fp = sym(condition.fp);
condition.sym.fpx0 = sym(model.px0);
condition.sym.fpx0 = mysubs(condition.sym.fpx0, condition.sym.p, condition.sym.fp);
condition.sym.fv = sym(model.fv);
condition.sym.fv = mysubs(condition.sym.fv, condition.sym.p, condition.sym.fp);
condition.sym.fu = sym(condition.fu);
condition.sym.fu = mysubs(condition.sym.fu, condition.sym.p, condition.sym.fp);
condition.sym.C = mysubs(model.sym.C, condition.sym.p, condition.sym.fp);
% predictor
condition.sym.fv = mysubs(condition.sym.fv, sym(model.t), sym('t'));
condition.sym.fu = mysubs(condition.sym.fu, sym(model.t), sym('t'));
% remaining initial conditions
qinitial = ismember(condition.p, model.px0); %R2013a compatible
varlist = cellfun(@symvar, condition.fp(qinitial), 'UniformOutput', false);
condition.px0 = union(vertcat(varlist{:}), [])'; %R2013a compatible
% remaining parameters
varlist = cellfun(@symvar, condition.fp, 'UniformOutput', false);
condition.pold = condition.p;
condition.p = setdiff(setdiff(union(vertcat(varlist{:}), [])', model.x), model.u); %R2013a compatible
if(doskip)
fprintf('skipped\n');
condition.ps = {};
condition.qfu_nonzero = [];
condition.qdvdx_nonzero = [];
condition.qdvdu_nonzero = [];
condition.qdvdp_nonzero = [];
condition.dvdx = {};
condition.dvdu = {};
condition.dvdp = {};
condition.qdfxdx_nonzero = [];
condition.dfxdx = {};
condition.su = {};
condition.sx = {};
condition.qfsv_nonzero = [];
condition.sv = {};
return;
end
% make short strings
condition.ps = {};
for j=1:length(condition.p)
condition.ps{j} = sprintf('p[%i]',j);
end
fprintf('p=%i, ', length(condition.p));
% make syms
condition.sym.p = sym(condition.p);
condition.sym.ps = sym(condition.ps);
condition.sym.px0s = mysubs(sym(condition.px0), ...
condition.sym.p, condition.sym.ps);
% make syms
condition.sym.fv = mysubs(condition.sym.fv, model.sym.x, model.sym.xs);
condition.sym.fv = mysubs(condition.sym.fv, model.sym.u, model.sym.us);
condition.sym.fv = mysubs(condition.sym.fv, condition.sym.p, condition.sym.ps);
condition.sym.fu = mysubs(condition.sym.fu, condition.sym.p, condition.sym.ps);
condition.sym.fpx0 = mysubs(condition.sym.fpx0, condition.sym.p, condition.sym.ps);
% remove zero inputs
condition.qfu_nonzero = logical(condition.sym.fu ~= 0);
if(~isempty(model.sym.us))
condition.sym.fv = mysubs(condition.sym.fv, model.sym.us(~condition.qfu_nonzero), ...
sym(zeros(1,sum(~condition.qfu_nonzero))));
end
% derivatives
if(~isempty(condition.sym.fv))
condition.sym.dfvdx = jacobian(condition.sym.fv, model.sym.xs);
if(~isempty(model.sym.us))
condition.sym.dfvdu = jacobian(condition.sym.fv, model.sym.us);
else
condition.sym.dfvdu = sym(ones(length(condition.sym.fv), 0));
end
condition.sym.dfvdp = jacobian(condition.sym.fv, condition.sym.ps);
else
condition.sym.dfvdx = sym(ones(0, length(model.sym.xs)));
condition.sym.dfvdu = sym(ones(0, length(model.sym.us)));
condition.sym.dfvdp = sym(ones(0, length(condition.sym.ps)));
end
% flux signs
condition.qdvdx_nonzero = logical(condition.sym.dfvdx~=0);
condition.qdvdu_nonzero = logical(condition.sym.dfvdu~=0);
condition.qdvdp_nonzero = logical(condition.sym.dfvdp~=0);
% short terms
condition.dvdx = cell(length(model.vs), length(model.xs));
for j=1:length(model.vs)
for i=1:length(model.xs)
if(condition.qdvdx_nonzero(j,i))
condition.dvdx{j,i} = sprintf('dvdx[%i]', j + (i-1)*length(model.vs));
else
condition.dvdx{j,i} = '0';
end
end
end
condition.sym.dvdx = sym(condition.dvdx);
fprintf('dvdx=%i, ', sum(condition.qdvdx_nonzero(:)));
condition.dvdu = cell(length(model.vs), length(model.us));
for j=1:length(model.vs)
for i=1:length(model.us)
if(condition.qdvdu_nonzero(j,i))
condition.dvdu{j,i} = sprintf('dvdu[%i]', j + (i-1)*length(model.vs));
else
condition.dvdu{j,i} = '0';
end
end
end
condition.sym.dvdu = sym(condition.dvdu);
fprintf('dvdu=%i, ', sum(condition.qdvdu_nonzero(:)));
condition.dvdp = cell(length(model.vs), length(condition.ps));
for j=1:length(model.vs)
for i=1:length(condition.ps)
if(condition.qdvdp_nonzero(j,i))
condition.dvdp{j,i} = sprintf('dvdp[%i]', j + (i-1)*length(model.vs));
else
condition.dvdp{j,i} = '0';
end
end
end
condition.sym.dvdp = sym(condition.dvdp);
fprintf('dvdp=%i, ', sum(condition.qdvdp_nonzero(:)));
% make equations
condition.sym.C = mysubs(condition.sym.C, condition.sym.p, condition.sym.ps);
condition.sym.fx = (model.N .* condition.sym.C) * transpose(model.sym.vs);
% Jacobian dfxdx
if(config.useJacobian)
condition.sym.dfxdx = (model.N .* condition.sym.C) * condition.sym.dvdx;
condition.qdfxdx_nonzero = logical(condition.sym.dfxdx~=0);
for j=1:length(model.xs)
for i=1:length(model.xs)
if(condition.qdfxdx_nonzero(j,i))
condition.dfxdx{j,i} = sprintf('dfxdx[%i]', j + (i-1)*length(model.xs));
else
condition.dfxdx{j,i} = '0';
end
end
end
fprintf('dfxdx=%i, ', sum(condition.qdfxdx_nonzero(:)));
end
% sx sensitivities
if(config.useSensis)
% su
condition.su = cell(length(model.us), length(condition.ps));
for j=1:length(model.us)
for i=1:length(condition.ps)
if(condition.qfu_nonzero(j))
condition.su{j,i} = sprintf('su[%i]', j + (i-1)*length(model.us));
else
condition.su{j,i} = '0';
end
end
end
condition.sym.su = sym(condition.su);
fprintf('su=%i, ', length(condition.ps)*sum(condition.qfu_nonzero(:)));
% input derivatives
if(~isempty(condition.sym.ps))
if(~isempty(condition.sym.fu))
condition.sym.dfudp = ...
jacobian(condition.sym.fu, condition.sym.ps);
else
condition.sym.dfudp = sym(ones(0,length(condition.sym.ps)));
end
% derivatives of step1 (DISABLED)
for j=1:length(model.u)
if(strfind(condition.fu{j}, 'step1('))
condition.sym.dfudp(j,:) = 0;
end
end
% derivatives of step2 (DISABLED)
for j=1:length(model.u)
if(strfind(condition.fu{j}, 'step2('))
condition.sym.dfudp(j,:) = 0;
end
end
% derivatives of spline3
for j=1:length(model.u)
if(strfind(condition.fu{j}, 'spline3('))
for j2=1:length(condition.sym.dfudp(j,:))
ustr = char(condition.sym.dfudp(j,j2));
if(strfind(ustr, 'D([3], spline3)('))
ustr = strrep(ustr, 'D([3], spline3)(', 'Dspline3(');
ustr = strrep(ustr, ')', ', 1)');
elseif(strfind(ustr, 'D([5], spline3)('))
ustr = strrep(ustr, 'D([5], spline3)(', 'Dspline3(');
ustr = strrep(ustr, ')', ', 2)');
elseif(strfind(ustr, 'D([7], spline3)('))
ustr = strrep(ustr, 'D([7], spline3)(', 'Dspline3(');
ustr = strrep(ustr, ')', ', 3)');
end
condition.sym.dfudp(j,j2) = sym(ustr);
end
end
end
% derivatives of spline_pos3
for j=1:length(model.u)
if(strfind(condition.fu{j}, 'spline_pos3('))
for j2=1:length(condition.sym.dfudp(j,:))
ustr = char(condition.sym.dfudp(j,j2));
if(strfind(ustr, 'D([3], spline_pos3)('))
ustr = strrep(ustr, 'D([3], spline_pos3)(', 'Dspline_pos3(');
ustr = strrep(ustr, ')', ', 1)');
elseif(strfind(ustr, 'D([5], spline_pos3)('))
ustr = strrep(ustr, 'D([5], spline_pos3)(', 'Dspline_pos3(');
ustr = strrep(ustr, ')', ', 2)');
elseif(strfind(ustr, 'D([7], spline_pos3)('))
ustr = strrep(ustr, 'D([7], spline_pos3)(', 'Dspline_pos3(');
ustr = strrep(ustr, ')', ', 3)');
end
condition.sym.dfudp(j,j2) = sym(ustr);
end
end
end
% derivatives of spline4
for j=1:length(model.u)
if(strfind(condition.fu{j}, 'spline4('))
for j2=1:length(condition.sym.dfudp(j,:))
ustr = char(condition.sym.dfudp(j,j2));
if(strfind(ustr, 'D([3], spline4)('))
ustr = strrep(ustr, 'D([3], spline4)(', 'Dspline4(');
ustr = strrep(ustr, ')', ', 1)');
elseif(strfind(ustr, 'D([5], spline4)('))
ustr = strrep(ustr, 'D([5], spline4)(', 'Dspline4(');
ustr = strrep(ustr, ')', ', 2)');
elseif(strfind(ustr, 'D([7], spline4)('))
ustr = strrep(ustr, 'D([7], spline4)(', 'Dspline4(');
ustr = strrep(ustr, ')', ', 3)');
elseif(strfind(ustr, 'D([9], spline4)('))
ustr = strrep(ustr, 'D([9], spline4)(', 'Dspline4(');
ustr = strrep(ustr, ')', ', 4)');
end
condition.sym.dfudp(j,j2) = sym(ustr);
end
end
end
% derivatives of spline_pos4
for j=1:length(model.u)
if(strfind(condition.fu{j}, 'spline_pos4('))
for j2=1:length(condition.sym.dfudp(j,:))
ustr = char(condition.sym.dfudp(j,j2));
if(strfind(ustr, 'D([3], spline_pos4)('))
ustr = strrep(ustr, 'D([3], spline_pos4)(', 'Dspline_pos4(');
ustr = strrep(ustr, ')', ', 1)');
elseif(strfind(ustr, 'D([5], spline_pos4)('))
ustr = strrep(ustr, 'D([5], spline_pos4)(', 'Dspline_pos4(');
ustr = strrep(ustr, ')', ', 2)');
elseif(strfind(ustr, 'D([7], spline_pos4)('))
ustr = strrep(ustr, 'D([7], spline_pos4)(', 'Dspline_pos4(');
ustr = strrep(ustr, ')', ', 3)');
elseif(strfind(ustr, 'D([9], spline_pos4)('))
ustr = strrep(ustr, 'D([9], spline_pos4)(', 'Dspline_pos4(');
ustr = strrep(ustr, ')', ', 4)');
end
condition.sym.dfudp(j,j2) = sym(ustr);
end
end
end
% derivatives of spline5
for j=1:length(model.u)
if(strfind(condition.fu{j}, 'spline5('))
for j2=1:length(condition.sym.dfudp(j,:))
ustr = char(condition.sym.dfudp(j,j2));
if(strfind(ustr, 'D([3], spline5)('))
ustr = strrep(ustr, 'D([3], spline5)(', 'Dspline5(');
ustr = strrep(ustr, ')', ', 1)');
elseif(strfind(ustr, 'D([5], spline5)('))
ustr = strrep(ustr, 'D([5], spline5)(', 'Dspline5(');
ustr = strrep(ustr, ')', ', 2)');
elseif(strfind(ustr, 'D([7], spline5)('))
ustr = strrep(ustr, 'D([7], spline5)(', 'Dspline5(');
ustr = strrep(ustr, ')', ', 3)');
elseif(strfind(ustr, 'D([9], spline5)('))
ustr = strrep(ustr, 'D([9], spline5)(', 'Dspline5(');
ustr = strrep(ustr, ')', ', 4)');
elseif(strfind(ustr, 'D([11], spline5)('))
ustr = strrep(ustr, 'D([11], spline5)(', 'Dspline5(');
ustr = strrep(ustr, ')', ', 5)');
end
condition.sym.dfudp(j,j2) = sym(ustr);
end
end
end
% derivatives of spline_pos5
for j=1:length(model.u)
if(strfind(condition.fu{j}, 'spline_pos5('))
for j2=1:length(condition.sym.dfudp(j,:))
ustr = char(condition.sym.dfudp(j,j2));
if(strfind(ustr, 'D([3], spline_pos5)('))
ustr = strrep(ustr, 'D([3], spline_pos5)(', 'Dspline_pos5(');
ustr = strrep(ustr, ')', ', 1)');
elseif(strfind(ustr, 'D([5], spline_pos5)('))
ustr = strrep(ustr, 'D([5], spline_pos5)(', 'Dspline_pos5(');
ustr = strrep(ustr, ')', ', 2)');
elseif(strfind(ustr, 'D([7], spline_pos5)('))
ustr = strrep(ustr, 'D([7], spline_pos5)(', 'Dspline_pos5(');
ustr = strrep(ustr, ')', ', 3)');
elseif(strfind(ustr, 'D([9], spline_pos5)('))
ustr = strrep(ustr, 'D([9], spline_pos5)(', 'Dspline_pos5(');
ustr = strrep(ustr, ')', ', 4)');
elseif(strfind(ustr, 'D([11], spline_pos5)('))
ustr = strrep(ustr, 'D([11], spline_pos5)(', 'Dspline_pos5(');
ustr = strrep(ustr, ')', ', 5)');
end
condition.sym.dfudp(j,j2) = sym(ustr);
end
end
end
% derivatives of spline10
for j=1:length(model.u)
if(strfind(condition.fu{j}, 'spline10('))
for j2=1:length(condition.sym.dfudp(j,:))
ustr = char(condition.sym.dfudp(j,j2));
if(strfind(ustr, 'D([3], spline10)('))
ustr = strrep(ustr, 'D([3], spline10)(', 'Dspline10(');
ustr = strrep(ustr, ')', ', 1)');
elseif(strfind(ustr, 'D([5], spline10)('))
ustr = strrep(ustr, 'D([5], spline10)(', 'Dspline10(');
ustr = strrep(ustr, ')', ', 2)');
elseif(strfind(ustr, 'D([7], spline10)('))
ustr = strrep(ustr, 'D([7], spline10)(', 'Dspline10(');
ustr = strrep(ustr, ')', ', 3)');
elseif(strfind(ustr, 'D([9], spline10)('))
ustr = strrep(ustr, 'D([9], spline10)(', 'Dspline10(');
ustr = strrep(ustr, ')', ', 4)');
elseif(strfind(ustr, 'D([11], spline10)('))
ustr = strrep(ustr, 'D([11], spline10)(', 'Dspline10(');
ustr = strrep(ustr, ')', ', 5)');
elseif(strfind(ustr, 'D([13], spline10)('))
ustr = strrep(ustr, 'D([13], spline10)(', 'Dspline10(');
ustr = strrep(ustr, ')', ', 6)');
elseif(strfind(ustr, 'D([15], spline10)('))
ustr = strrep(ustr, 'D([15], spline10)(', 'Dspline10(');
ustr = strrep(ustr, ')', ', 7)');
elseif(strfind(ustr, 'D([17], spline10)('))
ustr = strrep(ustr, 'D([17], spline10)(', 'Dspline10(');
ustr = strrep(ustr, ')', ', 8)');
elseif(strfind(ustr, 'D([19], spline10)('))
ustr = strrep(ustr, 'D([19], spline10)(', 'Dspline10(');
ustr = strrep(ustr, ')', ', 9)');
elseif(strfind(ustr, 'D([21], spline10)('))
ustr = strrep(ustr, 'D([21], spline10)(', 'Dspline10(');
ustr = strrep(ustr, ')', ', 10)');
end
condition.sym.dfudp(j,j2) = sym(ustr);
end
end
end
% derivatives of spline_pos10
for j=1:length(model.u)
if(strfind(condition.fu{j}, 'spline_pos10('))
for j2=1:length(condition.sym.dfudp(j,:))
ustr = char(condition.sym.dfudp(j,j2));
if(strfind(ustr, 'D([3], spline_pos10)('))
ustr = strrep(ustr, 'D([3], spline_pos10)(', 'Dspline_pos10(');
ustr = strrep(ustr, ')', ', 1)');
elseif(strfind(ustr, 'D([5], spline_pos10)('))
ustr = strrep(ustr, 'D([5], spline_pos10)(', 'Dspline_pos10(');
ustr = strrep(ustr, ')', ', 2)');
elseif(strfind(ustr, 'D([7], spline_pos10)('))
ustr = strrep(ustr, 'D([7], spline_pos10)(', 'Dspline_pos10(');
ustr = strrep(ustr, ')', ', 3)');
elseif(strfind(ustr, 'D([9], spline_pos10)('))
ustr = strrep(ustr, 'D([9], spline_pos10)(', 'Dspline_pos10(');
ustr = strrep(ustr, ')', ', 4)');
elseif(strfind(ustr, 'D([11], spline_pos10)('))
ustr = strrep(ustr, 'D([11], spline_pos10)(', 'Dspline_pos10(');
ustr = strrep(ustr, ')', ', 5)');
elseif(strfind(ustr, 'D([13], spline_pos10)('))
ustr = strrep(ustr, 'D([13], spline_pos10)(', 'Dspline_pos10(');
ustr = strrep(ustr, ')', ', 6)');
elseif(strfind(ustr, 'D([15], spline_pos10)('))
ustr = strrep(ustr, 'D([15], spline_pos10)(', 'Dspline_pos10(');
ustr = strrep(ustr, ')', ', 7)');
elseif(strfind(ustr, 'D([17], spline_pos10)('))
ustr = strrep(ustr, 'D([17], spline_pos10)(', 'Dspline_pos10(');
ustr = strrep(ustr, ')', ', 8)');
elseif(strfind(ustr, 'D([19], spline_pos10)('))
ustr = strrep(ustr, 'D([19], spline_pos10)(', 'Dspline_pos10(');
ustr = strrep(ustr, ')', ', 9)');
elseif(strfind(ustr, 'D([21], spline_pos10)('))
ustr = strrep(ustr, 'D([21], spline_pos10)(', 'Dspline_pos10(');
ustr = strrep(ustr, ')', ', 10)');
end
condition.sym.dfudp(j,j2) = sym(ustr);
end
end
end
end
% sx
condition.sx = cell(length(model.xs), length(condition.ps));
for j=1:length(model.xs)
for i=1:length(condition.ps)
condition.sx{j,i} = sprintf('sx[%i]', j);
end
end
condition.sym.sx = sym(condition.sx);
condition.sym.fsv = condition.sym.dvdx * condition.sym.sx + ...
condition.sym.dvdu * condition.sym.su + condition.sym.dvdp;
% sv
condition.qfsv_nonzero = logical(condition.sym.fsv ~= 0);
condition.sv = cell(length(model.vs), length(condition.ps));
for j=1:length(model.vs)
for i=1:length(condition.ps)
if(condition.qfsv_nonzero(j,i))
condition.sv{j,i} = sprintf('sv[%i]', j);
else
condition.sv{j,i} = '0';
end
end
end
condition.sym.sv = sym(condition.sv);
fprintf('sv=%i, ', sum(condition.qfsv_nonzero(:)));
if(config.useSensiRHS)
condition.sym.fsx = (model.N .* condition.sym.C) * condition.sym.sv;
fprintf('sx=%i... ', numel(condition.sym.fsx));
else
fprintf('sx=skipped ');
end
% sx initials
if(~isempty(condition.sym.fpx0))
condition.sym.fsx0 = jacobian(condition.sym.fpx0, condition.sym.ps);
else
condition.sym.fsx0 = sym(ones(0, length(condition.sym.ps)));
end
% steady state sensitivities
condition.sym.dfxdp = (model.N .* condition.sym.C) * (condition.sym.dvdp + ...
condition.sym.dvdx*condition.sym.fsx0 + ...
condition.sym.dvdu * condition.sym.dfudp);
end
fprintf('done\n');
% OBS
function data = arParseOBS(config, model, condition, data, m, c, d, doskip)
fprintf('parsing data m%i d%i -> c%i, %s (%s)...', m, d, c, data.name, data.checkstr);
% hard code conditions
data.sym.p = sym(data.p);
data.sym.fp = sym(data.fp);
data.sym.fy = sym(data.fy);
data.sym.fy = mysubs(data.sym.fy, data.sym.p, data.sym.fp);
data.sym.fystd = sym(data.fystd);
data.sym.fystd = mysubs(data.sym.fystd, data.sym.p, data.sym.fp);
data.sym.fu = sym(condition.fu);
data.sym.fu = mysubs(data.sym.fu, data.sym.p, data.sym.fp);
data.qfu_nonzero = logical(data.sym.fu ~= 0);
% predictor
data.sym.fu = mysubs(data.sym.fu, sym(model.t), sym('t'));
data.sym.fy = mysubs(data.sym.fy, sym(model.t), sym('t'));
data.sym.fystd = mysubs(data.sym.fystd, sym(model.t), sym('t'));
% remaining parameters
varlist = cellfun(@symvar, data.fp, 'UniformOutput', false);
data.pold = data.p;
data.p = setdiff(setdiff(union(vertcat(varlist{:}), [])', model.x), model.u); %R2013a compatible
if(doskip)
fprintf('skipped\n');
data.ps = {};
data.ys = [];
data.qu_measured = [];
data.qx_measured = [];
data.dfydxnon0 = {};
data.sx = {};
data.sy = {};
return;
end
% make short strings
for j=1:length(data.p)
data.ps{j} = sprintf('p[%i]',j);
end
data.ys = {};
for j=1:length(data.y)
data.ys{j} = sprintf('y[%i]',j);
end
fprintf('y=%i, p=%i, ', length(data.y), length(data.p));
% make syms
data.sym.p = sym(data.p);
data.sym.ps = sym(data.ps);
data.sym.y = sym(data.y);
data.sym.ys = sym(data.ys);
% substitute
data.sym.fy = mysubs(data.sym.fy, ...
model.sym.x, model.sym.xs);
data.sym.fy = mysubs(data.sym.fy, ...
model.sym.u, model.sym.us);
data.sym.fy = mysubs(data.sym.fy, ...
data.sym.p, data.sym.ps);
data.sym.fystd = mysubs(data.sym.fystd, ...
model.sym.x, model.sym.xs);
data.sym.fystd = mysubs(data.sym.fystd, ...
model.sym.u, model.sym.us);
data.sym.fystd = mysubs(data.sym.fystd, ...
data.sym.y, data.sym.ys);
data.sym.fystd = mysubs(data.sym.fystd, ...
data.sym.p, data.sym.ps);
% remove zero inputs
% data.sym.fy = mysubs(data.sym.fy, model.sym.us(~condition.qfu_nonzero), ...
% sym(zeros(1,sum(~condition.qfu_nonzero))));
% data.sym.fystd = mysubs(data.sym.fystd, model.sym.us(~condition.qfu_nonzero), ...
% sym(zeros(1,sum(~condition.qfu_nonzero))));
% derivatives fy
if(~isempty(data.sym.fy))
if(~isempty(model.sym.us))
data.sym.dfydu = jacobian(data.sym.fy, model.sym.us);
else
data.sym.dfydu = sym(ones(length(data.y), 0));
end
if(~isempty(model.x))
data.sym.dfydx = jacobian(data.sym.fy, model.sym.xs);
else
data.sym.dfydx = [];
end
data.sym.dfydp = jacobian(data.sym.fy, data.sym.ps);
else
data.sym.dfydu = [];
data.sym.dfydx = [];
data.sym.dfydp = [];
end
% what is measured ?
data.qu_measured = sum(logical(data.sym.dfydu~=0),1)>0;
data.qx_measured = sum(logical(data.sym.dfydx~=0),1)>0;
% derivatives fystd
if(~isempty(data.sym.fystd))
if(~isempty(model.sym.us))
data.sym.dfystddu = jacobian(data.sym.fystd, model.sym.us);
else
data.sym.dfystddu = sym(ones(length(data.y), 0));
end
if(~isempty(model.x))
data.sym.dfystddx = jacobian(data.sym.fystd, model.sym.xs);
else
data.sym.dfystddx = [];
end
data.sym.dfystddp = jacobian(data.sym.fystd, data.sym.ps);
data.sym.dfystddy = jacobian(data.sym.fystd, data.sym.ys);
else
data.sym.dfystddp = [];
data.sym.dfystddy = [];
data.sym.dfystddx = [];
end
% observed directly and exclusively
data.dfydxnon0 = logical(data.sym.dfydx ~= 0);
if(config.useSensis)
% sx sensitivities
data.sx = {};
for j=1:length(model.xs)
for i=1:length(condition.p)
data.sx{j,i} = sprintf('sx[%i]', j + (i-1)*length(model.xs));
end
end
data.sym.sx = sym(data.sx);
% sy sensitivities
data.sy = {};
for j=1:length(data.sym.fy)
for i=1:length(data.sym.ps)
data.sy{j,i} = sprintf('sy[%i]', j + (i-1)*length(data.sym.fy));
end
end
data.sym.sy = sym(data.sy);
% calculate sy
if(~isempty(data.sym.sy))
data.sym.fsy = data.sym.dfydp;
if(~isempty(condition.p))
qdynpara = ismember(data.p, condition.p); %R2013a compatible
else
qdynpara = false(size(data.p));
end
if(~isempty(condition.p))
tmpfsx = data.sym.dfydx * ...
data.sym.sx;
if(~isfield(condition, 'sym') || ~isfield(condition.sym, 'su'))
% su
condition.su = cell(length(model.us), length(condition.p));
for j=1:length(model.us)
for i=1:length(condition.p)
if(data.qfu_nonzero(j))
condition.su{j,i} = sprintf('su[%i]', j + (i-1)*length(model.us));
else
condition.su{j,i} = '0';
end
end
end
condition.sym.su = sym(condition.su);
fprintf('su=%i, ', length(condition.p)*sum(data.qfu_nonzero(:)));
end
tmpfsu = data.sym.dfydu * ...
condition.sym.su;
if(~isempty(model.x))
data.sym.fsy(:,qdynpara) = data.sym.fsy(:,qdynpara) + tmpfsx + tmpfsu;
else
data.sym.fsy(:,qdynpara) = data.sym.fsy(:,qdynpara) + tmpfsu;
end
end
else
data.sym.fsy = [];
end
fprintf('sy=%i, ', sum(logical(data.sym.fsy(:)~=0)));
% calculate systd
if(~isempty(data.sym.sy))
data.sym.fsystd = data.sym.dfystddp + ...
data.sym.dfystddy * data.sym.sy;
if(~isempty(condition.p))
qdynpara = ismember(data.p, condition.p); %R2013a compatible
else
qdynpara = false(size(data.p));
end
if(~isempty(condition.p))
tmpfsx = data.sym.dfystddx * ...
data.sym.sx;
if(~isfield(condition, 'sym') || ~isfield(condition.sym, 'su'))
% su
condition.su = cell(length(model.us), length(condition.p));
for j=1:length(model.us)
for i=1:length(condition.p)
if(data.qfu_nonzero(j))
condition.su{j,i} = sprintf('su[%i]', j + (i-1)*length(model.us));
else
condition.su{j,i} = '0';
end
end
end
condition.sym.su = sym(condition.su);
fprintf('su=%i, ', length(condition.p)*sum(data.qfu_nonzero(:)));
end
tmpfsu = data.sym.dfystddu * ...
condition.sym.su;
if(~isempty(model.x))
data.sym.fsystd(:,qdynpara) = data.sym.fsystd(:,qdynpara) + tmpfsx + tmpfsu;
else
data.sym.fsystd(:,qdynpara) = data.sym.fsystd(:,qdynpara) + tmpfsu;
end
end
else
data.sym.fsystd = [];
end
fprintf('systd=%i, ', sum(logical(data.sym.fsystd(:)~=0)));
end
fprintf('done\n');
% better subs
function out = mysubs(in, old, new)
if(~isnumeric(in) && ~isempty(old) && ~isempty(findsym(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
function checksum = addToCheckSum(str, checksum)
algs = {'MD2','MD5','SHA-1','SHA-256','SHA-384','SHA-512'};
if(nargin<2)
checksum = java.security.MessageDigest.getInstance(algs{2});
end
if(iscell(str))
for j=1:length(str)
checksum = addToCheckSum(str{j}, checksum);
end
else
if(~isempty(str))
checksum.update(uint8(str(:)));
end
end
function checkstr = getCheckStr(checksum)
h = typecast(checksum.digest,'uint8');
checkstr = dec2hex(h)';
checkstr = checkstr(:)';
clear checksum