Skip to content

Commit

Permalink
clang-format
Browse files Browse the repository at this point in the history
  • Loading branch information
@dweindl
dweindl committed Oct 21, 2024
1 parent 9a1564e commit 8b4393d
Show file tree
Hide file tree
Showing 98 changed files with 6,717 additions and 5,769 deletions.
202 changes: 113 additions & 89 deletions examples/parpeamici/steadystate/exampleSteadystateScaledTest.h
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,15 @@
#include <gtest/gtest.h>

#include <parpecommon/parpeConfig.h>
#include <parpecommon/misc.h>
#include <parpeamici/multiConditionProblem.h>
#include <parpeamici/multiConditionDataProvider.h>
#include <parpeamici/multiConditionProblem.h>
#include <parpecommon/misc.h>
#include <parpecommon/parpeConfig.h>

#include "../../../tests/parpecommon/testingMisc.h"

#ifdef PARPE_ENABLE_IPOPT
#include <parpeoptimization/localOptimizationIpopt.h>
#include <parpeamici/hierarchicalOptimization.h>
#include <parpeoptimization/localOptimizationIpopt.h>

#include "wrapfunctions.h"

Expand All @@ -19,8 +19,7 @@

class steadystateProblemTests : public ::testing::Test {

protected:

protected:
/*
const std::vector<double> t { 1.0e8 };
const std::vector<double> k { 0.1, 0.4, 0.7, 1.0 };
Expand All @@ -29,23 +28,26 @@ class steadystateProblemTests : public ::testing::Test {
0.437977375496898,
0.033333333333333};
*/
const int scalingParameterIdx = 5;
const int offsetParameterIdx = 6;
const int scaledObservableIdx = 3;
const int offsettedObservableIdx = 4;
const std::vector<double> t { 1.0e8 };
int const scalingParameterIdx = 5;
int const offsetParameterIdx = 6;
int const scaledObservableIdx = 3;
int const offsettedObservableIdx = 4;
std::vector<double> const t{1.0e8};
// const std::vector<double> k { };
//const std::vector<double> p { 1.0, 0.5, 0.4, 2.0, 0.1, 0.2, 0.2, 0.2, 2.0, 0.2, 3.0, 0.2, 0.2 };
const std::vector<double> x0 { 0.1, 0.4, 0.7 };
const std::vector<double> xSteadystateExp {0.456644592142075,
0.437977375496898,
0.033333333333333};
const std::vector<double> yExp {0.456644592142075,
0.437977375496898,
0.033333333333333,
2.0 * 0.456644592142075,
3.0 + 0.437977375496898,
0.456644592142075};
// const std::vector<double> p { 1.0, 0.5, 0.4, 2.0, 0.1, 0.2, 0.2,
// 0.2, 2.0, 0.2, 3.0, 0.2, 0.2 };
std::vector<double> const x0{0.1, 0.4, 0.7};
std::vector<double> const xSteadystateExp{
0.456644592142075,
0.437977375496898,
0.033333333333333};
std::vector<double> const yExp{
0.456644592142075,
0.437977375496898,
0.033333333333333,
2.0 * 0.456644592142075,
3.0 + 0.437977375496898,
0.456644592142075};
};

