Skip to content

Commit

Permalink
Implement biomass sampling experiments for metabolic networks in C++ (#…
Browse files Browse the repository at this point in the history
…197)

* implement biomass sampling experiments

* modify circleci/config.yml

* fix bug in log-concave tests

* change test function name
  • Loading branch information
TolisChal authored Feb 8, 2022
1 parent f9d3414 commit 6591ae8
Show file tree
Hide file tree
Showing 4 changed files with 263 additions and 32 deletions.
3 changes: 3 additions & 0 deletions test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,7 @@ else ()
include_directories (BEFORE ../include/lp_oracles)
include_directories (BEFORE ../include/nlp_oracles)
include_directories (BEFORE ../include/misc)
include_directories (BEFORE ../test)

#for Eigen
if (${CMAKE_VERSION} VERSION_LESS "3.12.0")
Expand Down Expand Up @@ -295,6 +296,8 @@ else ()
COMMAND logconcave_sampling_test -tc=hmc)
add_test(NAME logconcave_sampling_test_uld
COMMAND logconcave_sampling_test -tc=uld)
add_test(NAME logconcave_sampling_test_exponential_biomass_sampling
COMMAND logconcave_sampling_test -tc=exponential_biomass_sampling)

add_executable (simple_mc_integration simple_mc_integration.cpp $<TARGET_OBJECTS:test_main>)
add_test(NAME simple_mc_integration_over_limits
Expand Down
111 changes: 79 additions & 32 deletions test/logconcave_sampling_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -703,7 +703,7 @@ void benchmark_polytope_linear_program_optimization(
bool rounding=false,
bool centered=true,
unsigned int max_draws=80000,
unsigned int num_burns=20000) {
unsigned int num_burns=10000) {
typedef Cartesian<NT> Kernel;
typedef std::vector<Point> pts;
typedef boost::mt19937 RNGType;
Expand Down Expand Up @@ -756,6 +756,8 @@ void benchmark_polytope_linear_program_optimization(
GaussianRDHRWalk::Walk<Polytope, RandomNumberGenerator> gaussian_walk(P, x0, lp_params.L, rng);
int n_warmstart_samples = 100;

std::cout << "Warm start" << std::endl;

for (int i = 0; i < n_warmstart_samples; i++) {
gaussian_walk.apply(P, x0, lp_params.L, walk_length, rng);
}
Expand Down Expand Up @@ -788,14 +790,15 @@ void benchmark_polytope_linear_program_optimization(
hmc.apply(rng, walk_length);
}

hmc.solver->eta = hmc.solver->eta*NT(10);
hmc.disable_adaptive();

std::cout << std::endl;
std::cout << "Optimizing" << std::endl;
unsigned int num_phases = 0;
std::cout << "Sampling" << std::endl;
unsigned int num_phases = 1;

auto start = std::chrono::high_resolution_clock::now();
for (unsigned int j = 0; j < max_phases; j++) {
num_phases++;
std::cout << "Temperature " << opt_params.T << std::endl;
for (unsigned int j = 0; j < 1; j++) {
for (unsigned int i = 0; i < max_actual_draws; i++) {
hmc.apply(rng, walk_length);
if (f_lp(minimum) >= f_lp(hmc.x)) {
Expand All @@ -813,12 +816,12 @@ void benchmark_polytope_linear_program_optimization(

NT ETA = (NT) std::chrono::duration_cast<std::chrono::microseconds>(stop - start).count();

std::cerr << "Min LP Value: " << f_lp(minimum) << std::endl;
std::cerr << "Min biomass Value: " << f_lp(minimum) << std::endl;
std::cerr << "Argmin: " << minimum.getCoefficients().transpose() << std::endl;
std::cerr << "Average ETA: " << ETA / (NT) num_phases << std::endl;
std::cerr << "Average Time per Independent sample: " << ETA / total_min_ess << std::endl;
std::cerr << "Average Max PSRF: " << total_max_psrf / (NT) num_phases << std::endl;
std::cerr << "Average Min ESS: " << total_min_ess / (NT) num_phases << std::endl;
std::cerr << "ETA: " << ETA / (NT) num_phases << std::endl;
std::cerr << "Time per Independent sample: " << total_min_ess / ETA << std::endl;
std::cerr << "Max PSRF: " << total_max_psrf / (NT) num_phases << std::endl;
std::cerr << "Min ESS: " << total_min_ess / (NT) num_phases << std::endl;
std::cerr << "Average number of reflections: " <<
(1.0 * hmc.solver->num_reflections) / hmc.solver->num_steps << std::endl;
std::cerr << "Step size (final): " << hmc.solver->eta << std::endl;
Expand Down Expand Up @@ -893,10 +896,10 @@ void call_test_benchmark_polytopes_grid_search() {
std::make_tuple(generate_cube<Hpolytope>(100, false), "100_cube", false),
std::make_tuple(generate_prod_simplex<Hpolytope>(50, false), "50_prod_simplex", false),
std::make_tuple(generate_birkhoff<Hpolytope>(10), "10_birkhoff", false),
std::make_tuple(read_polytope<Hpolytope, NT>("./metabolic_full_dim/polytope_iAB_RBC_283.ine"), "iAB_RBC_283", true),
std::make_tuple(read_polytope<Hpolytope, NT>("./metabolic_full_dim/polytope_iAT_PLT_636.ine"), "iAT_PLT_636", true),
std::make_tuple(read_polytope<Hpolytope, NT>("./metabolic_full_dim/polytope_e_coli.ine"), "e_coli", true),
std::make_tuple(read_polytope<Hpolytope, NT>("./metabolic_full_dim/polytope_recon2.ine"), "recon2", true)
std::make_tuple(read_polytope<Hpolytope, NT>("metabolic_full_dim/polytope_iAB_RBC_283.ine"), "iAB_RBC_283", true),
std::make_tuple(read_polytope<Hpolytope, NT>("metabolic_full_dim/polytope_iAT_PLT_636.ine"), "iAT_PLT_636", true),
std::make_tuple(read_polytope<Hpolytope, NT>("metabolic_full_dim/polytope_e_coli.ine"), "e_coli", true),
std::make_tuple(read_polytope<Hpolytope, NT>("metabolic_full_dim/polytope_recon2.ine"), "recon2", true)
};

Hpolytope P;
Expand Down Expand Up @@ -998,30 +1001,74 @@ void call_test_benchmark_spectrahedra_grid_search() {

}

template <typename Point, typename NT>
Point load_biomass_function(std::string name)
{
std::vector<double> v1;
std::string line;
NT element;

std::ifstream myFile(name);
while(std::getline(myFile, line))
{
std::istringstream lineStream(line);

while(lineStream >> element) {
v1.push_back(element);
}
}
Point biomass_function = Point(v1.size(), v1);
NT length = biomass_function.length();
biomass_function *= (NT(1) / length);
return biomass_function;
}

inline bool exists_check (const std::string& name) {
std::ifstream f(name.c_str());
return f.good();
}

template <typename NT>
void call_test_optimization() {
void call_test_exp_sampling() {
typedef Cartesian<NT> Kernel;
typedef typename Kernel::Point Point;
typedef HPolytope<Point> Hpolytope;
std::string name;
std::vector<std::tuple<Hpolytope, Point, std::string, bool>> polytopes;


if (exists_check("metabolic_full_dim/e_coli_biomass_function.txt") && exists_check("metabolic_full_dim/polytope_e_coli.ine")){
Point biomass_function_e_coli = load_biomass_function<Point, NT>("metabolic_full_dim/e_coli_biomass_function.txt");
polytopes.push_back(std::make_tuple(read_polytope<Hpolytope, NT>("metabolic_full_dim/polytope_e_coli.ine"), biomass_function_e_coli, "e_coli", true));
}

std::vector<std::tuple<Hpolytope, std::string, bool>> polytopes{
std::make_tuple(generate_cube<Hpolytope>(100, false), "100_cube", false),
std::make_tuple(read_polytope<Hpolytope, NT>("./metabolic_full_dim/polytope_e_coli.ine"), "e_coli", true),
};
if (exists_check("metabolic_full_dim/iAT_PTL_636_biomass_function.txt") && exists_check("metabolic_full_dim/polytope_iAT_PTL_636.ine")){
Point biomass_function_iAT = load_biomass_function<Point, NT>("metabolic_full_dim/iAT_PTL_636_biomass_function.txt");
polytopes.push_back(std::make_tuple(read_polytope<Hpolytope, NT>("metabolic_full_dim/polytope_iAT_PTL_636.ine"), biomass_function_iAT, "iAT_PTL_636", true));
}

if (exists_check("metabolic_full_dim/recon1_function.txt") && exists_check("metabolic_full_dim/polytope_recon1.ine")){
Point biomass_function_recon1 = load_biomass_function<Point, NT>("metabolic_full_dim/recon1_biomass_function.txt");
polytopes.push_back(std::make_tuple(read_polytope<Hpolytope, NT>("metabolic_full_dim/polytope_recon1.ine"), biomass_function_recon1, "recon1", true));
}

Hpolytope P;
std::string name;
std::ofstream outfile;
typedef BoostRandomNumberGenerator<boost::mt19937, NT> RNGType;

for (std::tuple<Hpolytope, std::string, bool> polytope_tuple : polytopes) {
P = std::get<0>(polytope_tuple);
name = std::get<1>(polytope_tuple);
std::cerr << name << std::endl;
P.normalize();
Point coeffs = Point::all_ones(P.dimension());
for (unsigned int walk_length = 500; walk_length <= P.dimension(); walk_length += P.dimension() / 10) {
benchmark_polytope_linear_program_optimization<NT, Hpolytope>(coeffs, P, 0, NT(-1), walk_length, false);
}
for (std::tuple<Hpolytope, Point, std::string, bool> polytope_tuple : polytopes) {
P = std::get<0>(polytope_tuple);
name = std::get<2>(polytope_tuple);

std::cerr << "Model: " + name << std::endl;
std::cout<< "Dimension of " + name + ": " <<P.dimension()<<std::endl;

P.normalize();
RNGType rng(P.dimension());
Point coeffs = std::get<1>(polytope_tuple);
for (unsigned int walk_length = P.dimension(); walk_length <= 2*P.dimension(); walk_length += P.dimension()) {
benchmark_polytope_linear_program_optimization<NT, Hpolytope>(coeffs, P, 0, NT(-1), walk_length, false);
}
}

}
Expand Down Expand Up @@ -1145,8 +1192,8 @@ TEST_CASE("uld") {
call_test_uld<double>();
}

TEST_CASE("optimization") {
call_test_optimization<double>();
TEST_CASE("exponential_biomass_sampling") {
call_test_exp_sampling<double>();
}

TEST_CASE("benchmark_hmc") {
Expand Down
1 change: 1 addition & 0 deletions test/metabolic_full_dim/e_coli_biomass_function.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
1.057595362706991726e-01 6.268231060406626726e-04 5.417130007052790575e-01 -4.161378433635810614e-01 -2.634451558388141157e-01 -1.377859647362443540e-01 -1.239169700599200043e-02 -1.771478818636848418e+00 -1.474689492379328559e+00 3.575169148021269949e-01 -1.913772250344042547e+01 1.261758376970891285e-02 -1.425712605922182963e+00 2.234425519598324783e+01 2.942928129761246510e+02 8.469069143284598056e+01 1.059563561026396172e+01 -7.351985582038952316e+01 1.377496801177623098e+01 1.319707473523769181e+02 -3.196840752498849625e+01 -7.523019962729542875e+00 -1.761179600696947745e+01 -6.982659339018495892e+00
Loading

0 comments on commit 6591ae8

Please sign in to comment.