Skip to content

Commit

Permalink
improved logical structure and readability of nCVRestart loop
Browse files Browse the repository at this point in the history
  • Loading branch information
@niklasneubrand
niklasneubrand committed Nov 15, 2023
1 parent a5b5acf commit 95499dd
Showing 1 changed file with 106 additions and 66 deletions.
172 changes: 106 additions & 66 deletions arFramework3/Subfunctions/arCalcMerit.m
View file
Original file line number Diff line number Diff line change
@@ -1,91 +1,91 @@
% [ar2] = arCalcMerit([sensi], [pTrial], [dynamics], [doSimu])
% [ar2] = arCalcMerit([ar], [sensi], [pTrial], [dynamics], [doSimu])
%
%
% This function updates the fields representing objective functions (and
% parts of it) used for fitting for the currently chosen parameters.
%
%
% sensi propagate sensitivities [false]
% this argument is passed to arSimu
% pTrial trial parameter of fitting
% dynamics force evaluation of dynamics [false]
% doSimu should arSimu be called [true]
% ar d2d model/data structure
%
%
% ar2 d2d model/data structure with updated objective functions
%
%
% This function replaces arChi2 and is called by arFit. The fact that the
% function accepts serval calls is because of historical developments. It
% is not beautiful but required to not break compatibility.
%
% is not beautiful but required to not break compatibility.
%
% This function calls
% - arSimu which in turn calls arCalcRes
% - arCollectRes
%
%
% Possible calls:
% >> arCalcMerit then:
% qglobalar = true
% sensi = true
%
%
% >> arCalcMerit(ar) here, the global ar is overwritten by the argument
% qglobalar = false
% sensi = true
%
%
% >> arCalcMerit(ar,sensi)
% >> arCalcMerit(sensi)
% like arCalcMerit, but sensi can be set
%
% like arCalcMerit, but sensi can be set
%
% >> arCalcMerit(ar,sensi,ptrial)
% >> arCalcMerit(sensi,ptrial)
% like arCalcMerit(ar,sensi), but
% like arCalcMerit(ar,sensi), but
% silent = true
% ar.p is set to ptrial
% this is one possiblity to set 'silent' to true. The other
% one is:
% >> arCalcMerit(sensi,[])
%
%
% >> arCalcMerit(sensi,ptrial,dynamics)
% >> arCalcMerit(ar,sensi,ptrial,dynamics)
% dynamics is passed to arSimu, otherwise arSimu is called
% without this argument, i.e. with its default
%
% without this argument, i.e. with its default
%
% >> arCalcMerit(sensi,ptrial,dynamics,doSimu)
% >> arCalcMerit(ar,sensi,ptrial,dynamics,doSimu)
% doSimu can be set to false, then the residuals are calculated
% without updating the model trajectories (e.g. if ar.qFit or ar.
% model.data.qFit) has been changed.
%
%
% Example (call used by arFit)
% ar = arCalcMerit(ar, true, ar.p(ar.qFit==1));
%
%
% See also arGetMerit, arChi2


function varargout = arCalcMerit(varargin)

global ar
global ar %#ok<*GVMIS>

nargs = nargin;
% The possiblity providing ar as an argument and to use of qglobalar==0 is
% obsolete because the gloal "ar" is overwritten anyway in arSimu
% obsolete because the gloal "ar" is overwritten anyway in arSimu
% Implementation due to backwards compability:
if nargs>0 && isstruct(varargin{1})
ar = varargin{1};
varargin = varargin(2:end); % changing the meaning of varargin is not nicely implemented
nargs = nargs-1;
nargs = nargs-1;
qglobalar = false;
else
qglobalar = true;
end

if nargs>=1 && ~isempty(varargin{1})
sensi = varargin{1} && ar.config.useSensis;
else
else
sensi = true && ar.config.useSensis;
end

if nargs>=2 && ~isempty(varargin{2})
pTrial = varargin{2};
ar.p(ar.qFit==1) = pTrial + 0;
ar.p(ar.qFit==1) = pTrial + 0;
silent = true;
elseif nargs==2 && isempty(varargin{2})
silent = true;
Expand All @@ -106,15 +106,16 @@
end

if(~isfield(ar, 'fevals'))
ar.fevals = 0;
ar.fevals = 0;
end
ar.fevals = ar.fevals + 1;


atol = ar.config.atol;
rtol = ar.config.rtol;
maxsteps = ar.config.maxsteps;
qPositiveX = cell(1,length(ar.model));
qPositiveX = {ar.model(:).qPositiveX};
allowNegativeX = 0;

if(~isfield(ar.config, 'nCVRestart') || isnan(ar.config.nCVRestart))
nCVRestart = 10;
Expand All @@ -123,66 +124,105 @@
end

for i = 1:nCVRestart

if ~doSimu
has_error = false;
break
end