TEST_F(steadystateProblemTests, testSteadystate) {
Expand All @@ -61,22 +63,26 @@ TEST_F(steadystateProblemTests, testSteadystate) {
auto rdata = amici::runAmiciSimulation(*solver, nullptr, *model);

// verify steadystate concentrations
parpe::checkEqualArray(xSteadystateExp.data(),
rdata->x.data(),
static_cast<int>(xSteadystateExp.size()),
1e-5, 1e-5);
parpe::checkEqualArray(
xSteadystateExp.data(),
rdata->x.data(),
static_cast<int>(xSteadystateExp.size()),
1e-5,
1e-5);

// verify likelihood for matching measurement / simulation
amici::ExpData edata {*model};
amici::ExpData edata{*model};
edata.setObservedData(yExp);
edata.setObservedDataStdDev(std::vector<double>(yExp.size(), 1.0));
rdata = amici::runAmiciSimulation(*solver, &edata, *model);

EXPECT_EQ(rdata->status, amici::AMICI_SUCCESS);
EXPECT_NEAR(1e-5, rdata->chi2, 1e-5);

EXPECT_NEAR(parpe::getLogLikelihoodOffset(edata.nt() * edata.nytrue()),
rdata->llh, 1e-5);
EXPECT_NEAR(
parpe::getLogLikelihoodOffset(edata.nt() * edata.nytrue()),
rdata->llh,
1e-5);
}

TEST_F(steadystateProblemTests, testSteadystateMultiCond) {
Expand All @@ -90,57 +96,62 @@ TEST_F(steadystateProblemTests, testSteadystateMultiCond) {

model->setTimepoints(t);
model->setInitialStates(x0);
//model->setParameters(p);
// model->setParameters(p);

parpe::MultiConditionDataProviderDefault dp(std::move(model),
modelNonOwning->getSolver());
parpe::MultiConditionDataProviderDefault dp(
std::move(model), modelNonOwning->getSolver());

dp.edata_.emplace_back(amici::ExpData(*modelNonOwning));
dp.edata_[0].fixedParameters = modelNonOwning->getFixedParameters();
dp.edata_[0].setObservedData(yExp);
dp.edata_[0].setObservedDataStdDev(std::vector<double>(yExp.size(), 1.0));

//parpe::AmiciSummedGradientFunction<int>(&dp, nullptr);
// parpe::AmiciSummedGradientFunction<int>(&dp, nullptr);
parpe::MultiConditionProblem problem(&dp);
double cost;
problem.cost_fun_->evaluate(p, cost, gsl::span<double>());
EXPECT_NEAR(-parpe::getLogLikelihoodOffset(
dp.edata_[0].getObservedData().size()), cost, 1e-5);
EXPECT_NEAR(
-parpe::getLogLikelihoodOffset(dp.edata_[0].getObservedData().size()),
cost,
1e-5);
}


TEST_F(steadystateProblemTests, testSteadystateHierarchical) {
// introduce scaling parameters
auto model = amici::generic_model::getModel();
//model->setFixedParameters(k);
// model->setFixedParameters(k);
model->setInitialStates(x0);
//model->setParameters(p);
// model->setParameters(p);
model->setTimepoints(t);
auto modelNonOwning = model.get();

const double scalingExp = 2.0; // scaling parameter
const double offsetExp = 2.0; // offset parameter
const std::vector<double> pReduced { 1.0, 0.5, 0.4, 2.0, 0.1, /*2.0, 3.0,*/ 0.2, 4.0 };
double const scalingExp = 2.0; // scaling parameter
double const offsetExp = 2.0; // offset parameter
std::vector<double> const pReduced{
1.0, 0.5, 0.4, 2.0, 0.1, /*2.0, 3.0,*/ 0.2, 4.0};
auto yScaledExp = yExp;
yScaledExp[scaledObservableIdx] = scalingExp * yExp[0];
yScaledExp[offsettedObservableIdx] = offsetExp + yExp[1];
parpe::MultiConditionDataProviderDefault dp(std::move(model), modelNonOwning->getSolver());
parpe::MultiConditionDataProviderDefault dp(
std::move(model), modelNonOwning->getSolver());
// x0?
dp.edata_.emplace_back(amici::ExpData(*modelNonOwning));
dp.edata_[0].fixedParameters = modelNonOwning->getFixedParameters();
dp.edata_[0].setObservedData(yScaledExp);
dp.edata_[0].setObservedDataStdDev(std::vector<double>(yExp.size(), 1.0));

//parpe::MultiConditionProblem problem(&dp);
// parpe::MultiConditionProblem problem(&dp);

auto scalings = std::make_unique<parpe::AnalyticalParameterProviderDefault>();
auto scalings =
std::make_unique<parpe::AnalyticalParameterProviderDefault>();
scalings->conditionsForParameter.push_back({0});
scalings->optimizationParameterIndices.push_back(scalingParameterIdx);
// x[scalingIdx][conditionIdx] -> std::vector of observableIndicies
scalings->mapping.resize(1);
scalings->mapping[0][0] = {scaledObservableIdx};

auto offsets = std::make_unique<parpe::AnalyticalParameterProviderDefault>();
auto offsets =
std::make_unique<parpe::AnalyticalParameterProviderDefault>();
offsets->conditionsForParameter.push_back({0});
offsets->optimizationParameterIndices.push_back(offsetParameterIdx);
// x[scalingIdx][conditionIdx] -> std::vector of observableIndicies
Expand All @@ -149,52 +160,55 @@ TEST_F(steadystateProblemTests, testSteadystateHierarchical) {

auto sigmas = std::make_unique<parpe::AnalyticalParameterProviderDefault>();

auto gradFun = std::make_unique<parpe::AmiciSummedGradientFunction>(&dp, nullptr, nullptr);
parpe::HierarchicalOptimizationWrapper hier(gradFun.get(),
std::move(scalings),
std::move(offsets),
std::move(sigmas),
dp.getNumberOfSimulationConditions(),
modelNonOwning->nytrue,
parpe::ErrorModel::normal);
auto gradFun = std::make_unique<parpe::AmiciSummedGradientFunction>(
&dp, nullptr, nullptr);
parpe::HierarchicalOptimizationWrapper hier(
gradFun.get(),
std::move(scalings),
std::move(offsets),
std::move(sigmas),
dp.getNumberOfSimulationConditions(),
modelNonOwning->nytrue,
parpe::ErrorModel::normal);
// TODO: need to adapt to changed model
// double cost;
// hier.evaluate(pReduced, cost, gsl::span<double>(), nullptr, nullptr);
// EXPECT_NEAR(-parpe::getLogLikelihoodOffset(
// dp.edata[0].getObservedData().size()), cost, 1e-5);

// const std::vector<double> pFull { 1.0, 0.5, 0.4, 2.0,
// 0.1, scalingExp, offsetExp, 0.2, 4.0 };
// hier.fun->evaluate(pFull, {0}, cost, gsl::span<double>(), nullptr, nullptr);
// EXPECT_NEAR(-parpe::getLogLikelihoodOffset(
// dp.edata[0].getObservedData().size()), cost, 1e-5);
// double cost;
// hier.evaluate(pReduced, cost, gsl::span<double>(), nullptr, nullptr);
// EXPECT_NEAR(-parpe::getLogLikelihoodOffset(
// dp.edata[0].getObservedData().size()), cost, 1e-5);

// const std::vector<double> pFull { 1.0, 0.5, 0.4, 2.0,
// 0.1, scalingExp, offsetExp, 0.2, 4.0
// };
// hier.fun->evaluate(pFull, {0}, cost, gsl::span<double>(), nullptr,
// nullptr); EXPECT_NEAR(-parpe::getLogLikelihoodOffset(
// dp.edata[0].getObservedData().size()), cost, 1e-5);
}



//TEST(steadystateProblemTests, testOptimizationHierarchical) {
// /* setup model & solver */
// // introduce scaling parameters
// auto model = getModel();
// //model->setFixedParameters(k);
// //model->setInitialStates(x0);
// //model->setParameters(p);
// model->setTimepoints(t);
// model->requireSensitivitiesForAllParameters();
// auto modelNonOwning = model.get();
// TEST(steadystateProblemTests, testOptimizationHierarchical) {
// /* setup model & solver */
// // introduce scaling parameters
// auto model = getModel();
// //model->setFixedParameters(k);
// //model->setInitialStates(x0);
// //model->setParameters(p);
// model->setTimepoints(t);
// model->requireSensitivitiesForAllParameters();
// auto modelNonOwning = model.get();

// auto solver = model->getSolver();
// solver->setSensitivityMethod(amici::SensitivityMethod::adjoint);

// /* generate scaled data */
// const double scalingExp = 2.0; // scaling parameter
// const double offsetExp = 2.0; // offset parameter
// const std::vector<double> pReduced { 1.0, 0.5, 0.4, /*2.0, 0.1,*/ 2.0, 3.0, 0.2, 4.0 };
// const std::vector<double> pReduced { 1.0, 0.5, 0.4, /*2.0,
// 0.1,*/ 2.0, 3.0, 0.2, 4.0 };

// auto yScaledExp = yExp;
// yScaledExp[scaledObservableIdx] = scalingExp * yExp[0];
// yScaledExp[offsettedObservableIdx] = offsetExp + yExp[1];
// parpe::MultiConditionDataProviderDefault dp(std::move(model), std::move(solver));
// parpe::MultiConditionDataProviderDefault dp(std::move(model),
// std::move(solver));
// // x0?
// dp.edata.push_back(amici::ExpData(*modelNonOwning));
// dp.edata[0].fixedParameters = modelNonOwning->getFixedParameters();
Expand All @@ -204,26 +218,30 @@ TEST_F(steadystateProblemTests, testSteadystateHierarchical) {

// /* setup hierarchical optimization */
// // one scaling parameter
// auto scalings = std::make_unique<parpe::AnalyticalParameterProviderDefault>();
// auto scalings =
// std::make_unique<parpe::AnalyticalParameterProviderDefault>();
// scalings->conditionsForParameter.push_back({0});
// scalings->optimizationParameterIndices.push_back(scalingParameterIdx);
// // x[scalingIdx][conditionIdx] -> std::vector of observableIndicies
// scalings->mapping.resize(1);
// scalings->mapping[0][0] = {scaledObservableIdx};

// // analytical offset parameter
// auto offsets = std::make_unique<parpe::AnalyticalParameterProviderDefault>();
// auto offsets =
// std::make_unique<parpe::AnalyticalParameterProviderDefault>();
// offsets->conditionsForParameter.push_back({0});
// offsets->optimizationParameterIndices.push_back(offsetParameterIdx);
// // x[scalingIdx][conditionIdx] -> std::vector of observableIndicies
// offsets->mapping.resize(1);
// offsets->mapping[0][0] = {offsettedObservableIdx};

// auto sigmas = std::make_unique<parpe::AnalyticalParameterProviderDefault>();
// auto sigmas =
// std::make_unique<parpe::AnalyticalParameterProviderDefault>();

// // create wrapper
// auto gradFun = std::make_unique<parpe::AmiciSummedGradientFunction>(&dp, nullptr, nullptr);
// auto hier = std::make_unique<parpe::HierarchicalOptimizationWrapper>(std::move(gradFun),
// auto gradFun = std::make_unique<parpe::AmiciSummedGradientFunction>(&dp,
// nullptr, nullptr); auto hier =
// std::make_unique<parpe::HierarchicalOptimizationWrapper>(std::move(gradFun),
// std::move(scalings),
// std::move(offsets),
// std::move(sigmas),
Expand All @@ -234,14 +252,19 @@ TEST_F(steadystateProblemTests, testSteadystateHierarchical) {
// // evaluate and ensure scaling factor is computed so that y_mes = y_sim
// double cost;
// hier->evaluate(pReduced, cost, gsl::span<double>(), nullptr, nullptr);
// DOUBLES_EQUAL(-parpe::getLogLikelihoodOffset(dp.edata[0].getObservedData().size()), cost, 1e-5);
// DOUBLES_EQUAL(-parpe::getLogLikelihoodOffset(dp.edata[0].getObservedData().size()),
// cost, 1e-5);

// const std::vector<double> pFull { 1.0, 0.5, 0.4, 2.0,
// 0.1, scalingExp, offsetExp, 1.0, 0.2, 4.0 };
// hier->fun->evaluate(pFull, {0}, cost, gsl::span<double>(), nullptr, nullptr);
// DOUBLES_EQUAL(-parpe::getLogLikelihoodOffset(dp.edata[0].getObservedData().size()), cost, 1e-5);

// parpe::OptimizationProblemImpl problem(std::move(hier), std::make_unique<parpe::Logger>());
// 0.1, scalingExp, offsetExp, 1.0,
// 0.2, 4.0 };
// hier->fun->evaluate(pFull, {0}, cost, gsl::span<double>(), nullptr,
// nullptr);
// DOUBLES_EQUAL(-parpe::getLogLikelihoodOffset(dp.edata[0].getObservedData().size()),
// cost, 1e-5);

// parpe::OptimizationProblemImpl problem(std::move(hier),
// std::make_unique<parpe::Logger>());
// // std::vector<double> startingPoint = pReduced;
// // for(auto& pp : startingPoint)
// // pp += 1;
Expand All @@ -256,7 +279,8 @@ TEST_F(steadystateProblemTests, testSteadystateHierarchical) {
// //auto result = optimizer.optimize(&problem);
// // check status, cost, parameter
// //CHECK_EQUAL(1, std::get<0>(result));
// //DOUBLES_EQUAL(-parpe::getLogLikelihoodOffset(dp.edata[0].my.size()), std::get<1>(result), 1e-5);
// //DOUBLES_EQUAL(-parpe::getLogLikelihoodOffset(dp.edata[0].my.size()),
// std::get<1>(result), 1e-5);
// //DOUBLES_EQUAL(-1.0, std::get<2>(result).at(0), 1e-8);
// //std::cout<<std::get<2>(result);

Expand Down
Loading

0 comments on commit 8b4393d

Please sign in to comment.