From 6c2986f1727b31863520426f1938767997555c65 Mon Sep 17 00:00:00 2001 From: "fabian.froehlich" Date: Wed, 14 Mar 2018 23:33:31 +0100 Subject: [PATCH 1/4] quick'n'dirty fix for #293 --- @amimodel/generateM.m | 6 ++++++ .../example_dirac_secondorder.m | 13 ++++++------- src/returndata_matlab.cpp | 2 -- 3 files changed, 12 insertions(+), 9 deletions(-) diff --git a/@amimodel/generateM.m b/@amimodel/generateM.m index 89e20ca5c5..0c18381da6 100644 --- a/@amimodel/generateM.m +++ b/@amimodel/generateM.m @@ -357,7 +357,9 @@ function generateM(this, amimodelo2) switch(o2flag) case 1 fprintf(fid,[' sol.s2x = reshape(sol.sx(:,' num2str(nxtrue+1) ':end,:),length(tout),' num2str(nxtrue) ',' num2str(np) ',length(options_ami.sens_ind));\n']); + fprintf(fid,[' sol.s2x = sol.s2x(:,:,options_ami.sens_ind,:);\n']); fprintf(fid,[' sol.s2y = reshape(sol.sy(:,' num2str(nytrue+1) ':end,:),length(tout),' num2str(nytrue) ',' num2str(np) ',length(options_ami.sens_ind));\n']); + fprintf(fid,[' sol.s2y = sol.s2y(:,:,options_ami.sens_ind,:);\n']); fprintf(fid,[' sol.s2sigmay = reshape(sol.ssigmay(:,' num2str(nytrue+1) ':end,:),length(tout),' num2str(nytrue) ',' num2str(np) ',length(options_ami.sens_ind));\n']); case 2 fprintf(fid,[' sol.s2x = sol.sx(:,' num2str(nxtrue+1) ':end,:);\n']); @@ -386,6 +388,10 @@ function generateM(this, amimodelo2) fprintf(fid,[' sol.s2rz(:,iz,:) = reshape(sol.srz(:,2*(iz-1)+2,:),options_ami.nmaxevent,1,length(theta(options_ami.sens_ind)));\n']); end fprintf(fid,' end\n'); + if(o2flag && nz > 0) + fprintf(fid,[' sol.s2z = sol.s2z(:,:,options_ami.sens_ind,:);\n']); + fprintf(fid,[' sol.s2rz = sol.s2rz(:,:,options_ami.sens_ind,:);\n']); + end fprintf(fid,[' sol.sx = sol.sx(:,1:' num2str(nxtrue) ',:);\n']); fprintf(fid,[' sol.sy = sol.sy(:,1:' num2str(nytrue) ',:);\n']); fprintf(fid,[' sol.ssigmay = sol.ssigmay(:,1:' num2str(nytrue) ',:);\n']); diff --git a/examples/example_dirac_secondorder/example_dirac_secondorder.m b/examples/example_dirac_secondorder/example_dirac_secondorder.m index cba5ba56cc..ec1a375f32 100644 --- a/examples/example_dirac_secondorder/example_dirac_secondorder.m +++ b/examples/example_dirac_secondorder/example_dirac_secondorder.m @@ -35,7 +35,6 @@ function example_dirac_secondorder() optionsfd = options; optionsfd.sensi = 1; -optionsfd.sens_ind = 1:4; eps = 1e-4; xi = log10(p); @@ -53,9 +52,9 @@ function example_dirac_secondorder() if(usejava('jvm')) figure c_x = get(gca,'ColorOrder'); - for ip = 1:4 + for ip = 1:length(options.sens_ind) for jp = 1:length(options.sens_ind) - subplot(length(options.sens_ind),4,(jp-1)*4+ip) + subplot(length(options.sens_ind),length(options.sens_ind),(jp-1)*length(options.sens_ind)+ip) hold on for ix = 1:size(sol.x,2) plot(t,sol.s2x(:,ix,ip,jp),'.-','Color',c_x(ix,:)) @@ -64,7 +63,7 @@ function example_dirac_secondorder() ylim([-10,10]) legend('x1','x1_{fd}','x2','x2_{fd}','Location','NorthEastOutside') legend boxoff - title(['state sensitivity for p' num2str(ip) '-p' num2str(options.sens_ind(jp))]) + title(['state sensitivity for p' num2str(options.sens_ind(ip)) '-p' num2str(options.sens_ind(jp))]) xlabel('time t') ylabel('x') box on @@ -72,13 +71,13 @@ function example_dirac_secondorder() end set(gcf,'Position',[100 300 1200 500]) figure - for ip = 1:4 + for ip = 1:length(options.sens_ind) for jp = 1:length(options.sens_ind) - subplot(length(options.sens_ind),4,(jp-1)*4+ip) + subplot(length(options.sens_ind),length(options.sens_ind),(jp-1)*length(options.sens_ind)+ip) plot(t,abs(sol.s2x(:,:,ip,jp)-s2x_fd(:,:,ip,jp)),'r--') legend('error x1','error x2','Location','NorthEastOutside') legend boxoff - title(['state sensitivity for p' num2str(ip) '-p' num2str(options.sens_ind(jp))]) + title(['state sensitivity for p' num2str(options.sens_ind(ip)) '-p' num2str(options.sens_ind(jp))]) xlabel('time t') ylabel('error') ylim([1e-12,1e0]) diff --git a/src/returndata_matlab.cpp b/src/returndata_matlab.cpp index 090f8fc075..002a468989 100644 --- a/src/returndata_matlab.cpp +++ b/src/returndata_matlab.cpp @@ -179,7 +179,6 @@ void writeMatlabField0(mxArray *matlabStruct, const char *fieldName, double *array = initAndAttachArray(matlabStruct, fieldName, dim); - /* transform rowmajor (c++) to colmajor (matlab) */ array[0] = static_cast(fielddata); } @@ -200,7 +199,6 @@ void writeMatlabField1(mxArray *matlabStruct, const char *fieldName, double *array = initAndAttachArray(matlabStruct, fieldName, dim); - /* transform rowmajor (c++) to colmajor (matlab) */ for(int i = 0; i < dim0; i++) array[i] = static_cast(fieldData[i]); } From e5807f6551d677182c446b654f1dae3fed5171d8 Mon Sep 17 00:00:00 2001 From: "fabian.froehlich" Date: Thu, 15 Mar 2018 10:16:28 +0100 Subject: [PATCH 2/4] revert previous changes and properly compute augcoefficient --- @amimodel/generateM.m | 6 ------ .../example_dirac_secondorder.m | 13 +++++++------ src/rdata.cpp | 19 +++++++++++++------ 3 files changed, 20 insertions(+), 18 deletions(-) diff --git a/@amimodel/generateM.m b/@amimodel/generateM.m index 0c18381da6..89e20ca5c5 100644 --- a/@amimodel/generateM.m +++ b/@amimodel/generateM.m @@ -357,9 +357,7 @@ function generateM(this, amimodelo2) switch(o2flag) case 1 fprintf(fid,[' sol.s2x = reshape(sol.sx(:,' num2str(nxtrue+1) ':end,:),length(tout),' num2str(nxtrue) ',' num2str(np) ',length(options_ami.sens_ind));\n']); - fprintf(fid,[' sol.s2x = sol.s2x(:,:,options_ami.sens_ind,:);\n']); fprintf(fid,[' sol.s2y = reshape(sol.sy(:,' num2str(nytrue+1) ':end,:),length(tout),' num2str(nytrue) ',' num2str(np) ',length(options_ami.sens_ind));\n']); - fprintf(fid,[' sol.s2y = sol.s2y(:,:,options_ami.sens_ind,:);\n']); fprintf(fid,[' sol.s2sigmay = reshape(sol.ssigmay(:,' num2str(nytrue+1) ':end,:),length(tout),' num2str(nytrue) ',' num2str(np) ',length(options_ami.sens_ind));\n']); case 2 fprintf(fid,[' sol.s2x = sol.sx(:,' num2str(nxtrue+1) ':end,:);\n']); @@ -388,10 +386,6 @@ function generateM(this, amimodelo2) fprintf(fid,[' sol.s2rz(:,iz,:) = reshape(sol.srz(:,2*(iz-1)+2,:),options_ami.nmaxevent,1,length(theta(options_ami.sens_ind)));\n']); end fprintf(fid,' end\n'); - if(o2flag && nz > 0) - fprintf(fid,[' sol.s2z = sol.s2z(:,:,options_ami.sens_ind,:);\n']); - fprintf(fid,[' sol.s2rz = sol.s2rz(:,:,options_ami.sens_ind,:);\n']); - end fprintf(fid,[' sol.sx = sol.sx(:,1:' num2str(nxtrue) ',:);\n']); fprintf(fid,[' sol.sy = sol.sy(:,1:' num2str(nytrue) ',:);\n']); fprintf(fid,[' sol.ssigmay = sol.ssigmay(:,1:' num2str(nytrue) ',:);\n']); diff --git a/examples/example_dirac_secondorder/example_dirac_secondorder.m b/examples/example_dirac_secondorder/example_dirac_secondorder.m index ec1a375f32..cba5ba56cc 100644 --- a/examples/example_dirac_secondorder/example_dirac_secondorder.m +++ b/examples/example_dirac_secondorder/example_dirac_secondorder.m @@ -35,6 +35,7 @@ function example_dirac_secondorder() optionsfd = options; optionsfd.sensi = 1; +optionsfd.sens_ind = 1:4; eps = 1e-4; xi = log10(p); @@ -52,9 +53,9 @@ function example_dirac_secondorder() if(usejava('jvm')) figure c_x = get(gca,'ColorOrder'); - for ip = 1:length(options.sens_ind) + for ip = 1:4 for jp = 1:length(options.sens_ind) - subplot(length(options.sens_ind),length(options.sens_ind),(jp-1)*length(options.sens_ind)+ip) + subplot(length(options.sens_ind),4,(jp-1)*4+ip) hold on for ix = 1:size(sol.x,2) plot(t,sol.s2x(:,ix,ip,jp),'.-','Color',c_x(ix,:)) @@ -63,7 +64,7 @@ function example_dirac_secondorder() ylim([-10,10]) legend('x1','x1_{fd}','x2','x2_{fd}','Location','NorthEastOutside') legend boxoff - title(['state sensitivity for p' num2str(options.sens_ind(ip)) '-p' num2str(options.sens_ind(jp))]) + title(['state sensitivity for p' num2str(ip) '-p' num2str(options.sens_ind(jp))]) xlabel('time t') ylabel('x') box on @@ -71,13 +72,13 @@ function example_dirac_secondorder() end set(gcf,'Position',[100 300 1200 500]) figure - for ip = 1:length(options.sens_ind) + for ip = 1:4 for jp = 1:length(options.sens_ind) - subplot(length(options.sens_ind),length(options.sens_ind),(jp-1)*length(options.sens_ind)+ip) + subplot(length(options.sens_ind),4,(jp-1)*4+ip) plot(t,abs(sol.s2x(:,:,ip,jp)-s2x_fd(:,:,ip,jp)),'r--') legend('error x1','error x2','Location','NorthEastOutside') legend boxoff - title(['state sensitivity for p' num2str(options.sens_ind(ip)) '-p' num2str(options.sens_ind(jp))]) + title(['state sensitivity for p' num2str(ip) '-p' num2str(options.sens_ind(jp))]) xlabel('time t') ylabel('error') ylim([1e-12,1e0]) diff --git a/src/rdata.cpp b/src/rdata.cpp index 91b41a8ee2..c9d8fceac1 100644 --- a/src/rdata.cpp +++ b/src/rdata.cpp @@ -162,21 +162,28 @@ void ReturnData::applyChainRuleFactorToSimulationResults(const Model *model) { std::vector unscaledParameters(np); model->unscaleParameters(unscaledParameters.data()); std::vector augcoefficient(np); + + if (sensi == AMICI_SENSI_ORDER_SECOND && o2mode == AMICI_O2MODE_FULL) { + for (int ip = 0; ip < np; ++ip) { + switch (pscale[ip]) { + case AMICI_SCALING_LOG10: + augcoefficient.at(ip) = unscaledParameters.at(ip) * log(10); + case AMICI_SCALING_LN: + augcoefficient.at(ip) = unscaledParameters.at(ip); + case AMICI_SCALING_NONE: + break; + } + } + } for (int ip = 0; ip < nplist; ++ip) { switch (pscale[model->plist(ip)]) { case AMICI_SCALING_LOG10: coefficient.at(ip) = log(10.0); pcoefficient.at(ip) = unscaledParameters.at(model->plist(ip)) * log(10); - if (sensi == AMICI_SENSI_ORDER_SECOND && o2mode == AMICI_O2MODE_FULL) - augcoefficient.at(ip) = unscaledParameters.at(ip) * log(10); - break; case AMICI_SCALING_LN: coefficient.at(ip) = 1.0; pcoefficient.at(ip) = unscaledParameters.at(model->plist(ip)); - if (sensi == AMICI_SENSI_ORDER_SECOND && o2mode == AMICI_O2MODE_FULL) - augcoefficient.at(ip) = unscaledParameters.at(ip); - break; case AMICI_SCALING_NONE: coefficient.at(ip) = 1.0; break; From 0accc3a137e0fe47f274ebbbe4dccb77ad13de58 Mon Sep 17 00:00:00 2001 From: "fabian.froehlich" Date: Thu, 15 Mar 2018 10:27:15 +0100 Subject: [PATCH 3/4] add breaks --- src/rdata.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/rdata.cpp b/src/rdata.cpp index c9d8fceac1..fd06f24e4a 100644 --- a/src/rdata.cpp +++ b/src/rdata.cpp @@ -168,9 +168,12 @@ void ReturnData::applyChainRuleFactorToSimulationResults(const Model *model) { switch (pscale[ip]) { case AMICI_SCALING_LOG10: augcoefficient.at(ip) = unscaledParameters.at(ip) * log(10); + break; case AMICI_SCALING_LN: augcoefficient.at(ip) = unscaledParameters.at(ip); + break; case AMICI_SCALING_NONE: + augcoefficient.at(ip) = 1; break; } } @@ -181,9 +184,11 @@ void ReturnData::applyChainRuleFactorToSimulationResults(const Model *model) { case AMICI_SCALING_LOG10: coefficient.at(ip) = log(10.0); pcoefficient.at(ip) = unscaledParameters.at(model->plist(ip)) * log(10); + break; case AMICI_SCALING_LN: coefficient.at(ip) = 1.0; pcoefficient.at(ip) = unscaledParameters.at(model->plist(ip)); + break; case AMICI_SCALING_NONE: coefficient.at(ip) = 1.0; break; From 1fc548aef403669cfb66144a6ab32b990334662d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Thu, 15 Mar 2018 10:40:01 +0100 Subject: [PATCH 4/4] use default initialization for chain rule coeffs coefficient, pcoefficient and augcoefficient --- src/rdata.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/rdata.cpp b/src/rdata.cpp index fd06f24e4a..18973e7b79 100644 --- a/src/rdata.cpp +++ b/src/rdata.cpp @@ -161,7 +161,7 @@ void ReturnData::applyChainRuleFactorToSimulationResults(const Model *model) { std::vector pcoefficient(nplist, 1.0); std::vector unscaledParameters(np); model->unscaleParameters(unscaledParameters.data()); - std::vector augcoefficient(np); + std::vector augcoefficient(np, 1.0); if (sensi == AMICI_SENSI_ORDER_SECOND && o2mode == AMICI_O2MODE_FULL) { for (int ip = 0; ip < np; ++ip) { @@ -173,7 +173,6 @@ void ReturnData::applyChainRuleFactorToSimulationResults(const Model *model) { augcoefficient.at(ip) = unscaledParameters.at(ip); break; case AMICI_SCALING_NONE: - augcoefficient.at(ip) = 1; break; } } @@ -186,11 +185,9 @@ void ReturnData::applyChainRuleFactorToSimulationResults(const Model *model) { pcoefficient.at(ip) = unscaledParameters.at(model->plist(ip)) * log(10); break; case AMICI_SCALING_LN: - coefficient.at(ip) = 1.0; pcoefficient.at(ip) = unscaledParameters.at(model->plist(ip)); break; case AMICI_SCALING_NONE: - coefficient.at(ip) = 1.0; break; } }