try
if doSimu
if(qglobalar) % since ar is overwritten anyway in arSimu, the possiblity to use of qglobalar obsolete
arSimu(sensi, false, dynamics);
else
ar = arSimu(ar, sensi, false, dynamics);
end
if(qglobalar) % since ar is overwritten anyway in arSimu, the possiblity to use of qglobalar obsolete
arSimu(sensi, false, dynamics);
else
ar = arSimu(ar, sensi, false, dynamics);
end

has_error = false;
break

catch error_id

has_error = true;
if(~silent)
arFprintf(1, 'Integration error.\n');
disp(error_id.message);
end
if nCVRestart > 1
if strcmp(error_id.identifier,'MATLAB:UndefinedFunction')

% in the following cases we have no fix and exit the loop:
% undefined function
if strcmp(error_id.identifier,'MATLAB:UndefinedFunction')
break
end
% no typical CVODES error flag -> there probably was another error
% (there are also the CVODES flags 1, 2, and 99
% but currently they are not handled specifically)
for m=1:length(ar.model)
if all([ar.model(m).condition(:).status] >= 0)
break
else
error_printed = 0;
for m=1:length(ar.model)
for c=1:length(ar.model(m).condition)
if(ar.model(m).condition(c).status==-1)
% CV_TOO_MUCH_WORK
if(isempty(qPositiveX{m}))
qPositiveX{m} = ar.model(m).qPositiveX;
ar.model(m).qPositiveX(:) = 0;
else
ar.config.maxsteps = (1+.2*(i-1))*maxsteps;
if(~error_printed)
arFprintf(1, 'Integration error. New attempt (%d/%d) with 20%% increased maxsteps.\n', i, nCVRestart)
error_printed = 1;
end
end
elseif(ar.model(m).condition(c).status<-1)
ar.config.atol = (1+.05*i)*atol;
ar.config.rtol = (1+.05*i)*rtol;
if(~error_printed)
arFprintf(1, 'Integration error. New attempt (%d/%d) with 5%% increased tolerances.\n', i, nCVRestart)
error_printed = 1;
end
else
error( error_id.message );
end
end
end

% heuristic methods to try another time
if i < nCVRestart

increaseMaxsteps = 0;
increaseTols = 0;

for m=1:length(ar.model)

% CV_TOO_MUCH_WORK
if any([ar.model(m).condition(:).status] == -1)

if any(ar.model(m).qPositiveX(:))
% some states have to positive?
% -> temporarily allow negative states
ar.model(m).qPositiveX(:) = 0;
allowNegativeX = 1;

else
% no positivity restraints?
% -> increase max number of steps
increaseMaxsteps = 1;
end

end


% other CVODES error, e.g., -2 = CV_TOO_MUCH_ACC (to much accuracy)
if any([ar.model(m).condition(:).status] < -1)
increaseTols = 1;
end

end
end

% console output
arFprintf(1, 'New Attempt (%d/%d).\n', i+1, nCVRestart)
if allowNegativeX
arFprintf(1, 'Allow negative states.\n')
end
if increaseMaxsteps
% why 20%?
ar.config.maxsteps = (1.2)*ar.config.maxsteps;
arFprintf(1, 'Increase maxsteps by 20%% (now: %0.2e).\n', ar.config.maxsteps)
end
if increaseTols
% why 5%?
ar.config.atol = (1.05)*ar.config.atol;
ar.config.rtol = (1.05)*ar.config.rtol;
arFprintf(1, 'Increase tolerances by 5%% (now: atol=%0.2e, rtol=%0.2e).\n', ...
ar.config.atol, ar.config.rtol)
end

end
end
end

% reset parameters
ar.config.atol = atol;
ar.config.rtol = rtol;
ar.config.maxsteps = maxsteps;

for m=1:length(ar.model)
if(~isempty(qPositiveX{m}))
ar.model(m).qPositiveX = qPositiveX{m};
end
end
for m=1:length(ar.model)
ar.model(m).qPositiveX = qPositiveX{m};
for c=1:length(ar.model(m).condition)
if(isfield(ar.model(m).condition(c), 'xExpSimu'))
if(sum((min(ar.model(m).condition(c).xExpSimu(:,qPositiveX{m}==1),[],1) ./ arRange(ar.model(m).condition(c).xExpSimu(:,qPositiveX{m}==1),1) < -ar.config.rtol) & (min(ar.model(m).condition(c).xExpSimu(:,qPositiveX{m}==1),[],1) < -ar.config.atol)) > 0)
Expand Down Expand Up @@ -218,7 +258,7 @@
exbounds = [g>0; g<0];
qred = sum(onbound(:) & exbounds(:),1)>0;
ar.firstorderopt = norm(g(~qred));
% fprintf('first order optimality criterion %f (%i)\n', ar.firstorderopt, -sum(qred));
% fprintf('first order optimality criterion %f (%i)\n', ar.firstorderopt, -sum(qred));
else
ar.firstorderopt = nan;
end
Expand Down

0 comments on commit 95499dd

Please sign in to comment.