From e288d1cccc480038128cfe1db4dc94f353b5ab5f Mon Sep 17 00:00:00 2001 From: Timo Rachel Date: Mon, 4 Sep 2023 22:00:35 +0200 Subject: [PATCH 01/12] update pleAsssessIdentifibility --- arFramework3/PLE/pleAssessIdentifiability.m | 76 ++++++++++++--------- 1 file changed, 44 insertions(+), 32 deletions(-) diff --git a/arFramework3/PLE/pleAssessIdentifiability.m b/arFramework3/PLE/pleAssessIdentifiability.m index 3fb18ec8..25ad26ea 100644 --- a/arFramework3/PLE/pleAssessIdentifiability.m +++ b/arFramework3/PLE/pleAssessIdentifiability.m @@ -3,7 +3,13 @@ % % Used to set ar.ple.IDstatus_point_inBounds % -function [id,idR,idL] = pleAssessIdentifiability(alpha,parsToCheck) +% INPUTS +% parsToCheck: cell array of parameter indices to check for identifiability [all parameters] +% alpha: confidence level for chi2 test [0.95] +% verbose: if true, print out which parameters are close to the bounds [false] +% + +function [id,idR,idL] = pleAssessIdentifiability(parsToCheck,alpha,verbose) global ar if ~exist("alpha",'var') || isempty(alpha) @@ -12,50 +18,56 @@ if ~exist("parsToCheck",'var') || isempty(parsToCheck) parsToCheck =1:length(ar.ple.p); end +if ~exist("verbose",'var') || isempty(verbose) + verbose = false; +end thresh = icdf('chi2',alpha,1); closeVal = 0.001; idL = NaN(size(ar.ple.chi2s)); idR = NaN(size(ar.ple.chi2s)); -for i=1:length(ar.ple.chi2s) - [chi2min,imin] = min(ar.ple.chi2s{i}); - pi = ar.ple.ps{i}(imin(1),i); - - left = find(ar.ple.ps{i}(:,i)pi); +for i=1:length(ar.p) + if ar.p(i) ~= 0 + [chi2min,imin] = min(ar.ple.chi2s{i}); + pi = ar.ple.ps{i}(imin(1),i); - aboveL = ar.ple.chi2s{i}(left)'>(chi2min+thresh); - aboveR = ar.ple.chi2s{i}(right)'>(chi2min+thresh); - for j=parsToCheck - aboveL = aboveL & (ar.ple.ps{i}(left,j) < (ar.ple.ub(j)-closeVal)); - aboveL = aboveL & (ar.ple.ps{i}(left,j) > (ar.ple.lb(j)+closeVal)); - aboveR = aboveR & (ar.ple.ps{i}(right,j) < (ar.ple.ub(j)-closeVal)); - aboveR = aboveR & (ar.ple.ps{i}(right,j) > (ar.ple.lb(j)+closeVal)); - if(sum(ar.ple.ps{i}(right,j) > (ar.ple.ub(j)-closeVal))) - fprintf('Profile %i: parameter %i is close to upper bound (right)\n',i,j) - end - if(sum(ar.ple.ps{i}(right,j) < (ar.ple.lb(j)+closeVal))) - fprintf('Profile %i: parameter %i is close to lower bound (right)\n',i,j) + left = find(ar.ple.ps{i}(:,i)pi); + + aboveL = ar.ple.chi2s{i}(left)'>(chi2min+thresh); + aboveR = ar.ple.chi2s{i}(right)'>(chi2min+thresh); + for j=parsToCheck{i} + aboveL = aboveL & (ar.ple.ps{i}(left,j) < (ar.ple.ub(j)-closeVal)); + aboveL = aboveL & (ar.ple.ps{i}(left,j) > (ar.ple.lb(j)+closeVal)); + aboveR = aboveR & (ar.ple.ps{i}(right,j) < (ar.ple.ub(j)-closeVal)); + aboveR = aboveR & (ar.ple.ps{i}(right,j) > (ar.ple.lb(j)+closeVal)); + if ~verbose + if(sum(ar.ple.ps{i}(right,j) > (ar.ple.ub(j)-closeVal))) + fprintf('Profile %i: parameter %i is close to upper bound (right)\n',i,j) + end + if(sum(ar.ple.ps{i}(right,j) < (ar.ple.lb(j)+closeVal))) + fprintf('Profile %i: parameter %i is close to lower bound (right)\n',i,j) + end + if(sum(ar.ple.ps{i}(left,j) > (ar.ple.ub(j)-closeVal))) + fprintf('Profile %i: parameter %i is close to upper bound (left)\n',i,j) + end + if(sum(ar.ple.ps{i}(left,j) < (ar.ple.lb(j)+closeVal))) + fprintf('Profile %i: parameter %i is close to lower bound (left)\n',i,j) + end + end end - if(sum(ar.ple.ps{i}(left,j) > (ar.ple.ub(j)-closeVal))) - fprintf('Profile %i: parameter %i is close to upper bound (left)\n',i,j) - end - if(sum(ar.ple.ps{i}(left,j) < (ar.ple.lb(j)+closeVal))) - fprintf('Profile %i: parameter %i is close to lower bound (left)\n',i,j) - end - % if(sum(aboveR)==0) - % i - % j - % end + idL(i) = sum(aboveL)>0; + idR(i) = sum(aboveR)>0; + else + idL(i) = 0; + idR(i) = 0; end - idL(i) = sum(aboveL)>0; - idR(i) = sum(aboveR)>0; - end id = idR & idL; ar.ple.IDstatus_point_inBounds_alpha = alpha; ar.ple.IDstatus_point_inBounds = id; + \ No newline at end of file From 4173e95b2631fa6c619cbe76569e32bcec13fb4d Mon Sep 17 00:00:00 2001 From: nneubr Date: Tue, 14 Nov 2023 13:21:46 +0100 Subject: [PATCH 02/12] console output easier to understand --- arFramework3/Subfunctions/arCalcMerit.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/arFramework3/Subfunctions/arCalcMerit.m b/arFramework3/Subfunctions/arCalcMerit.m index 74dcac41..2dc69006 100644 --- a/arFramework3/Subfunctions/arCalcMerit.m +++ b/arFramework3/Subfunctions/arCalcMerit.m @@ -161,7 +161,7 @@ ar.config.atol = (1+.05*i)*atol; ar.config.rtol = (1+.05*i)*rtol; if(~error_printed) - arFprintf(1, 'Integration error, restarting %d / %d with 5%% increased precision.\n',i,nCVRestart) + arFprintf(1, 'Integration error, restarting %d / %d with 5%% increased tolerances.\n',i,nCVRestart) error_printed = 1; end else From a5b5acffea9092a1d22257b3f48f15559d1b0f10 Mon Sep 17 00:00:00 2001 From: nneubr Date: Tue, 14 Nov 2023 14:26:09 +0100 Subject: [PATCH 03/12] Improved concole output --- arFramework3/Subfunctions/arCalcMerit.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/arFramework3/Subfunctions/arCalcMerit.m b/arFramework3/Subfunctions/arCalcMerit.m index 2dc69006..acb1cc5c 100644 --- a/arFramework3/Subfunctions/arCalcMerit.m +++ b/arFramework3/Subfunctions/arCalcMerit.m @@ -153,7 +153,7 @@ else ar.config.maxsteps = (1+.2*(i-1))*maxsteps; if(~error_printed) - arFprintf(1, 'Integration error, restarting %d / %d with 20%% increased maxsteps.\n',i-1,nCVRestart) + arFprintf(1, 'Integration error. New attempt (%d/%d) with 20%% increased maxsteps.\n', i, nCVRestart) error_printed = 1; end end @@ -161,7 +161,7 @@ ar.config.atol = (1+.05*i)*atol; ar.config.rtol = (1+.05*i)*rtol; if(~error_printed) - arFprintf(1, 'Integration error, restarting %d / %d with 5%% increased tolerances.\n',i,nCVRestart) + arFprintf(1, 'Integration error. New attempt (%d/%d) with 5%% increased tolerances.\n', i, nCVRestart) error_printed = 1; end else From 95499ddafa78e5af661de82c5373581f50cfce58 Mon Sep 17 00:00:00 2001 From: nneubr Date: Wed, 15 Nov 2023 14:23:17 +0100 Subject: [PATCH 04/12] improved logical structure and readability of nCVRestart loop --- arFramework3/Subfunctions/arCalcMerit.m | 172 +++++++++++++++--------- 1 file changed, 106 insertions(+), 66 deletions(-) diff --git a/arFramework3/Subfunctions/arCalcMerit.m b/arFramework3/Subfunctions/arCalcMerit.m index acb1cc5c..6060f168 100644 --- a/arFramework3/Subfunctions/arCalcMerit.m +++ b/arFramework3/Subfunctions/arCalcMerit.m @@ -1,77 +1,77 @@ % [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; @@ -79,13 +79,13 @@ 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; @@ -106,7 +106,7 @@ end if(~isfield(ar, 'fevals')) - ar.fevals = 0; + ar.fevals = 0; end ar.fevals = ar.fevals + 1; @@ -114,7 +114,8 @@ 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; @@ -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) @@ -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 From 0f83333cdfb044255e95f7b8b935ff559a889098 Mon Sep 17 00:00:00 2001 From: ckreutz Date: Wed, 6 Dec 2023 14:29:52 +0100 Subject: [PATCH 05/12] Additional example for RTF, dose-dependency --- .../Data/LDH_WT.def | 16 ++ .../Data/LDH_WT.xlsx | Bin 0 -> 19109 bytes .../Models/TransientDoseFunction.def | 45 +++++ .../Setup.m | 154 ++++++++++++++++++ 4 files changed, 215 insertions(+) create mode 100644 arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Data/LDH_WT.def create mode 100644 arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Data/LDH_WT.xlsx create mode 100644 arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Models/TransientDoseFunction.def create mode 100644 arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Setup.m diff --git a/arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Data/LDH_WT.def b/arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Data/LDH_WT.def new file mode 100644 index 00000000..cf98399d --- /dev/null +++ b/arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Data/LDH_WT.def @@ -0,0 +1,16 @@ +DESCRIPTION + +PREDICTOR +time T min time 0 10 + +INPUTS + +OBSERVABLES + +ERRORS + +CONDITIONS +isKO "0" + +RANDOM + diff --git a/arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Data/LDH_WT.xlsx b/arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Data/LDH_WT.xlsx new file mode 100644 index 0000000000000000000000000000000000000000..f78efa45703dded659fcf44457668f39673494ba GIT binary patch literal 19109 zcmeIagL`Gowmlr9V>{{CwmTg=opfy5cE`4D+qP}nR>%H!zvte2&Uw%M{(|3KE6-j# zsa2yi=d4kAYLuKLC>Sab1P~Mu5D+0yCXI!)IWQ0qA~+Bb3J?^Cx}dd{gQ1m!_D@$E zLwij+7fXw;xnLj^*+3wG_y6zuf4B$A<9lTK7+?dh#M%WDsLkqZzd=>3>2@QLG2(wE zxx}WG&gH;aTiY|z|6W}F{jA^z&g*AavlCgJrhY9g*!*Oog+{mvgc4W5MXABn*wI_2 zA!dMLF9!_i8ffjx-HTc6x_ftcEp1|XWUj^QuMx_{4 z%%j--dt8bH+ishB7_Sl}Q0jdx2)}zE!RcOS>}1iE^QWt090_w#XeydQkB&$SqESpk zB7QOdwVX^;AxvEVe#p0%?`l6d8!`04Y)c7lZ-$sJ5M9xE;T58r$X8h%{Zh(`P_=`0}DM)=d~m-sqDLqr}4&PEqGF5*vYKTuKc z4zIV-wGEz#(_zAgU6%4NWK?eAM(2vaq!(L9NGdYBL=oHaoqi;zh5Lnv6j2FR3a5@J zsOY+%7Bc=E43NSP3Iv1= z1O@D3LH`dwakjQI*R!@Z|1*{Sx6c3r<}<*%|8HOAv14Za46wqF0Uv%dowiJ#a)Gu+ zGxEunXI+6N-JM}&8Rb}*MA`!Ee6{|CyA40GnWvB#gpF`EZMH@v#zbp%u&B>Ru zA(76m$&3QljcRLoXnd-U(fYGzxNiA~A29U)X+txnjmBbtAX@+i1cU&117Jh{YCYvY zHLd3uP(4L=`Bl8Ge$UjO<&w=cLNG%2uq0drlmF2SF%$u+#uYk{ zPBLas;>a>$6z(eoDF@Y-Qf>{gT1FRGB}AQhnt0_|;8(WqU!2bcILmIIoXYBbBJvv% zdMtzrudy6x?eZ$zPebALX1?NsN)p$G+@Ti#??m zx|=Lpa$?g&SuH4V=3uCS2&m>>-rZ(+zuo7T9LCllHy^5!&}$e_3t1& zd?r-kzZTBs!?z7hRWOJZ9-Bxe-39OjMs$EG*DMbE8{a56{c4TD^nA${lMACJ8ACPb z+?)G?%ad^pH#P7EGe~%)_E`Lw(@FRd7v!+g@+`2u`Jq>zR@;{+7qCg5emQtUviUd` zm20Ee+gCy`tBAF@nY*k+w-a_NJb%RJbNMby`4Us2^~h==Cjon0>e#f`h?2yN0`q0u zaz^qf0@6O3lX3&@ebChfNs5zrC|PSk>Br#*6#hE**(_;9f2-Kj(%0WrlkZ_ByT)@X zBx)M|al$Ncj-RBf4_uU#VKH+we~$D@_q%?&;tqV(guiAy>oPl7=|$xWG+ zE+UF%ul`q~1GdXA8UUp50|FB2Ka9@aMAy#HK*7Pz)XLcYPuy~d7BK5$Km@qzy1+PJ zOx!nx(0rON8xBZNC&2!^4Z`IyrC3Q14~`54+azO;_}=LH4_h_L+A16*wbfdR+0m!- zU}m)KSz1Z5*(yA&owEMuUeUT&=m>$T^Ub7<(;@rD6)l$xZdwK5EQwJ+bCZxH%Q!La z;GBbz-8)K3Me#hp$~j_-(VTySL*;|?c-A-O;2@UyE|fQb<+(G$8d9~6#7W+IU{5fq zwRlS`n%;Y;wLB@f-PG4Lw+FiZzmppOBS(T*2q2(9>^}>SzkJQX#L&`^{;yBQKdnAj zl@7yU$Lz#-AqZ^Y{z%1IlWgz7UlUK;ATh{@$2C-cpeVyiiD}0{0U_nQm#4@tNcwKg z8$bIE6t3fLh=!#8BnQTNQ&MtUk`w>OG0`d#Qe+M2JxZtFsufAUN0G*% zK845|NZKA)r?e_IF>uvMIgj`Z<8p*)i&6>L0cK3UBQmPr%Gn$xmcFyIVWVu16Hzk_ z=G@&k?18l95nV;1r`RAC$4ErWw&`AfJcJwG9WGQ`V1iXE;rEt!B=J-}!k}Hx>bz@C0TWX*i^LH7#5`b1H!v zMxn2Dh*W8@UZ=jlKpaKVdqZgXFuiU4~=+4z-iZGlo)9eM#x9@f*v2ZJ!E27O#@1}F-N!ImxPd2KX!K;J+pkx6X( z?IMdzV4rb_Bt$4}ABJgs`_uCLZ7C|}sr)4HKE|?&GUG+3;YHLeM_`h#BZx+ZcJ|`O zh1cFu%UiMcpX-dt5wbFQ-tPBj^5TN{k1VMbL^`M`Qcmu*KCTzEw>&T3p2r^P3A{5{ zJ+Y@8-@>vwU2fK{s+bUX=`?QM(#YvE`CLCPlDyxZmi0Q~a%12LVlnNHu=$*7drWlC^@E_w+wX2 zEYe2Pz_nd!A$fak6El7uQ3QkMjFLmFzD6h!4)lm&AFZs%j3u(1?dq2+l#L}*ZL`nr z6kJo_>D)+#$lxo1I@&fR}Pq~|yi9Ud-3MWso|*l%b`43)#?QRlsO!en;lh(M4c zhKU}mHjKc?y6%ES3l?fg7DQ3l+V>b8Y(avMc%m|N-||ariosCE+n2Lxp(XRXcSE)| zctA~%XBM44O3D{Ile(55%()m|CrlBFOm)sUlG&Q9YeBIe`|CQr5L?EFP%&4Cknt99 zkCwg>Hd(Y_P2V&eC4_y4>BHuc$HtOV_H4zApr{joFW5)TlkVWh;vN=^0DV~~{C+{Rgwx^*7WK^5(>)=df676Ll7)nXCu9N- z^H#cAJLzKOHkiS=SAD5GU`8D-yQ0(^h*G)TbqAaGd?) zyl3%nq0*`FVU!xA_A)Jc8k<>so|Xi=?0ax;nOJdV225CDfUrWUc7Pe$Km+x!RW7P( z8#|qr1t(8>ZpdG>Xt^##d+~VC<+;>;_LVW4hf<}+K;=YASaL$Gq%Bq#ndh2`W>rY0 ze*TG8WYDy5#4J6oPN$IMM_6Y{Nb-J}4g+Q^s4UnMk4WZ}L`C&JWUcbk{hU^90vkaJ zHH@iE8V_f;l{$r@K>U#BRkbd{c3BP6n}fb@DI zN=IfOY38#d-9Zo>|+>SRyqdhyD>5UpH^UM9@j87ouijw2YJJ!qW zNeHeyDQ;>&sBdbY+Mbrih%~bE^*mC2JTFmmlR=r#n`HSZJelgm=;njK+4b{1qF;SERyZJev&5dbE1}ynR03DG~-`0;LDNtKtZ%jC|AS zFvzW_PDuu_%&|bi!|@GiXZ zQatuMfgY)}+PndaivOgGftyG>Er1j;I3OUv&mX#I@8D`-_}7#@ue@fp$A;)hKkbe9 zRC%m7W{#7HVT%h4V-W@(3*uQoB0-jj;gc`#8usBH6mLL1VYDQUQy<_XhHD#mG#iri`%O>;o2bYoygw#_ z6jVm8AL9fUl!NV8+%xJGL||>I5)QWNpN`o$!Zsxp$i@P=fnmS&spPD6C&DDdZOZwY zBDBRebp-84i#fT}hxZ0v1ZneoIYEmxS%w0ZgkjL4$T?=>1UYEq=WCgiby%Vmg<*c9 zo>8fQH*=_i)K%4?k_uox_BW+Te97zh-q@RO0>4Omd~ZzXjgs;X#whS!c_HF`gU3pCDt@#P|l>^*bRaeJ0{Ak~4NqJJ< zF*BDYttD&0EeupS3Qn6kbVb>f$Cszqv!szZ*N4OK95pcIo&%~8XP7}SepqJ*?&@&} z;*B;N)2*w;b8AJD+9P~--wvD>i$piia#c6VX6(vx0#}0+-#E}SdUmVA6^5F zz3*^ALp?2Wz`vrJ>F1BdJHdxGqdvK=hjibkx}Uubx}z`Uc))efCkvme7D{~e#2zYt z3ym>b@x@d6(wWka-=LLRyIhnNXn*^bsSx_jTh4dWK?#J@3ft+e<$&`?ih0H*3wei8 zJ$m}G3@}joaa?R5>LBaYv?i>=8*dtPjB^#yuH%9;7H~uil1^IdPW}dQqvpgznHOu0 z!7FZkMP`n3bf}k4MdQ}WYrpMT>5(RvJLNBH|D%D%4^yV8gCINu4H}Wpi*?1%aDMHW zU)^t(2-U=sVVGcSIfMb^3D3oJvd@<*rX<2k6_vI!WJCM3hJjzwdTseRBS`s**IS;I z67nLYx9cXgN5A}v-*`5&(T>35*eI~|g?eR|AEcmj*d6^567Wt23<2r`lM=7U@M zAJG4skN0|&ID=CH0daHz0b%~b$Nww|?M(~~9qj4-AKT(PB3nK0Wba`3h$+e^F0ph%-mvU*n~VbZx+G z-&=ovJkOP7$wtXOB4}@>F9e;QFPw~R^u4~)^LKE*UVU?o+Ji{I0wlgIPLKy7Jb zh5A<7%fscV%=DEzpGxbtzMiqFrNh1R^Lc+it@EJoBg@^9S@Z4cW^{e;{%Gg@_f{wU z(1G8_(OwyWb(TyLLe|I8$NQ})SAVyw_lwPIlkAkY+tb08%pb5`eO>WrdxGXJ%*DWFZ0_3u48cZcWD zuwGv39ekfDU5Asr?#8IIRut8wo;_6D!XI}$NC{# z&tBIjiTrO{_Yt4d>|R$1FJ#`;>~6Pc=y=}s?16?n*NnXwrVk!5gI-~Rc_04HB4=rJ zPG<_5=$V=oQ+)03JM>H5FFS`X+LdWxX@}a{S>6m=o*Z|Y?C_457N?*Z>Tr$veY&Tu zt=8oCG4WE02%ZaW5ijoIJ0A?P3hXj-=ar678Ip1nvNvr(ErJ)>jI0v&$$x=^6v5-%{8VJw%YcVrioIn6BIUP95p_=nA&(3*nrXTcU?g?*%K@I zna3Y?)AF>bWkK7tg{4jUapRYko)p^JEES*6HZ17$;33R0DD`o?0vC!_qT<5~@*?Q- zFTrp8P)${7wl1y7EL(}VZNNcle#f_Y1z+dFY7sGoU+=9eyH5dX>S{h&is7F0WM{)) z7w|J@^^f6}epdls%jt;-I4Ute9r+|Ng0rGoU%HqcS1=PKvI*kdhKFfBEqb0T9v4t| zacba0nTE$;@B3oOtvtuW8V9QU?od?$0b}W1{9?AxeS9MKMLMV@bG7g$P+#e3s~*b)+6yuKCr>?Bg>qZA9w`}@YV}V=ze24flm`vsFB$jkKefI zr%QLP|D|D*80>XHiJa*x=I5`Q-?`1My3MYqBSkFA{eto}s#KaY)o*8XBF74H5_D6Y)^{7|OG@sK)NBexwRE-CaxYpgt`w3`lEM(1DQUoTczFL2;9A zikb}C>TTY}MWo(FO0-rS4!*Ii7j)ups6a$6@Zy~U!)WeNFPl5Y=|w-@G;@jKn2pJo zL#88gJj>m5QP0TlYWUe-*k&To^3ml@q-=?-X&3eDGe%i%v3a#CO1u>3+vXhU9@}MF zk!UIhIh_w$UC)-agxYQnOq$6~)z^*rEvS$%Bt4~=Q0fNeEcoOu(0U-^`sNj6FctMb zL>16hZZh9n|u z4XfcTt7#Hs*Gu!imdLmkF3E)WHSyI;Ik*gS$q|-d>>(WiSFDoGwZ62Nol(NBWitml zVh9q`7)aTS5=pgFNZF1Q)5fsO$oGDcE(4X-$9=G8Z<>Ta@bleWrG(CF97?`RXo3$x zO588IM37#9o7?u}B$a0OM=B&&T3(StndcUOVYMnjb}AtR zp)+$ntkh7OqE+8;LQNDSc{x-G3DN##d!0VV>(UYoKGMZA?;it;`5`#bB*SHV8!CtZ{2$^J6PL zne!7XNyER})q>3;>z1)&slyqTff}3Rkfpm`u1&%zO9y2jb8DB)I(VG6MYh1-0@}f6DUiYgv-j3mu=1|qI`6$ zb&=l<)*>r&4J0GP(5fJIAr-XmHQL{NVb#q+)S2P2#D3W$*EmtEv4aI!6sNdz@n*Cr zSIF=ETq2Y&8>&pa?6xs5HJCm~8SLl6}7$Wh=I)S`aZkh>Uc#gMNd zjb1-bo2j3dvJVuot1`$>3l(6aM%)E_ ze&BvPrEhi=8uSo4qF|Qpo7}C+T;$DUTXXlIBr94W6%h_7McLGNz+AL=QQ9L@omhL8 zZi`|Ni8>m&#<>Okm~M+;Zi`Nkv5?ZdlbqqC0cIny6SQ+$p?|dD;GGw z6!q^;vxiIR9m5*BpS;P^L|DZN)A|UxlVVDyipWxpE3&$|W$d6%;HQf0RnQFNGIJLt zn53FeG6oj1`V`)W@0FB3C#>3*n|oDeAO|tRj?0x&rp)C^rDJu`9uA4ZB3mS|Y`F{q zk%A-KJZe$b+#mNME3i@&rC2HaamI+(xdIQHn9p zDe+q+@u$8FeOuxps-9n!61+I~kGyPMYAcBV!+rBbH73(Q7k6%677SQB$bh?m7@Fn4 zA;~fNVSuFH7S!30*l8Nta{A>7Gkm@$VLU>|PJu4k+(9J1%?jC-FkG*6y$6EnV+y zgi5ch=R>-HjipblI4`a?(J|_Z8Lp0;_*0T(GOX1Q{l__879n2Pkn0%(UgDY8Bry&-NYVFV=w~Nk zCtyuEL}o#I(lT3X7Uwqm{i5$8wjgxL;;re|Is!3QD z+*Xs_My6MhuLoO4_7*H4ct_TFQpb5xpP(o36EBc1Yz7#Q*1k(O=%CumwhBmPzLCfy zEWIhSVuAzlAvSTm+jzjQHeRX7!8Kg{nN*zC89`o$0YM)t$hb3*xV|t3Eu;Xb6oY7m zT^<&@1EFxrrj^XL$QqWN>`rrvuPgA7GjQT@5f!CUyZIIY|G?1z;0U+$?+ZzWFiBOy zNJVz5tVMU`l$*A=Ww1F(r2CkZB=6>JyIDXCM#ea>VT<7UqqVHbuuG)_8d;*A8vhwLDpP zyMz%?VA>nwms1Tbe1_m?9oZ*%p8EaiW^L(h=&%>k&?+Iyw#D4C)e4Goi7=fz z;J?ORiYgiW#G@?p_Y##6Y6=>haZ$CzCP&($<*Msz6511zgl1%^gkiBt&&E~|gk9DQ zF^!HlFD6Vq^Lb!hsrt!E8?`CxCJ!|T@FuHk!NPVY^B|}J%^zCQ3(sjr%b{MJEOL3Hxqbt~#Yw6+Sg3`o* z(&UYaHU@#|K3T3-_)X0jl$dD?kw%FkUWgNGU)@+(22s0?rLBjMcV-tJ2^qksq34iW zfZNB=YBe$ZV}Kh0j6SXghjb-)CmfWJC<}j)lJl{&Tgr`w=OCK166vIHv4Qk}%V6BP zD`;fKB~713tpj*bpUth9*3X4P?n2k@U5#5jv4SXrvDI;)BNCFEv?nGzH+^5!NcW5} zmMn9zRvzWoLTR>2C@QZ4;Xr^P#L|_n4a?>(de|gQBk2c#pIPi!dv03`4XRW^Xz&SF z-Db#y^fhYk3eaBGcH9j|6xg0zw-i?QN{i3I=^_N@{ooXDXm|d~q^((V2N&r`Pl6jX;oZC*-V`))pEzQ9jz~J924JrvId*>Jj|?q*3@tN4 zhN6j8Pd5Hq=ezM7g1Wy2P_F|A{J31O#G^KRK+YYHGF0hFgA1U`b(MqVlWxcQ@QpXx%aNxL5+1)fSx37Tn>Z>@!}Ewfo7i zzq{W6=KQj2WW6odV!3us0S*;#35?ym;44}Ui&WZs&-SVeSKn12?_1=c`Z^#4;SG^` zkposTJI;dJ&csS{=RXEz-C7)gUmFyVY3%r`&lxO|^OZhR5C~?nk;pg~H>?%%0)FsWD)rh1{fC?GtI{O(o3oR}? ze~p`Kcam64fk#~yCyVQe$`_Yq^!q*HA*Aw(Qe)*J5fPHyahV>l+;<5vcFC3uDXWyM zC`z9f6hf$_E|E;=BKXyX`K{8B&{ie3nA!0X4*f92um;cF4dO6ogN-zzXyST(S$`4E zFr_4IXBtiWwE)jFh!-qsaINzucOfkV>_EO!raqn`4UYn}EANSb0G7W8atfK55G!pZ z*RQstPzF|*WT3JHR-%r5$SD{w57l5L5+p~&6iS&`F83Xoqzfp@H-ar7XQ{u$F+>9M zv6q){N=T>#_C|VDDycw*WFQ|2Vg5|&b%9KB3e-3K&NUBKlj_>05nweePvule*yTT^ zGP#f7?t&N5G>03uc)m<~!!;uxsp+!u)fKNtvJ^MtIgRDOMq(XxC{cDe_2Q)QjI^Ew z#-!hCf@=dwQqZ^%+nH=5oGkPF0kN$InOjq0;KCP;ZUz9N2U~#yQjxce&9vYca+APW zAM%D917ijm*Rtu_6T|Wl-;S6@e85vLhF#_7%NBRLgl2pV<tN`UT>Z_Y|bF)oaqDGlV)S@JseD!_^~GjjNL|HgdW(&t`a_he3VF3f~Y?~ zr$HjUz)k0`TGsXcejI3k4VDueLFe8;56Zb$7U8~(8_|>upbXrRk`{mg+rwmImtY%` z8A45Jw#yO;3eJet2`t$lme{XLO=U;j@O}x&8@Dgk!a0fx3z0p?#XfGR>am7*kFpb%2X4gD#+!(z^ zeLG<@!}{OOo)oF5@z z=-d1zEgZxO+E+hG5}DjxA}okfZ|CakR}=yrKAdc)ngIdLWerBF-mNb>TTNmHuuhon z-44)a1p6@xZbgx#H7 z-HO!Nglv> zVfO6}cJGnVkg~F=Xuk$6jO>$*eHE=!MjM$$M3*hEjxPlAU3jcUu?|2HRS}Q-+zyXm zu|txzb9Y%caZQ&ow`u!j`xjT#cUf;nSCfALi%RR2RZEftv`wwr9wBB^L%K21D2GqW zrd4B1WKPc%-C37znVmTVUKmr>9SY!GP#d|-#?-k%t6_);;DadsR%#g26Bf$7p^}S` zE0ROA>BypNDq{Va`2BE=Cw+{GWVDJgz;N8YX6i`d3i%`TV}VzuKj{HW>6|p090i3w zZmXjq+dqLiZdw|#)s)mBnyyZ4&Z1bilCMW4md@tuH<=Li;{6}^X=hRPX*b`N^Bh;j zrq*l$-A0jr9F#BK%JL61iz$|sZ71M8*1CGFH>9gOQC@4l348T(EC)-g9yB6T) z0n4(zLc8Y$c}X=aIYfqv&abJAxJiD(VP$p*eHYupme>>RT_H^A4N9HIEVu?4onKryve8#9Vm{>Uc=v=O$G zOZ_5+5MHZ&CnoI|`1EQWj55=-h9>}H2+na($acu5!mfj%?ecU#!ImoOk3rC9n9x;e zse`CzIf*cXDru7;51+)ro0%D_c-TtI*Zi?G=BNop!t5jI_3Gyv8`Uq9kZS=;jgtU#P)|+ z8;Cy$PS?VR#9}%9F7YMURV9i>Aevy76QKZ;SxvySG{z;jKcejCkL zP-IG#Ta!%K{mA_;xcpm*JVFC9qfJ*S7UiagvPpUJ<&w}r;V_zGWCi(5j>v&xy%;=w zT)_5B+|~$1T`V5DIHgQHokhHL8LPZ6i`X|E`mI@6ZJw(xg!ri4zdE&2A*|q;Xuhk# zev<2)&~A<`X3t6YEio8Ci-zchtkj}H_o{gws4&HsmL*X*0sRaM-=i0UI2ZMRZ4eoe z!5_SdBkTGKbM2;!B4q=$*T2?e+1nn}JmW&`2OuKPUuFa^LI81l4r`aKDuv^NXb(KV zlgyE9k`fx8{*vy@idpm&kJwbnw1;BY`LMK5Z=>0u|Ky36B%gr@?G}GmMBlH53=P1u z2MZkS|17KlF#Qp$#&I&vdFWm_zCAMH-I&@KZs|rO6+?)U%-X5DSd=!c_>8MQ; zrw#jPVQOeXYG*4Gr|l}G?YVCa!WpS+VR>~TQ_+{7EnrYp7)qF4f*z1l>=~^3k~KkX z6P%ZWe#q^lzsXW1TnMdiTXRSSBS5GW-rF@YA_4#9pH!({hhLX^j%Z^UkBwyX)YXuN zC5mf|a=ieO|uVjz0C%*4(SLfMAlkfYp=}cqgH`mw0|SInMy}KDpLCS zQ=x&+uyip>EpQP=m-ew*(9e znZ1f+(jr}An~KA%l?$+**2iq?uuwTTyK_LLjU6n}X(+Ctu- zNJ(po3g2@wB}rONH=z;;oR~Pv0TXgvHtat4GN-?Q1GET-&HW9o`Af?COk_>HYDGZ8 z76;GJAfcYA?+_YfIYbphU98S#-=Rr+!lH`9p{QNRB%oOa$LsBvrlyU1`UtJnBxG14 zNE!MoNR!LODYfnj+ob2CepzX2C3J~wO1}2Mf=K3tVtRjJjPQ`?O5%keXQ;Qvfz2A( zUq7ZM5##nX#wn_drN~7EyQSDsB=#EEi(Dp5+Sn_hIwio`2QCWF?lWhz5aHIW#t-Lw z)f|w0aew4usAw7vT1ds~yV1YE7%S34wO_*p?wsXm9LX#;hXPMpExUG07iS|MmA4HM zG@|dg6bli|9Wk%oJ9Y{PI0;FAE9$0J$tq*$0k%Ck_R%*Y&f{Kwsr6iRYCVxV6DP!RvEW3@V!qjXwI-@;zblobH!}M zpUgvXSE0y%MC3y^FebV=6h}Ii@04MS3~c0tZ~szWg*74Kcz%8KYjJ{MZ@b&?gy_#HXSwieR2t&73rC{A!lnP*pA1T=M z6-3^UF?wX-oQ3unC`j(OF+)0LZ}39m>X^fD{fsO296OUP^+DLj?uP?gqEf{3?IZ>& zZ!$PlJZ98>wF#kGCb5Y$acZ`!nR8AcWb$(L1`4DOy6{fA57y+d`SF@HYDVNCvm4W= z?NQmQYf;)hb*{iY7kK8J+|V^M8;8^XsS!oq#Vo}7VfKxgVFW5}>TaWriL#lLfm+m_ zAq7j|n(J9bON)2KJJ)${#F|O&Qwop*&19efqUd6-QN%F*pxCxq++y;pa8bSxw&C7z za5D2ynWHDZFvzlbXw}c+ttJF7wgEw0x1_ZVMs(v_APEF)w#=zs>Jy4GoM2T2>F5Aw z1XF|5$?23}wn9F!5O(YW98M)G92ezNF*2O6u^t?FsEm;ZH{9|3-z$}Vcw3DFjM%Ye z@`~S1cp66dq)G*)B-Yd6O%uPl74e|Z7jU^n&15^0Rp_BXw&heehZB%H7@a7Zv6~5g zJ-aXGhG0>DWWF#yWqMNw$Y?fRO_B%ew;2 zn!5caly<^S)U$urhT7?tcqX>Q$7k>gglh~3qXRrUto;W|0&qEVscf;?+H2CY_E|=|gWlcaX zF&uOms|`#JRMKE~Jb2-KfcbeQJrc3?mEg zK^gWTpLHXv(!)XO2Hv{YlK?9g4b7nMRm(kRG*8u_nltu< z+ck`t7!2>^Hs~;g?#r=6%CYTNYoUY~*;3(Hy&974u8$5J>*eLeO6Lj)y{X;kT1WpH2FQ?!OGl#+oE5 zo$wSuR*=*r3CXnuYa$GB{2H-7tq@VJ#ZtRh+V*F0iX0dLGq{YKzQv zHcD!!XIn?<#B6QJQz`X8H>+Yi#c$l>E2J?npu;GtPq zv0DPx&RAaqU2nUCL*nYm8Ygq}!pFfR$XKchT^9AkH>gLEF|k@q7n0E`176vO4L8X> zIUtp$NfWcksgZZC#}u8k^y11q;K<#!jG4GRoQG*cHmu5>nj9_khA6kA%o-cd*OEdTpjwD}!_8r# zdngo#Oga=8zcxWl|rh%ts`oN7r>Y;HmgBLpQZ;~@}Mi`3u_iXD(ic+GE zF;sM#5_|X&IiVA&ax4vE%`S*=?gb(wqll@!NRtA0;9o{2>x&a(h!Z{n!u9wPh)VQ?^*#!%hxd#_^-Dqif&4N#DJ0p#j#fE1}>+M~h4eE#W2hmZgeJ32sbpMhxJP{uE?LbP$K0ctgz&Jlq>R z4_%)+K5p+UO{B~}x_HBWR~2o&pua!R3`@>fdxm+Psdvd*zY})c z4;?GLK{in@JZJK|n0_2eee_%=zR8VLm@C8s9XA1<1GEcQ3+$hjUg}rur%ldn!d>^YFBcm|e-#s^Ub-tt0U#=v?*gQ&&yp|_)dK*cWV<1@h^ntf|rg~a6o;z zdjR#d@BFKpo?uFUy|kNn`p7lqrQ9_a(s=Kq^pDO^y-(iE-C7B0++d%fn^oPH&K2Ckf#cRYPWxX{d|7X|ZM3N?#cB}f2uDRz~ zaq%FU7U8bCTf-*nPK{LR_CJ3cvE2peKxB?2DGRPQ*|nWcJWi~HU5aiC;zD0CCX?hx zA70;1_D=3R9A78hhDSzUcCAFJyELcy@ zET^>IpRV8Dt~qb7_f8((PiA=br(UkFug^dA-LA*__Oz4u_b)zS02?F!b9;1ssH>p{ zK-6RjP}Cs-6qyFr`f_&GHum)T)^>(}^gRIa^8X1j0bS=GyCu=bfH`mlzLm4zDZM(u zDY)hnw{iTfn(2m)WUJ9EEs*5($nP1XGugCl0Un$_&c0K7@>z^$B`Kf@-14^s83ioV z?{B{)_U#g~k1P(OVvK3xe9;p8!SK4X#ob(y?1fzmc77;-sh{X^i6ntLpb{`QM;lna zW{n>QnniZi}7_9YnTxRMZ^^ZCKWTc z!SsI`G8afXobnMl;T&7Xk3Vw=Uo(c;49#%2{yo{ecHz5`9kTN~=*jTodj99Vz3Npc zsDyqG22jamlyvU88 zrUXDcEjINPcJkY<}_7AC7FVeQW}G*m<@*E?Z)v zPcpTL=P^}Gs7aGax+(a<-|77t|Ffa|!o}>>lqTXu5&hod?P+Ja!SnH8U2D)x5;I98 zW9j~-r84PKBhbk;wbF3i>()Ls$Z5kXUHkNaE9s#=T~alulT9tiZ5d?#W?IocuM_*n zf}#3%#rLeG7thRzu4hZ4+KcDMnYI}(z3cFiSL1W1%}-9cOLD$*@}4X5nG5m}J;8Uc zt&Q1rDC@p^^dQxK4%|aqm2O?6AX3~!0^En}ubZQM(X_*C?f4_``0)7n_`vu`_|W(% z+~&Se0#={vO}}FRslkn6P5`&RzJ|IEz81S4UROo&r~P6s)cAGn>uPUbPhTH#H{k9i z?jg1}mPTx(C`VfHweT~6L4XA8I0lM>0^DFx;~A_1%emDe1~9d7xwJgEJUl!gJVZRm zfSbvA5llC%7yBQDUEBP#sQ)R<<9Z~*#?Jx<5fTw08Ymhn8Z3G|if*C`Cza&gmNi4%G)` zvYi@OfB+*@I)rL{G4R-m!u&9qR1BCg%`&g7^=w!|td6=&+o_e-+ff=JC4N-J1Xr+0 zJN2?hVjXPyLxb&_ujOb`8$FGne5KMNH^nk^H1B8>K|j?CPzogp@Xy6{MZ`eoXniOa*wv@fw&))4}3Bv@tO3!MhoN^$e?vA3{!l{=pYwJag9 zz^1vKzn68|=7{yq1<$J&491p+z}`j^80EAq-of&*;h&%rh@K%D@l KI3xP!xBm}g@sK0{ literal 0 HcmV?d00001 diff --git a/arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Models/TransientDoseFunction.def b/arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Models/TransientDoseFunction.def new file mode 100644 index 00000000..cc0deb15 --- /dev/null +++ b/arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Models/TransientDoseFunction.def @@ -0,0 +1,45 @@ +DESCRIPTION + +PREDICTOR +t T h time 0 10 + +COMPARTMENTS + +STATES + +INPUTS +Signal_sus C au conc "signum_TF * amp_sus*(1-exp(-(log10(10^(10*t/maxt)+10^(toffset_TF_sus))-log10(1+10^(toffset_TF_sus)))/timescale_sus))" + +REACTIONS + +DERIVED + + +OBSERVABLES +Response C au conc 0 0 "offset_TF + Signal_sus" + +ERRORS +// the following error model assumes the same error for each expriment: +Response "sd_TF" + + +SUBSTITUTIONS + + +CONDITIONS +timescale_sus "timescale_sus*timeUnitFactor/maxt" +toffset_TF_sus "toffset_TF_sus*timeUnitFactor/maxt" +signum_TF "1" +timeUnitFactor "100" +maxt "10" +amp_sus "(Max_amp_sus*dose^(hill_amp_sus))/(KD_amp_sus^(hill_amp_sus)+dose^(hill_amp_sus))*(1-isKO) + isKO*fold_as*(Max_amp_sus*dose^(hill_amp_sus))/(KD_amp_sus^(hill_amp_sus)+dose^(hill_amp_sus))" +timescale_sus "(Max_timescale_sus*KD_timescale_sus^(hill_timescale_sus))/(KD_timescale_sus^(hill_timescale_sus)+dose^(hill_timescale_sus))*(1-isKO) + isKO*fold_ts*(Max_timescale_sus*KD_timescale_sus^(hill_timescale_sus))/(KD_timescale_sus^(hill_timescale_sus)+dose^(hill_timescale_sus))" +toffset_TF_sus "(Max_toffset_TF_sus*KD_toffset_TF_sus^(hill_toffset_TF_sus))/(KD_toffset_TF_sus^(hill_toffset_TF_sus)+dose^(hill_toffset_TF_sus))*(1-isKO) + isKO*fold_to*(Max_toffset_TF_sus*KD_toffset_TF_sus^(hill_toffset_TF_sus))/(KD_toffset_TF_sus^(hill_toffset_TF_sus)+dose^(hill_toffset_TF_sus))" + + +PARAMETERS +//name value qFit dolog lb ub +//offset 0 1 0 -100 100 +//signum 1 0 0 -1 1 + + diff --git a/arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Setup.m b/arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Setup.m new file mode 100644 index 00000000..4f090e87 --- /dev/null +++ b/arFramework3/Examples/ToyModels/TransientFunction/Example4_FitDoseDependentFunctionToData/Setup.m @@ -0,0 +1,154 @@ +% addpath('../TransientFunction_library/') % functions for fitting the transient model (e.g. fitting both signums) + +nLHS = 500; + +%% LDH +arInit +arLoadModel('TransientDoseFunction') +% arLoadData('LDH',1) +arLoadData('LDH_WT',1) +arLoadData('LDH_KO',1) +arCompileAll(true) + +% Initialize_FitTransient % setting bounds, identification of parameter-types from name + +% arLoadPars(1) + +arSetPars('KD_toffset_TF_sus',log10(2), 1,1,-5,3); +arSetPars('Max_toffset_TF_sus',10, 1,0,-1,100); % unlog +arSetPars('hill_toffset_TF_sus',1, 1,0,1,5); +arSetPars('hill_timescale_sus',1, 1,0,1,5); +arSetPars('hill_amp_sus',1, 1,0,1,5); + +arFitLHS(nLHS) % multi-start fitting +% arFitTransient % single fit +arQplot('y') +arPlot +arPrint + +%% +arSave('LDH') + +arPLEInit +ple +pleExtend +pleSmooth +arSave('current') + +%% +% arLoad('20220815T133125_PLEs_Gut') +PlotHillCurves(2) +print -dpng HillCurves_LDH + +PlotHillCurves(3) +print -dpng HillCurves2_LDH + +PlotHillCurves_3D +PrintAllToPng(1:2,'HillCurves_LDH_3D') + +% PLE +PlotProfile('fold_as','max. amplitude ratio KO/WT') +title('LDH') +print -dpng Profile_fold_Amax_LDH + +PlotProfile('fold_ts','velocity ratio KO/WT') +title('LDH') +print -dpng Profile_fold_timescale_LDH + +PlotProfile('fold_to','delay ratio KO/WT') +title('LDH') +print -dpng Profile_fold_toffset_LDH + +PlotFits +suptitle('LDH') +print -dpng PlotFits_LDH + + +Plot3D +print -dpng 3D_LDH +% arSave('LDH_foldAufHill_LHS500_gut') + +%% IL1b +arInit +arLoadModel('TransientDoseFunction') +% arLoadData('LDH',1) +arLoadData('IL1b_WT',1) +arLoadData('IL1b_KO',1) +arCompileAll(true) + +% Initialize_FitTransient % setting bounds, identification of parameter-types from name + +% arLoadPars(1) + +arSetPars('KD_toffset_TF_sus',log10(2), 1,1,-5,3); +arSetPars('Max_toffset_TF_sus',10, 1,0,-1,100); % unlog +arSetPars('hill_toffset_TF_sus',1, 1,0,1,5); +arSetPars('hill_timescale_sus',1, 1,0,1,5); +arSetPars('hill_amp_sus',1, 1,0,1,5); + +arFitLHS(nLHS) % multi-start fitting +% arFitTransient % single fit +arQplot('y') +arPlot +arPrint + +%% +arSave('IL1b') + +arPLEInit +ple +pleExtend +arSave('current') + + +%% without t_offset +% indFix = [arPrint('fold_to'),arPrint('_toffset_')] +% ar.qFit(indFix) = 0; + +arFit +arSave('IL1b_toffsetFixed') + +arPLEInit +ple +pleExtend +pleSmooth +arSave('current') + +%% Fits IL1b +% arLoad('20220815T150340_IL1b_toffsetFixed') + +PlotFits +suptitle('IL1b') +print -dpng PlotFits_IL1b + +PlotProfile('fold_as','max. amplitude ratio KO/WT') +title('IL1b') +print -dpng Profile_fold_Amax_IL1b + +PlotProfile('fold_ts','velocity ratio KO/WT') +title('IL1b') +print -dpng Profile_fold_timescale_IL1b + +PlotProfile('fold_to','delay ratio KO/WT') +title('IL1b') +print -dpng Profile_fold_toffset_IL1b + + +PlotHillCurves(2) +print -dpng HillCurves_IL1b + +PlotHillCurves(3) +print -dpng HillCurves2_IL1b + +PlotHillCurves_3D +PrintAllToPng(1:2,'HillCurves_IL1b_3D') + +PlotHillCurves_3D('1./((Max_timescale_sus*KD_timescale_sus^(hill_timescale_sus))./(KD_timescale_sus^(hill_timescale_sus)+dosesFine.^(hill_timescale_sus))*(1-isKO) + isKO*fold_ts*(Max_timescale_sus*KD_timescale_sus^(hill_timescale_sus))./(KD_timescale_sus^(hill_timescale_sus)+dosesFine.^(hill_timescale_sus)));',... + 'Velocity'); +PlotHillCurves_3D('((Max_toffset_TF_sus*KD_toffset_TF_sus^(hill_toffset_TF_sus))./(KD_toffset_TF_sus^(hill_toffset_TF_sus)+dosesFine.^(hill_toffset_TF_sus))*(1-isKO) + isKO*fold_to*(Max_toffset_TF_sus*KD_toffset_TF_sus^(hill_toffset_TF_sus))./(KD_toffset_TF_sus^(hill_toffset_TF_sus)+dosesFine.^(hill_toffset_TF_sus)));',... + 'Max. Delay'); + +Plot3D +print -dpng 3D_IL1b + +% arSave('IL1b_foldAufHill_LHS500_gut') From 65434b3f9a0c5acb4a55ea384b3f0803a9022770 Mon Sep 17 00:00:00 2001 From: ckreutz Date: Tue, 12 Dec 2023 15:17:49 +0100 Subject: [PATCH 06/12] Modification arCreateDataStruct --- arFramework3/Advanced/arCreateDataStruct.m | 8 ++++++-- arFramework3/Advanced/arFindInputs.m | 2 +- arFramework3/Advanced/arRecompile.m | 7 +++++++ 3 files changed, 14 insertions(+), 3 deletions(-) diff --git a/arFramework3/Advanced/arCreateDataStruct.m b/arFramework3/Advanced/arCreateDataStruct.m index 0271419f..060cd82a 100644 --- a/arFramework3/Advanced/arCreateDataStruct.m +++ b/arFramework3/Advanced/arCreateDataStruct.m @@ -161,8 +161,12 @@ D.fp = transpose(D.p); % execute substitutions from ar.model.p/ar.model.fp -[int, iA, iB] = intersect(D.fp, ar.model.p); -D.fp(iA) = ar.model.fp(iB); +% try + [int, iA, iB] = intersect(D.fp, ar.model(m).p); +% catch +% ar.model.p +% end +D.fp(iA) = ar.model(m).fp(iB); %% now replace fp which was provided as function argument: [~,ia,ib] = intersect(D.p,pold); diff --git a/arFramework3/Advanced/arFindInputs.m b/arFramework3/Advanced/arFindInputs.m index 1c8dbacb..205ef78b 100644 --- a/arFramework3/Advanced/arFindInputs.m +++ b/arFramework3/Advanced/arFindInputs.m @@ -61,7 +61,7 @@ for c = 1 : length( step1 ) ar.model(m).condition(a).fu{b}(step1(c):end); - chk = strsplit(ar.model(m).condition(a).fu{b}(step1(c):end),','); + chk = strsplit(ar.model(m).condition(a).fu{b}(step1(c):end),',') stepLocations{end+1} = chk{3}; %#ok end for c = 1 : length( step2 ) diff --git a/arFramework3/Advanced/arRecompile.m b/arFramework3/Advanced/arRecompile.m index fcb092ec..958b6129 100644 --- a/arFramework3/Advanced/arRecompile.m +++ b/arFramework3/Advanced/arRecompile.m @@ -24,6 +24,13 @@ function arRecompile(sameParameterSettings, varargin) if isfield(ar,'setup') % only available in higher verions arIn = arDeepCopy(ar); + + if(~isfield(arIn.setup,'backup_model_folder_local')) + error('arRecompile: ar.setup.backup_model_folder_local does not exist. Save first via arSave.') + end + if(~isfield(arIn.setup,'backup_data_folder')) + error('arRecompile: ar.setup.backup_data_folder does not exist. Save first via arSave.') + end try for i=1:length(arIn.setup.commands) From 0ec58a977aef07e65f9f0eb2680e4c42c2ea5fea Mon Sep 17 00:00:00 2001 From: ckreutz Date: Tue, 12 Dec 2023 15:20:27 +0100 Subject: [PATCH 07/12] Updates transient function --- .../WriteTextSummary.m | 82 +++++++++++++++++++ ...proximateTimeCoursesByTransientFunction2.m | 17 ++-- .../arSetParsTransient.m | 42 ++++++++++ 3 files changed, 135 insertions(+), 6 deletions(-) create mode 100644 arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/WriteTextSummary.m create mode 100644 arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arSetParsTransient.m diff --git a/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/WriteTextSummary.m b/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/WriteTextSummary.m new file mode 100644 index 00000000..a38a5746 --- /dev/null +++ b/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/WriteTextSummary.m @@ -0,0 +1,82 @@ +% WriteTextSummary(fit, folder) +% +% fit a fit or cell array of fits as generated by +% arConditon2NewDataStruct*.m +% +% folder [WriteTextSummary] +% +% This function write matlab code for simulating the dynamics based on +% fitted transient functions in the specified folder. +% +% Example: +% arLoad +% fits = arApproximateTimeCoursesByTransientFunction2; +% WriteTextSummary(fits) +% t = linspace(0,100,101); +% plot(t,feval(fits{1}.label,t)); + +function WriteTextSummary(fit,folder) +if ~exist('folder','var') || isempty(folder) + folder = 'WriteTextSummary'; +end +if ~isdir(folder) + mkdir(folder); +end + +if iscell(fit) + for i=1:length(fit) + if isfield(fit{i},'pLabel') % otherwise fitting failed + WriteTextSummary(fit{i}); + else + fprintf('WriteTextSummary.m: It seems that fit{%i} fitting was failed.\n',i) + end + end + +else + if isfield(fit,'label') + name = sprintf('%s',fit.label); + else + name = sprintf('Model%i_%s_condition%i',fit.m, fit.x, fit.c); + end + name = strrep(name,' ','_'); + filename = [folder,filesep,name,'.m']; + fid = fopen(filename,'w'); + + fprintf(fid,'%s This function can be used to approximate model%i, condition%i, dynamic variable %s\n\n','%',fit.m, fit.c, fit.x); + fprintf(fid,'function [out, outSD] = %s(t)\n',name); + for i=1:length(fit.data.p) + if strcmp('init____dummy___',fit.data.p{i})~=1 % omit this D2D dummy parameter + if strcmp(fit.data.p{i},fit.data.fp{i})~=1 + fprintf(fid,'%s = %s;\n',fit.data.p{i},fit.data.fp{i}); + else + ind = strmatch(fit.data.p{i},fit.pLabel,'exact'); + fprintf(fid,'%s = %f;\n',fit.data.p{i},fit.p(ind)); + end + end + end + fprintf(fid,'\n'); + for i=1:length(fit.data.fu) + fprintf(fid,'%s = %s;\n',fit.u{i},Replacements(fit.data.fu{i})); + end + fprintf(fid,'\n'); + for i=1:length(fit.data.fy) + fprintf(fid,'out = %s;\n',Replacements(fit.data.fy{i})); + end + fprintf(fid,'\n'); + fprintf(fid,'outSD = NaN(size(out));\n'); + fprintf(fid,'outSD(:) = %f;\n',fit.approxErr); + fprintf(fid,'\n'); + + + + + + fclose(fid); +end + + +function out = Replacements(in) +out = strrep(in,'*','.*'); +out = strrep(out,'/','./'); +out = strrep(out,'^','.^'); + diff --git a/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arApproximateTimeCoursesByTransientFunction2.m b/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arApproximateTimeCoursesByTransientFunction2.m index 12193e5e..1990df3c 100644 --- a/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arApproximateTimeCoursesByTransientFunction2.m +++ b/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arApproximateTimeCoursesByTransientFunction2.m @@ -100,7 +100,7 @@ % load tmp fits % if specific fits should only be fitted %% for d=1:length(fits) -% try + try arInit; arLoadModel('TransientFunction_ForConditionFit2'); % fits{d}.data.fp{4} = '(100)' % for bachmann and conditions with too short tExp @@ -156,17 +156,22 @@ plotFitTransient(fits(d),plotfolder); end -% catch ERR -% pcatch{end+1} = struct('p',ar.p,'lb',ar.lb,'ub',ar.ub); -% save error ok pcatch -% end + catch ERR + pcatch{end+1} = struct('p',ar.p,'lb',ar.lb,'ub',ar.ub); + save error ok pcatch + end end +save arApproximateTimeCoursesByTransientFunction2 arPlot %% calculate the approximation error from the fitted SD for i=1:length(fits) - fits{i}.approxErr = 10.^fits{i}.p(4)./range(fits{i}.data.yExp); + if(isfield(fits{i},'')) + fits{i}.approxErr = 10.^fits{i}.p(4)./range(fits{i}.data.yExp); + else + fits{i}.approxErr = NaN*range(fits{i}.data.yExp); + end end diff --git a/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arSetParsTransient.m b/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arSetParsTransient.m new file mode 100644 index 00000000..f3407bc0 --- /dev/null +++ b/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arSetParsTransient.m @@ -0,0 +1,42 @@ +% This function calls arSetPars but also translates info (e.g. bounds) into +% ar.fit_transient +function arSetParsTransient(pLabel, varargin) +global ar +if length(varargin)>0 + p = varargin{1}; +else + p = []; +end +if length(varargin)>1 + qFit = varargin{2}; +else + qFit = []; +end +if length(varargin)>2 + qLog10 = varargin{3}; +else + qLog10 = []; +end +if length(varargin)>3 + lb = varargin{4}; +else + lb = []; +end +if length(varargin)>4 + ub = varargin{5}; +else + ub = []; +end + +% translate bounds to fit_transient.bounds +if(isfield(ar,'fit_transient')) + ind = strmatch(pLabel,ar.fit_transient.bounds.pLabel,'exact'); + if ~isemtpy(lb) + ar.fit_transient.bounds.lb(ind) = lb; + end + if ~isemtpy(ub) + ar.fit_transient.bounds.ub(ind) = ub; + end +end +% arSetPars(pLabel, [p], [qFit], [qLog10], [lb], [ub], [type], [meanp], [stdp]) +arSetPars(pLabel_or_ar, varargin{:}) \ No newline at end of file From e9dbb6851f4c4a0db8af4b3983f3497d7809390b Mon Sep 17 00:00:00 2001 From: Timo Rachel Date: Mon, 8 Jan 2024 17:23:00 +0100 Subject: [PATCH 08/12] typos and minot changes --- .../arApproximateTimeCoursesByTransientFunction2.m | 11 +++++++---- .../arConditon2NewDataStruct2.m | 2 +- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arApproximateTimeCoursesByTransientFunction2.m b/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arApproximateTimeCoursesByTransientFunction2.m index 1990df3c..9739cc66 100644 --- a/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arApproximateTimeCoursesByTransientFunction2.m +++ b/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arApproximateTimeCoursesByTransientFunction2.m @@ -61,9 +61,12 @@ end %% Step 1: Create data structs from conditions: -def_file = [fileparts(which('arInit')),filesep,'Examples',filesep,'ToyModels',filesep,'TransientFunction',filesep,'TransientFunction_library',filesep,'TransientFunction_ForConditionFit2.def']; +def_name = 'TransientFunction_ForConditionFit2'; +def_file = [fileparts(which('arInit')),filesep,'Examples',filesep,'ToyModels',filesep,'TransientFunction',filesep,'TransientFunction_library',filesep,def_name,'.def']; copyfile(def_file,'Models'); -arLoadModel('TransientFunction_ForConditionFit2'); +if ~strcmp(def_name, ar.model(end).name) + arLoadModel('TransientFunction_ForConditionFit2'); +end fits = cell(0); maxTF = []; for m=1:length(ar.model) @@ -128,8 +131,8 @@ SDest = 10.^SDest; end if sum(~isnan(ar.model.data.yExp))>1 && range(ar.model.data.yExp)/101e-10 - fprintf('Fit %i out of %i has large SD: LHS(10) started ...\n',d,length(fits)); - arFitLHS(20) % 10 fits should be enough for 7 parameters! + fprintf('Fit %i out of %i has still large SD: LHS(20) started ...\n',d,length(fits)); + arFitLHS(20) end end diff --git a/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arConditon2NewDataStruct2.m b/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arConditon2NewDataStruct2.m index 6378e235..77b20622 100644 --- a/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arConditon2NewDataStruct2.m +++ b/arFramework3/Examples/ToyModels/TransientFunction/TransientFunction_library/arConditon2NewDataStruct2.m @@ -52,7 +52,7 @@ mtarget = strmatch('TransientFunction_ForConditionFit2',{ar.model.name},'exact'); if length(mtarget)~=1 disp('In this implementation, the transient function has to be added to have ar.model.fy, ar.model.fystd etc '); - error('Model ''TransientFunction_ForConditionFit2'' not available. Please create a model structure and add ''TransientFunction'' it via arLoadModel.'); + error('Model ''TransientFunction_ForConditionFit2'' not available. Please create a model structure and add ''TransientFunction'' to it via arLoadModel.'); end end From d558983661c9c635bef67928680042e30504bae9 Mon Sep 17 00:00:00 2001 From: Franz-Georg Wieland Date: Wed, 17 Jan 2024 10:21:29 +0100 Subject: [PATCH 09/12] Option to hide windows in ple calculation Use flag hideWindows to hide calculation of ple and waitbar. Use flag ar.config.noWaitBarWindow to move waitbar to command line output in general. Use flag ar.config.noWaitBar to hide waitbar completely. --- arFramework3/MatlabTools/arWaitbar.m | 12 ++++++++---- arFramework3/PLE/ple.m | 14 ++++++++++++-- 2 files changed, 20 insertions(+), 6 deletions(-) diff --git a/arFramework3/MatlabTools/arWaitbar.m b/arFramework3/MatlabTools/arWaitbar.m index 3c9c468c..24b8b383 100644 --- a/arFramework3/MatlabTools/arWaitbar.m +++ b/arFramework3/MatlabTools/arWaitbar.m @@ -20,12 +20,16 @@ function arWaitbar(j, n, text) global arWaitbarGlobal; global arOutputLevel; +global ar; % disable waitbar window and use command line output instead -if ~isfield(arWaitbarGlobal,'showWindow') - arWaitbarGlobal.showWindow = 1; +if ( isfield( ar, 'config' ) && isfield( ar.config, 'noWaitBarWindow' ) && ar.config.noWaitBarWindow == 1 ) + arWaitbarGlobal.showWindow = 0; +else + if ~isfield(arWaitbarGlobal,'showWindow') + arWaitbarGlobal.showWindow = 1; + end end -global ar; if ( isfield( ar, 'config' ) && isfield( ar.config, 'noWaitBar' ) && ar.config.noWaitBar == 1 ) return; end @@ -116,7 +120,7 @@ function arWaitbar(j, n, text) sprintf('%s\n%i/%i %2i%% %s -> %s%s', ... strrep(text, '_','\_'), j, n, round(j/n*100), secToHMS(timeleft), secToHMS(timeelapsed), funtext)); else - fprintf('%s %i/%i %2i%% %s -> %s%s\n', ... + fprintf('%s %i/%i %2i%% Estimated: %s -> Elapsed: %s%s\n', ... strrep(text, '_','\_'), j, n, round(j/n*100), secToHMS(timeleft), secToHMS(timeelapsed), funtext); drawnow; end diff --git a/arFramework3/PLE/ple.m b/arFramework3/PLE/ple.m index 4907ff3f..42ed9d73 100755 --- a/arFramework3/PLE/ple.m +++ b/arFramework3/PLE/ple.m @@ -1,4 +1,4 @@ -% ple([i], [samplesize], [relchi2stepincrease], [maxstepsize], [minstepsize], [breakonbounds], [doLeftRightBranch]) +% ple([i], [samplesize], [relchi2stepincrease], [maxstepsize], [minstepsize], [breakonbounds], [doLeftRightBranch],[hideWindows]) % % Profile Likelihood Exploit % @@ -18,6 +18,8 @@ % break regardless for the profile parameter if it % reaches the boundary. % doLeftRightBranch calculate left/right branch [true true] +% hideWindows hide windows of calculation and [false] +% wait bar % % The profile likelihood calculation by the functions ple* was intended % as running independent of D2D, i.e. it was intended to be also used by @@ -29,7 +31,7 @@ % See also: arPLEInit, pleExtend, pleSmooth function ple(jk, samplesize, relchi2stepincrease, ... - maxstepsize, minstepsize, breakonlb, breakonub, doLeftRightBranch) + maxstepsize, minstepsize, breakonlb, breakonub, doLeftRightBranch,hideWindows) global ar @@ -89,6 +91,14 @@ function ple(jk, samplesize, relchi2stepincrease, ... doLeftRightBranch = [true true]; end +if(~exist('hideWindows','var') || isempty(hideWindows)) + hideWindows = false; +end +if hideWindows == true + ar.ple.showCalculation = false; + ar.config.noWaitBarWindow = true; +end + % Specifications on what to do if there are more than 1 parameter. The % function ple will evaluate itself for every parameter separately. if(nargin<1) From 8ee14b796df684fb4b4277d5276c420bd38fbeee Mon Sep 17 00:00:00 2001 From: Timo Rachel Date: Wed, 31 Jan 2024 16:08:01 +0100 Subject: [PATCH 10/12] updates on arLRT and arWaterfallPlot --- arFramework3/MatlabTools/arLRT.m | 16 +++++++++++++++- arFramework3/arWaterfallPlot.m | 8 ++++---- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/arFramework3/MatlabTools/arLRT.m b/arFramework3/MatlabTools/arLRT.m index 03d02f9c..e17e77cd 100644 --- a/arFramework3/MatlabTools/arLRT.m +++ b/arFramework3/MatlabTools/arLRT.m @@ -1,12 +1,26 @@ % [pval,q] = arLRT(DeltaChi2, DeltaDF, alpha) % % Performs Likelihood Ratio Test -% % DeltaDF degrees of freedom +% alpha significance level % % q = TRUE: accept H0, the smaller model can not be dismissed % q = FALSE: reject H0, dismiss smaller model in favour of larger model +% +% Example: (ar2 ist the larger model, m1 must be nested in m2, i.e. a special case) +% ar = ar1; +% arCalcMerit +% [~,m1] = arGetMerit; +% ar = ar2; +% arCalcMerit +% [~,m2] = arGetMerit; +% DeltaChi2 = m1.loglik_all - m2.loglik_all % must be positive +% DeltaDF = m2.npara - m1.npara; % must be positive + function [pval,q] = arLRT(DeltaChi2, DeltaDF, alpha) +if(~exist('alpha','var') || isempty(alpha)) + alpha = 0.05; +end pval = 1-chi2cdf(DeltaChi2, DeltaDF); q = pval > alpha; \ No newline at end of file diff --git a/arFramework3/arWaterfallPlot.m b/arFramework3/arWaterfallPlot.m index 18fd98ee..7b357ae7 100644 --- a/arFramework3/arWaterfallPlot.m +++ b/arFramework3/arWaterfallPlot.m @@ -13,17 +13,17 @@ % % see also arPlotChi2s -function arWaterfallPlot +function [fig, axis] = arWaterfallPlot global ar +fig = figure(); if(min(ar.chi2s)>0) - semilogy(1:length(ar.chi2s),sort(ar.chi2s),'.','MarkerSize',12) + axis = semilogy(1:length(ar.chi2s),sort(ar.chi2s),'.','MarkerSize',12); else - plot(1:length(ar.chi2s),sort(ar.chi2s),'.','MarkerSize',12) + axis = plot(1:length(ar.chi2s),sort(ar.chi2s),'.','MarkerSize',12); end xlabel('Fit index [sorted]') ylabel('Merit after optimization') xlim([0,length(ar.chi2s)+1]) title(['Waterfall plot: ',num2str(length(ar.chi2s)),' fits, ',num2str(sum(isinf(ar.chi2s))), ' failed']) grid on - From ceb53e30548e78b51afd2bc92ee0cf0d78029420 Mon Sep 17 00:00:00 2001 From: Andreas Raue Date: Tue, 16 Apr 2024 17:07:56 +0200 Subject: [PATCH 11/12] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 3615ae90..34c33b53 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # Data2Dynamics Software -**Contact:** Andreas Raue - +**Contact:** Andreas Raue - (for support issues, please use the issues and forum, thanks!) **Cite:** From 1a82e5bad57e30a8390d6caf3cae9f4cfb7d05d7 Mon Sep 17 00:00:00 2001 From: Niklas Neubrand <38315848+niklasneubrand@users.noreply.github.com> Date: Tue, 30 Apr 2024 13:04:47 +0200 Subject: [PATCH 12/12] Minor fixes regarding console output --- arFramework3/Advanced/arFits.m | 463 +++++++++++++------------- arFramework3/Subfunctions/arFprintf.m | 18 +- arFramework3/arCompileAll.m | 6 +- 3 files changed, 256 insertions(+), 231 deletions(-) diff --git a/arFramework3/Advanced/arFits.m b/arFramework3/Advanced/arFits.m index 8f12ecfc..139f894d 100755 --- a/arFramework3/Advanced/arFits.m +++ b/arFramework3/Advanced/arFits.m @@ -1,4 +1,4 @@ -% arFits(ps, [log_fit_history], [backup_save], [prefunc], [postfunc]) +% arFits(ps, [log_fit_history], [backup_save], [prefunc], [postfunc], [output_level]) % % Performing a multistart fit sequence with given initial guesses using arFit % @@ -12,6 +12,9 @@ % [false] % prefunc function to be called before each fitting % postfunc function to be called after each fitting +% output_level integer specifying the level of output ranging from 0-3 +% if empty, the default from arFprintf is used +% [empty] % % Performs a fitting sequence of the provided matrix of initial values for % the parameter values using arFit. ps has to have the right dimensions @@ -43,252 +46,266 @@ % % See also arFit, arFitLHS, arChi2LHS -function arFits(ps, log_fit_history, backup_save, prefunc, postfunc) +function arFits(ps, log_fit_history, backup_save, prefunc, postfunc, output_level) -global ar - -if(~exist('log_fit_history','var') || isempty(log_fit_history)) - log_fit_history = false; -end -if(~exist('backup_save','var') || isempty(backup_save)) - backup_save = false; -end -if(~exist('prefunc','var') || isempty(prefunc)) - prefunc = []; -end -if(~exist('postfunc','var') || isempty(postfunc)) - postfunc = []; -end - -if(~isfield(ar.config,'useFitErrorMatrix')) - ar.config.useFitErrorMatrix = false; -end - -dop = find(sum(~isnan(ps),2)==size(ps,2)); - -n = length(dop); -if ~isfield(ar,'ps') || size(ps,1) 1 - if isfield(ar.config,'showLiveWaterfall') && ar.config.showLiveWaterfall ==1 - arPlotChi2s - fig = gcf; - fig.Children(3).XLim = [0 length(ar.ps)]; - end + if(~exist('postfunc','var') || isempty(postfunc)) + postfunc = []; end - if j == n - close(figB) + if(~exist('output_level','var') || isempty(output_level)) + output_level = []; end - if handles.stop_now==1 - arFprintf(1,'\nThe fitting was stopped after %i/%i iterations manually by the user.\n',j-1,n); - break; + if(~isfield(ar.config,'useFitErrorMatrix')) + ar.config.useFitErrorMatrix = false; end - ar.p = ps(dop(j),:); - if(isfield(ar.config,'useDouble') && ar.config.useDouble==1) - ar.p(ar.iref) = ar.p(ar.iprimary); + dop = find(sum(~isnan(ps),2)==size(ps,2)); + + n = length(dop); + if ~isfield(ar,'ps') || size(ps,1) + ar.fit_hist = []; + end + + + arWaitbar(0); + + global handles + handles.stop_now = 0; + figB = buttonStop(); + + for j=1:n + arWaitbar(j, n); + if j > 1 + if isfield(ar.config,'showLiveWaterfall') && ar.config.showLiveWaterfall ==1 + arPlotChi2s + fig = gcf; + fig.Children(3).XLim = [0 length(ar.ps)]; + end + end + if j == n + % try statement allows to close the button manually without crash + try + close(figB) + end end - ar.fit_hist(dop(j)).hist = ar.fit; - ar.fit_hist(dop(j)).optimizer = ar.config.optimizer; - if(ar.config.optimizer==5) - ar.fit_hist(dop(j)).optimizerStep = ar.config.optimizerStep; - else - ar.fit_hist(dop(j)).optimizerStep = nan; + if handles.stop_now==1 + arFprintf(1,'\nThe fitting was stopped after %i/%i iterations manually by the user.\n',j-1,n); + break; + end + + ar.p = ps(dop(j),:); + if(isfield(ar.config,'useDouble') && ar.config.useDouble==1) + ar.p(ar.iref) = ar.p(ar.iprimary); end - ar.fit_hist(dop(j)).config = ar.config.optim; - ar.fit_hist(dop(j)).name = [name '_' sprintf('run%i', dop(j))]; - [~,imin] = min(ar.fit.chi2_hist + ar.fit.constr_hist); - ar.fit_hist(dop(j)).p = ar.fit.p_hist(imin,:); + tic; + try + arCalcMerit(true,ar.p(ar.qFit==1)); + ar.chi2s_start(dop(j)) = arGetMerit('chi2fit'); + ar.chi2sconstr_start(dop(j)) = arGetMerit('chi2constr'); + if ~isempty(prefunc) + try + feval( prefunc ); + catch + arFprintf(1, 'Error: Failure calling pre-fitting function'); + end + end + arFit(true); + if ~isempty(postfunc) + try + feval( postfunc ); + catch + arFprintf(1, 'Error: Failure calling post-fitting function'); + end + end + ar.ps(dop(j),:) = ar.p; + ar.chi2s(dop(j)) = arGetMerit('chi2fit'); + ar.chi2sconstr(dop(j)) = arGetMerit('chi2constr'); + ar.exitflag(dop(j)) = ar.fit.exitflag; + ar.fun_evals(dop(j)) = ar.fit.fevals; + ar.iter(dop(j)) = ar.fit.iter; + ar.optim_crit(dop(j)) = ar.firstorderopt; + catch exception + ar.chi2s(dop(j)) = inf; + ar.ps_errors(dop(j),:) = ar.p; + fprintf('fit #%i: %s\n', dop(j), exception.message); + end + + ar.timing(dop(j)) = toc; + if(isfield(ar, 'fit')) + if(isfield(ar.fit,'optimLog')) % coincides with ar.config.logFitting + ar.optimLogs{dop(j)} = ar.fit.optimLog; + end + end + + if(log_fit_history) + name = ar.config.optimizers{ar.config.optimizer}; + if(ar.config.optimizer==5) + tmpnames = arNLS; + name = [name '_' tmpnames{ar.config.optimizerStep+1}]; %#ok + end + + ar.fit_hist(dop(j)).hist = ar.fit; + ar.fit_hist(dop(j)).optimizer = ar.config.optimizer; + if(ar.config.optimizer==5) + ar.fit_hist(dop(j)).optimizerStep = ar.config.optimizerStep; + else + ar.fit_hist(dop(j)).optimizerStep = nan; + end + ar.fit_hist(dop(j)).config = ar.config.optim; + ar.fit_hist(dop(j)).name = [name '_' sprintf('run%i', dop(j))]; + + [~,imin] = min(ar.fit.chi2_hist + ar.fit.constr_hist); + ar.fit_hist(dop(j)).p = ar.fit.p_hist(imin,:); + end + if(backup_save) + save('arFits_backup.mat','ar'); + end end - if(backup_save) - save('arFits_backup.mat','ar'); - end -end - -arFprintf([],'total fitting time: %s\n', secToHMS(sum(ar.timing(~isnan(ar.timing))))); -arFprintf([],'mean fitting time: %s\n', secToHMS(10^mean(log10(ar.timing(~isnan(ar.timing)))))); -arWaitbar(-1); - -if(chi2Reset>min(ar.chi2s + ar.chi2sconstr)) - [chi2min,imin] = min(ar.chi2s + ar.chi2sconstr); - ar.p = ar.ps(imin,:); - if ar.config.fiterrors == -1 || (ar.config.fiterrors==0 && sum(ar.qFit(ar.qError==1)<2)==0) % if no error parameters fitted - arFprintf([],'selected best fit #%i with %f (old = %f)\n', imin, 2*(ar.ndata+ar.nconstr)*log(sqrt(2*pi)) + chi2min, 2*(ar.ndata+ar.nconstr)*log(sqrt(2*pi)) + chi2Reset); + + arFprintf(output_level,'total fitting time: %s\n', secToHMS(sum(ar.timing(~isnan(ar.timing))))); + arFprintf(output_level,'mean fitting time: %s\n', secToHMS(10^mean(log10(ar.timing(~isnan(ar.timing)))))); + arWaitbar(-1); + + if(chi2Reset>min(ar.chi2s + ar.chi2sconstr)) + [chi2min,imin] = min(ar.chi2s + ar.chi2sconstr); + ar.p = ar.ps(imin,:); + if ar.config.fiterrors == -1 || (ar.config.fiterrors==0 && sum(ar.qFit(ar.qError==1)<2)==0) % if no error parameters fitted + arFprintf(output_level,'selected best fit #%i with %f (old = %f)\n', ... + imin, 2*(ar.ndata+ar.nconstr)*log(sqrt(2*pi)) + chi2min, ... + 2*(ar.ndata+ar.nconstr)*log(sqrt(2*pi)) + chi2Reset); + else + arFprintf(output_level,'selected best fit #%i with %f (old = %f)\n', ... + imin, 2*ar.ndata*log(sqrt(2*pi)) + chi2min, ... + 2*ar.ndata*log(sqrt(2*pi)) + chi2Reset); + end else - arFprintf([],'selected best fit #%i with %f (old = %f)\n', ... - imin, 2*ar.ndata*log(sqrt(2*pi)) + chi2min, 2*ar.ndata*log(sqrt(2*pi)) + chi2Reset); + arFprintf(output_level,'did not find better fit\n'); + ar.p = pReset; end -else - arFprintf([],'did not find better fit\n'); - ar.p = pReset; -end -try - arCalcMerit(true,[]); -catch ERR - if strcmp(ERR.identifier,'d2d:arCollectRes:NaN_in_res')~=1 && strcmp(ERR.identifier,'d2d:arCollectRes:NaN_in_sres')~=1 - rethrow(ERR) - else % do not throw an error if NaNs are in residuals for the initial guess - warning(ERR.message) + try + arCalcMerit(true,[]); + catch ERR + if strcmp(ERR.identifier,'d2d:arCollectRes:NaN_in_res')~=1 && strcmp(ERR.identifier,'d2d:arCollectRes:NaN_in_sres')~=1 + rethrow(ERR) + else % do not throw an error if NaNs are in residuals for the initial guess + warning(ERR.message) + end end -end - - -function figButton = buttonStop() -figButton = []; -persistent batchmode -if isempty(batchmode) - ss = get(0, 'ScreenSize'); - if max(ss)<2 - batchmode = true; - else - batchmode = false; + + + function figButton = buttonStop() + figButton = []; + persistent batchmode + if isempty(batchmode) + ss = get(0, 'ScreenSize'); + if max(ss)<2 + batchmode = true; + else + batchmode = false; + end end -end - -if batchmode - return % do not show Waitbar in batch mode -end - - -% Create a figure window -figButton = uifigure; -figButton.Position = [608 258 321 118]; - -% Create a UI axes -% ax = uiaxes('Parent',fig,... -% 'Units','pixels',... -% 'Position', [104, 123, 300, 201]); - -% Create a push button -btn = uibutton(figButton,'push',... - 'Position',[100 50 100 22],... - 'ButtonPushedFcn', @(btn,event) stopButtonPushed(figButton)); -btn.Text = 'Stop fits'; - - -% Create the function for the ButtonPushedFcn callback -function stopButtonPushed(fig2) - global handles - handles.stop_now = 1; - close(fig2) + + if batchmode + return % do not show Waitbar in batch mode + end + + + % Create a figure window + figButton = uifigure; + figButton.Position = [608 258 321 118]; + + % Create a UI axes + % ax = uiaxes('Parent',fig,... + % 'Units','pixels',... + % 'Position', [104, 123, 300, 201]); + + % Create a push button + btn = uibutton(figButton,'push',... + 'Position',[100 50 100 22],... + 'ButtonPushedFcn', @(btn,event) stopButtonPushed(figButton)); + btn.Text = 'Stop fits'; + + + % Create the function for the ButtonPushedFcn callback + function stopButtonPushed(fig2) + global handles + handles.stop_now = 1; + % try statement allows to manually close the figure without causing a crash + try %#ok + close(fig2) + end + + \ No newline at end of file diff --git a/arFramework3/Subfunctions/arFprintf.m b/arFramework3/Subfunctions/arFprintf.m index f490e567..c358b9ec 100644 --- a/arFramework3/Subfunctions/arFprintf.m +++ b/arFramework3/Subfunctions/arFprintf.m @@ -1,19 +1,19 @@ % arFprintf([output_level], [varargin] ) -% +% % arFprintf wraps around fprintf. It uses global variable arOutputLevel to -% throw an output only, if output_level <= global arOutputLevel -% -% output_level <= arOutputLevel, then output is done +% throw an output only, if output_level <= global arOutputLevel +% +% output_level <= arOutputLevel, then output is done % 0 - All printed output except errors suppressed % 1 - Minimal output (only errors and requested info (arPrint etc)) % 2 - Regular output level (updates) % 3 - Verbose output (everything) % varargin is passed to fprintf and contains format and output arguments -% +% % If arOutputLevel is not set, it is set to arOutputLevel = 2. % A global variable is used to speed up the code. This function might be % called very frequently. -% +% % Examples: % >> global arOutputLevel % >> arOutputLevel = 2; @@ -22,7 +22,7 @@ % >> arFprintf(2,'Parameter %s = %.2f\n','kon',pi) % Parameter kon = 3.14 % >> arFprintf(3,'Parameter %s = %.2f\n','kon',pi) % no output -% +% @@ -32,6 +32,10 @@ function arFprintf( output_level, varargin ) arOutputLevel = 2; end +if isempty( output_level ) + output_level = 2; +end + if ( output_level <= arOutputLevel ) try % try because fprintf does surprisingly sometimes cause an error if matlab is started in -batch mode (observed for R2018b under Windows) fprintf( varargin{:} ); diff --git a/arFramework3/arCompileAll.m b/arFramework3/arCompileAll.m index 70909302..f200cd18 100644 --- a/arFramework3/arCompileAll.m +++ b/arFramework3/arCompileAll.m @@ -120,7 +120,11 @@ function arCompileAll(forcedCompile, debug_mode, source_dir, simplifyEquations) end % enable timedebug mode, use with debug_mode = true! -timedebug = false; +if debug_mode + timedebug = true; +else + timedebug = false; +end usePool = exist('gcp','file')>0 && ~isempty(gcp('nocreate'));