From 165101c6375830242a916a605c5f4efe08af43d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vincent=20No=C3=ABl?= Date: Wed, 4 Sep 2024 10:12:40 -0300 Subject: [PATCH] Added first draft of custom pop output --- .../popmaboss/ICD_phenomenologicalPM.pbnd | 76 +++++++++++++++ .../ICD_phenomenological_TDC_ratio.cfg | 74 +++++++++++++++ engine/src/BooleanNetwork.h | 94 +++++++++++++++++++ engine/src/Cumulator.h | 4 + engine/src/PopMaBEstEngine.cc | 62 ++++++++---- engine/src/PopMaBEstEngine.h | 8 +- engine/src/PopMaBoSS.cc | 17 +++- engine/src/RunConfig.cc | 1 + engine/src/RunConfig.h | 6 +- engine/src/RunConfigGrammar.l | 3 + engine/src/RunConfigGrammar.y | 33 ++++++- engine/src/plot_dist.py | 25 +++++ 12 files changed, 381 insertions(+), 22 deletions(-) create mode 100644 engine/examples/popmaboss/ICD_phenomenologicalPM.pbnd create mode 100644 engine/examples/popmaboss/ICD_phenomenological_TDC_ratio.cfg create mode 100644 engine/src/plot_dist.py diff --git a/engine/examples/popmaboss/ICD_phenomenologicalPM.pbnd b/engine/examples/popmaboss/ICD_phenomenologicalPM.pbnd new file mode 100644 index 000000000..2c82a03e4 --- /dev/null +++ b/engine/examples/popmaboss/ICD_phenomenologicalPM.pbnd @@ -0,0 +1,76 @@ +node TumorCell +{ + rate_up = 0.0; + rate_down = 0.0; +} + +node DyingTumorCell +{ + rate_up = TumorCell ? $chemotherapy : 0.0 ; + rate_down = 0.0; +} + +node CALR +{ + rate_up = DyingTumorCell ? $u_CALR : 0.0; + // rate_down = Death ? $d_CALR : 0.0 ; + rate_down = 0.0; +} + +node ATP +{ + rate_up = DyingTumorCell ? $u_ATP : 0.0 ; + rate_down = DyingTumorCell ? 0.0 : $d_ATP ; +} + +node HMGB1 +{ + rate_up = DyingTumorCell ? $u_HMGB1 : 0.0 ; + rate_down = DyingTumorCell ? 0.0 : $d_HMGB1 ; +} + +node DC +{ + rate_up = 0.0; + rate_down = 0.0; +} + +node ActDC +{ + rate_up = DC ? ((#cell(CALR)/#cell(1)) * (#cell(ATP)/#cell(1)) *(#cell(HMGB1)/#cell(1))) : 0.0 ; + rate_down = 0.0; +} + +node MigrDC +{ + rate_up = ActDC ? $migrDC : 0.0; + rate_down = 0.0; +} + +node LNodeDC +{ + rate_up = MigrDC ? $lNodeDC : 0.0; + rate_down = 0.0; +} + +node TCell +{ + rate_up = LNodeDC ? $primingTCell : 0.0 ; + rate_down = 0.0 ; +} + +node CTL +{ + rate_up = TCell ? $activateCTL : 0.0 ; + rate_down = 0.0 ; +} + +death +{ + rate = (TumorCell ? (-log(1- (#cell(CTL)/#cell(1)))) : 0.0) + (CTL ? (-log(1- (#cell(TumorCell)/#cell(1)))) : 0.0) + (DyingTumorCell ? $deathChemotherapy : 0.0); +} + +division +{ + rate = (TumorCell ? $growthFactor : 0.0) + ((TCell & !(CTL)) ? (#cell(TumorCell)/#cell(1)*$clonalExpansion) : 0.0) ; + } diff --git a/engine/examples/popmaboss/ICD_phenomenological_TDC_ratio.cfg b/engine/examples/popmaboss/ICD_phenomenological_TDC_ratio.cfg new file mode 100644 index 000000000..1bbae5f00 --- /dev/null +++ b/engine/examples/popmaboss/ICD_phenomenological_TDC_ratio.cfg @@ -0,0 +1,74 @@ +// +// MaBoSS 2.0 configuration template generated at Wed Nov 20 11:08:24 2019 +// + +// global configuration variables +time_tick = 1; +//max_time = 720; +max_time = 720; +//sample_count = 100000; +sample_count = 160000; +discrete_time = 0; +use_physrandgen = 0; +seed_pseudorandom = 0; +display_traj = 0; +statdist_traj_count = 0; +statdist_cluster_threshold = 1; +thread_count = 16; +statdist_similarity_cache_max_size = 20000; +// init_pop = 20; + +// variables to be set in the configuration file or by using the --config-vars option +$u_CALR = 1/4; +$d_CALR = 50; +$u_ATP = 1/4; +$d_ATP = 50; +$u_HMGB1 = 1/6; +$d_HMGB1 = 50; +$migrDC = 1/(2*24); +$lNodeDC = 1/(3*24); +$activateCTL = 1/24; +$primingTCell = 1/(24); +$chemotherapy = .5; +$deathChemotherapy = 1/100; +$clonalExpansion = 1/12; +// $expCALR = 0; // define in upp +$growthFactor = 1/(240); +// $secrATP = 0; // define in upp +// $secrHMGB1 = 0; // defined in upp +$popCTL = 0; // defined in upp +// $popTumorCell = 0; // defined in upp +// set is_internal attribute value to 1 if node is an internal node +TumorCell.is_internal = 1; +DyingTumorCell.is_internal = 1; +CALR.is_internal = 1; +ATP.is_internal = 1; +HMGB1.is_internal = 1; +DC.is_internal = 0; +ActDC.is_internal = 1; +MigrDC.is_internal = 1; +LNodeDC.is_internal = 1; +TCell.is_internal = 0; +CTL.is_internal = 1; +// if NODE initial state is: +// - equals to 1: NODE.istate = 1; +// - equals to 0: NODE.istate = 0; +// - random: NODE.istate = -1; OR [NODE].istate = 0.5 [0], 0.5 [1]; OR skip NODE.istate declaration +// - weighted random: [NODE].istate = P0 [0], P1 [1]; where P0 and P1 are arithmetic expressions + +[TumorCell,DC].pop_istate = 1.0 [{[1,0]:18} , {[0,1]:2}]; + +TumorCell.istate = 0; +DC.istate = 0; +DyingTumorCell.istate = 0; +CALR.istate = 0; +ATP.istate = 0; +HMGB1.istate = 0; +ActDC.istate = 0; +MigrDC.istate = 0; +LNodeDC.istate = 0; +TCell.istate = 0; +CTL.istate = 0; + +custom_pop_output = #cell(DC) > 0 ? (#cell(DC -- TCell)+#cell(TCell))/#cell(DC) : 0.0; +// custom_pop_output = #cell(TCell); \ No newline at end of file diff --git a/engine/src/BooleanNetwork.h b/engine/src/BooleanNetwork.h index 5a2ef0d9f..ca02e52df 100644 --- a/engine/src/BooleanNetwork.h +++ b/engine/src/BooleanNetwork.h @@ -1331,6 +1331,64 @@ namespace std { }; } +class PopSize { + unsigned int size; +public: + PopSize(unsigned int size) : size(size) { } + PopSize() : size(0) { } + void set() {size = 0;} + unsigned int getSize() const {return size;} + static bool isPopState() {return false;} + + // & operator for applying the mask + PopSize operator&(const NetworkState_Impl& mask) const { + return PopSize(size); + } + + // & operator for applying the mask + PopSize operator&(const NetworkState& mask) const { + return PopSize(size); + } + + PopSize applyMask(const PopSize& mask, std::map& scale) const { + return PopSize(size); + } + + bool operator==(const PopSize& pop_size) const { + return pop_size.getSize() == size; + } + + bool operator<(const PopSize& pop_size) const { + return size < pop_size.getSize(); + } + + int hamming(Network* network, const NetworkState& state) const { + return 0; + } + std::set* getNetworkStates() const { + return new std::set(); + } + + std::string getName(Network * network, const std::string& sep=" -- ") const { + return std::to_string(size); + } + + void displayOneLine(std::ostream& os, Network* network, const std::string& sep = " -- ") const { + os << getName(network, sep); + } + + +}; + +namespace std { + template <> struct hash + { + size_t operator()(const PopSize & x) const + { + return x.getSize(); + } + }; +} // abstract base class used for expression evaluation class Expression { @@ -1426,6 +1484,42 @@ class NodeExpression : public Expression { } }; + +class StateExpression: public Expression { + NetworkState state; + Network* network; +public: + StateExpression(NetworkState state, Network* network) : state(state), network(network) { } + + Expression* clone() const {return new StateExpression(state, network);} + + double eval(const Node* this_node, const NetworkState& network_state) const { + return state.getState() == network_state.getState() ? 1.0 : 0.0; + } + + double eval(const Node* this_node, const NetworkState& network_state, const PopNetworkState& pop_state) const { + return state.getState() == network_state.getState() ? 1.0 : 0.0; + } + + bool hasCycle(Node* node) const { + return false; + } + + void display(std::ostream& os) const { + state.displayOneLine(os, network); + } + + bool isLogicalExpression() const {return true;} + std::vector getNodes() const{ + return std::vector(); + } + + void generateLogicalExpression(LogicalExprGenContext& genctx) const {} + + ~StateExpression() { + } +}; + class PopExpression : public Expression { Expression* expr; diff --git a/engine/src/Cumulator.h b/engine/src/Cumulator.h index 72865a709..12b556a2a 100644 --- a/engine/src/Cumulator.h +++ b/engine/src/Cumulator.h @@ -694,6 +694,10 @@ class Cumulator { this->output_mask = output_mask; } + S getOutputMask() const { + return output_mask; + } + void setRefnodeMask(const NetworkState_Impl& refnode_mask) { this->refnode_mask = refnode_mask; } diff --git a/engine/src/PopMaBEstEngine.cc b/engine/src/PopMaBEstEngine.cc index 185a8a179..d4c188676 100644 --- a/engine/src/PopMaBEstEngine.cc +++ b/engine/src/PopMaBEstEngine.cc @@ -46,13 +46,15 @@ */ #include "PopMaBEstEngine.h" +#include "BooleanNetwork.h" +#include "Cumulator.h" #include "MetaEngine.h" +#include "ProbTrajDisplayer.h" #include "Probe.h" #include #include #include #include -#include "Utils.h" #ifndef WINDOWS #include #else @@ -105,7 +107,9 @@ PopMaBEstEngine::PopMaBEstEngine(PopNetwork *pop_network, RunConfig *runconfig) } merged_cumulator = NULL; - cumulator_v.resize(thread_count); + custom_pop_cumulator = NULL; + custom_pop_cumulator_v.resize(thread_count, NULL); + cumulator_v.resize(thread_count, NULL); unsigned int count = sample_count / thread_count; unsigned int firstcount = count + sample_count - count * thread_count; @@ -114,15 +118,17 @@ PopMaBEstEngine::PopMaBEstEngine(PopNetwork *pop_network, RunConfig *runconfig) for (unsigned int nn = 0; nn < thread_count; ++nn) { + Cumulator *custom_cumulator = new Cumulator(runconfig, runconfig->getTimeTick(), runconfig->getMaxTime(), (nn == 0 ? firstcount : count), (nn == 0 ? first_scount : scount)); + custom_pop_cumulator_v[nn] = custom_cumulator; Cumulator *cumulator = new Cumulator(runconfig, runconfig->getTimeTick(), runconfig->getMaxTime(), (nn == 0 ? firstcount : count), (nn == 0 ? first_scount : scount )); if (has_internal) { -#ifdef USE_STATIC_BITSET + #ifdef USE_STATIC_BITSET NetworkState_Impl state2 = ~internal_state.getState(); cumulator->setOutputMask(PopNetworkState(state2, 1)); -#else + #else cumulator->setOutputMask(PopNetworkState(~internal_state.getState(), 1)); -#endif + #endif } cumulator_v[nn] = cumulator; } @@ -194,13 +200,15 @@ struct ArgWrapper unsigned int start_count_thread; unsigned int sample_count_thread; Cumulator *cumulator; + Cumulator* custom_cumulator; RandomGeneratorFactory *randgen_factory; long long int* elapsed_time; int seed; FixedPoints *fixpoint_map; + PopExpression* custom_pop_expression; std::ostream *output_traj; - ArgWrapper(PopMaBEstEngine *mabest, unsigned int start_count_thread, unsigned int sample_count_thread, Cumulator *cumulator, RandomGeneratorFactory *randgen_factory, long long int * elapsed_time, int seed, FixedPoints *fixpoint_map, std::ostream *output_traj) : mabest(mabest), start_count_thread(start_count_thread), sample_count_thread(sample_count_thread), cumulator(cumulator), randgen_factory(randgen_factory), elapsed_time(elapsed_time), seed(seed), fixpoint_map(fixpoint_map), output_traj(output_traj) {} + ArgWrapper(PopMaBEstEngine *mabest, unsigned int start_count_thread, unsigned int sample_count_thread, Cumulator *cumulator, Cumulator* custom_cumulator, RandomGeneratorFactory *randgen_factory, long long int * elapsed_time, int seed, FixedPoints *fixpoint_map, PopExpression* custom_pop_expression, std::ostream *output_traj) : mabest(mabest), start_count_thread(start_count_thread), sample_count_thread(sample_count_thread), cumulator(cumulator), custom_cumulator(custom_cumulator), randgen_factory(randgen_factory), elapsed_time(elapsed_time), seed(seed), fixpoint_map(fixpoint_map), custom_pop_expression(custom_pop_expression), output_traj(output_traj) {} }; void *PopMaBEstEngine::threadWrapper(void *arg) @@ -211,7 +219,7 @@ void *PopMaBEstEngine::threadWrapper(void *arg) ArgWrapper *warg = (ArgWrapper *)arg; try { - warg->mabest->runThread(warg->cumulator, warg->start_count_thread, warg->sample_count_thread, warg->randgen_factory, warg->seed, warg->fixpoint_map, warg->output_traj); + warg->mabest->runThread(warg->cumulator, warg->custom_cumulator, warg->start_count_thread, warg->sample_count_thread, warg->randgen_factory, warg->seed, warg->fixpoint_map, warg->custom_pop_expression, warg->output_traj); } catch (const BNException &e) { @@ -223,7 +231,7 @@ void *PopMaBEstEngine::threadWrapper(void *arg) return NULL; } -void PopMaBEstEngine::runThread(Cumulator *cumulator, unsigned int start_count_thread, unsigned int sample_count_thread, RandomGeneratorFactory *randgen_factory, int seed, FixedPoints *fixpoint_map, std::ostream *output_traj) +void PopMaBEstEngine::runThread(Cumulator *cumulator, Cumulator* custom_cumulator, unsigned int start_count_thread, unsigned int sample_count_thread, RandomGeneratorFactory *randgen_factory, int seed, FixedPoints *fixpoint_map, PopExpression* custom_pop_expression, std::ostream *output_traj) { const std::vector &nodes = pop_network->getNodes(); PopNetworkState pop_network_state; @@ -233,6 +241,8 @@ void PopMaBEstEngine::runThread(Cumulator *cumulator, unsigned { random_generator->setSeed(seed + start_count_thread + nn); cumulator->rewind(); + custom_cumulator->rewind(); + if (verbose > 1) std::cout << "> New simulation" << std::endl; @@ -409,8 +419,17 @@ void PopMaBEstEngine::runThread(Cumulator *cumulator, unsigned (*output_traj) << '\t' << TH << '\n'; } - cumulator->cumul(pop_network_state, tm, TH); - + if (runconfig->hasCustomPopOutput()) + { + NetworkState s; + PopNetworkState t_popstate(pop_network_state & cumulator->getOutputMask().getMap().begin()->first); + PopSize pop_size((unsigned int) custom_pop_expression->eval(NULL, s, t_popstate)); + // std::cout << "PopNetworkState = " << t_popstate.getName(network) << ", PopSize = " << pop_size.getSize() << std::endl; + custom_cumulator->cumul(pop_size, tm, TH); + } else { + cumulator->cumul(pop_network_state, tm, TH); + } + if (tm >= max_time) { break; @@ -432,6 +451,8 @@ void PopMaBEstEngine::runThread(Cumulator *cumulator, unsigned std::cout << std::endl; cumulator->trajectoryEpilogue(); + custom_cumulator->trajectoryEpilogue(); + } delete random_generator; } @@ -461,9 +482,9 @@ void PopMaBEstEngine::run(std::ostream *output_traj) fixpoint_map_v.push_back(fixpoint_map); #ifdef MPI_COMPAT - ArgWrapper* warg = new ArgWrapper(this, start_sample_count, cumulator_v[nn]->getSampleCount(), cumulator_v[nn], randgen_factory, &(thread_elapsed_runtimes[world_rank][nn]), seed, fixpoint_map, output_traj); + ArgWrapper* warg = new ArgWrapper(this, start_sample_count, cumulator_v[nn]->getSampleCount(), cumulator_v[nn], custom_pop_cumulator_v[nn], randgen_factory, &(thread_elapsed_runtimes[world_rank][nn]), seed, fixpoint_map, runconfig->getCustomPopOutputExpression(), output_traj); #else - ArgWrapper* warg = new ArgWrapper(this, start_sample_count, cumulator_v[nn]->getSampleCount(), cumulator_v[nn], randgen_factory, &(thread_elapsed_runtimes[nn]), seed, fixpoint_map, output_traj); + ArgWrapper* warg = new ArgWrapper(this, start_sample_count, cumulator_v[nn]->getSampleCount(), cumulator_v[nn], custom_pop_cumulator_v[nn], randgen_factory, &(thread_elapsed_runtimes[nn]), seed, fixpoint_map, runconfig->getCustomPopOutputExpression(), output_traj); #endif pthread_create(&tid[nn], NULL, PopMaBEstEngine::threadWrapper, warg); @@ -523,11 +544,13 @@ void PopMaBEstEngine::mergePairOfFixpoints(FixedPoints* fixpoints_1, FixedPoints struct PopMergeWrapper { Cumulator* cumulator_1; Cumulator* cumulator_2; + Cumulator* custom_cumulator_1; + Cumulator* custom_cumulator_2; FixedPoints* fixpoints_1; FixedPoints* fixpoints_2; - PopMergeWrapper(Cumulator* cumulator_1, Cumulator* cumulator_2, FixedPoints* fixpoints_1, FixedPoints* fixpoints_2) : - cumulator_1(cumulator_1), cumulator_2(cumulator_2), fixpoints_1(fixpoints_1), fixpoints_2(fixpoints_2) { } + PopMergeWrapper(Cumulator* cumulator_1, Cumulator* cumulator_2, Cumulator* custom_cumulator_1, Cumulator* custom_cumulator_2, FixedPoints* fixpoints_1, FixedPoints* fixpoints_2) : + cumulator_1(cumulator_1), cumulator_2(cumulator_2), custom_cumulator_1(custom_cumulator_1), custom_cumulator_2(custom_cumulator_2), fixpoints_1(fixpoints_1), fixpoints_2(fixpoints_2) { } }; void* PopMaBEstEngine::threadMergeWrapper(void *arg) @@ -538,6 +561,7 @@ void* PopMaBEstEngine::threadMergeWrapper(void *arg) PopMergeWrapper* warg = (PopMergeWrapper*)arg; try { Cumulator::mergePairOfCumulators(warg->cumulator_1, warg->cumulator_2); + Cumulator::mergePairOfCumulators(warg->custom_cumulator_1, warg->custom_cumulator_2); PopMaBEstEngine::mergePairOfFixpoints(warg->fixpoints_1, warg->fixpoints_2); } catch(const BNException& e) { std::cerr << e; @@ -569,7 +593,7 @@ void PopMaBEstEngine::mergeResults(std::vector*> cumu for(unsigned int i=0; i < size; i+=(step_lvl*2)) { if (i+step_lvl < size) { - PopMergeWrapper* warg = new PopMergeWrapper(cumulator_v[i], cumulator_v[i+step_lvl], fixpoint_map_v[i], fixpoint_map_v[i+step_lvl]); + PopMergeWrapper* warg = new PopMergeWrapper(cumulator_v[i], cumulator_v[i+step_lvl], custom_pop_cumulator_v[i], custom_pop_cumulator_v[i+step_lvl], fixpoint_map_v[i], fixpoint_map_v[i+step_lvl]); pthread_create(&tid[nb_threads], NULL, PopMaBEstEngine::threadMergeWrapper, warg); nb_threads++; wargs.push_back(warg); @@ -751,6 +775,7 @@ void PopMaBEstEngine::epilogue() mergeResults(cumulator_v, fixpoint_map_v); merged_cumulator = cumulator_v[0]; + custom_pop_cumulator = custom_pop_cumulator_v[0]; fixpoints = fixpoint_map_v[0]; #ifdef MPI_COMPAT @@ -761,7 +786,7 @@ void PopMaBEstEngine::epilogue() { #endif merged_cumulator->epilogue(pop_network, reference_state); - + custom_pop_cumulator->epilogue(pop_network, reference_state); #ifdef MPI_COMPAT } #endif @@ -776,6 +801,7 @@ PopMaBEstEngine::~PopMaBEstEngine() delete t_arg_wrapper; delete merged_cumulator; + delete custom_pop_cumulator; } void PopMaBEstEngine::displayFixpoints(FixedPointDisplayer *displayer) const @@ -854,6 +880,10 @@ if (getWorldRank() == 0) { #endif } +void PopMaBEstEngine::displayCustomPopProbTraj(ProbTrajDisplayer* displayer) const +{ + custom_pop_cumulator->displayProbTraj(pop_network, refnode_count, displayer); +} const std::map > PopMaBEstEngine::getFixPointsDists() const { std::map > res; diff --git a/engine/src/PopMaBEstEngine.h b/engine/src/PopMaBEstEngine.h index 6f0c62c84..089a6cad8 100644 --- a/engine/src/PopMaBEstEngine.h +++ b/engine/src/PopMaBEstEngine.h @@ -53,6 +53,7 @@ #include #include +#include "CustomPopProbTrajDisplayer.h" #include "MetaEngine.h" #include "FixedPointEngine.h" #include "BooleanNetwork.h" @@ -80,7 +81,8 @@ class PopMaBEstEngine : public MetaEngine { Cumulator* merged_cumulator; std::vector* > cumulator_v; - + Cumulator* custom_pop_cumulator; + std::vector* > custom_pop_cumulator_v; public: static const std::string VERSION; static int verbose; @@ -116,9 +118,9 @@ class PopMaBEstEngine : public MetaEngine { double computeTH(const MAP& nodeTransitionRates, double total_rate) const; void epilogue(); static void* threadWrapper(void *arg); - void runThread(Cumulator* cumulator, unsigned int start_count_thread, unsigned int sample_count_thread, RandomGeneratorFactory* randgen_factory, int seed, FixedPoints* fixpoint_map, std::ostream* output_traj); + void runThread(Cumulator* cumulator, Cumulator* custom_cumulator, unsigned int start_count_thread, unsigned int sample_count_thread, RandomGeneratorFactory* randgen_factory, int seed, FixedPoints* fixpoint_map, PopExpression* simple_pop_expression, std::ostream* output_traj); void displayRunStats(std::ostream& os, time_t start_time, time_t end_time) const; - + void displayCustomPopProbTraj(ProbTrajDisplayer* displayer) const; static void mergePairOfFixpoints(FixedPoints* fixpoints_1, FixedPoints* fixpoints_2); static void* threadMergeWrapper(void *arg); void mergeResults(std::vector*> cumulator_v, std::vector fixpoint_map_v); diff --git a/engine/src/PopMaBoSS.cc b/engine/src/PopMaBoSS.cc index eaff6e835..70317b6b5 100644 --- a/engine/src/PopMaBoSS.cc +++ b/engine/src/PopMaBoSS.cc @@ -46,9 +46,12 @@ */ #include +#include "BooleanNetwork.h" +#include "CustomPopProbTrajDisplayer.h" #include "PopMaBEstEngine.h" #include "Function.h" #include +#include #include #include "Utils.h" #include "RandomGenerator.h" @@ -333,6 +336,7 @@ int main(int argc, char* argv[]) std::ostream* output_fp = NULL; std::ostream* output_pop_probtraj = NULL; std::ostream* output_simple_pop_probtraj = NULL; + std::ostream* output_custom_pop_probtraj = NULL; #ifdef HDF5_COMPAT hid_t hdf5_file; @@ -412,6 +416,9 @@ int main(int argc, char* argv[]) output_fp = new std::ofstream((std::string(output) + "_fp" + format_extension(format)).c_str()); output_pop_probtraj = new std::ofstream((std::string(output) + "_pop_probtraj" + format_extension(format)).c_str()); output_simple_pop_probtraj = new std::ofstream((std::string(output) + "_simple_pop_probtraj" + format_extension(format)).c_str()); + if (runconfig->hasCustomPopOutput()) { + output_custom_pop_probtraj = new std::ofstream((std::string(output) + "_custom_pop_probtraj" + format_extension(format)).c_str()); + } #ifdef HDF5_COMPAT } else if (format == HDF5_FORMAT) { hdf5_file = H5Fcreate((std::string(output) + format_extension(format)).c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); @@ -431,10 +438,15 @@ int main(int argc, char* argv[]) ProbTrajDisplayer* pop_probtraj_displayer; FixedPointDisplayer* fp_displayer; + ProbTrajDisplayer* custom_pop_probtraj_displayer = NULL; if (format == CSV_FORMAT) { pop_probtraj_displayer = new CSVSimplePopProbTrajDisplayer(pop_network, *output_pop_probtraj, *output_simple_pop_probtraj, hexfloat); fp_displayer = new CSVFixedPointDisplayer(pop_network, *output_fp, hexfloat); + if (runconfig->hasCustomPopOutput()){ + custom_pop_probtraj_displayer = new CSVCustomPopProbTrajDisplayer(pop_network, *output_custom_pop_probtraj, hexfloat); + } + } else if (format == JSON_FORMAT) { pop_probtraj_displayer = new JSONSimpleProbTrajDisplayer(pop_network, *output_pop_probtraj, *output_simple_pop_probtraj, hexfloat); // Use CSV displayer for fixed points as the Json one is not fully implemented @@ -448,9 +460,12 @@ int main(int argc, char* argv[]) pop_probtraj_displayer = NULL; fp_displayer = NULL; } - + mabest.display(pop_probtraj_displayer, fp_displayer); + if (runconfig->hasCustomPopOutput()){ + mabest.displayCustomPopProbTraj(custom_pop_probtraj_displayer); + } time(&end_time); mabest.displayRunStats(*output_run, start_time, end_time); diff --git a/engine/src/RunConfig.cc b/engine/src/RunConfig.cc index cc456b1c6..d2e56d1c2 100644 --- a/engine/src/RunConfig.cc +++ b/engine/src/RunConfig.cc @@ -77,6 +77,7 @@ RunConfig::RunConfig() statdist_similarity_cache_max_size = 20000; init_pop = 1; pop_base = 1.0; + custom_pop_output_expression = NULL; } RunConfig::~RunConfig() diff --git a/engine/src/RunConfig.h b/engine/src/RunConfig.h index 2b5b2cd28..dac88e6d1 100644 --- a/engine/src/RunConfig.h +++ b/engine/src/RunConfig.h @@ -75,6 +75,8 @@ class RunConfig { void dump_perform(Network* network, std::ostream& os, bool is_template, std::string version) const; unsigned int init_pop; double pop_base; + PopExpression* custom_pop_output_expression; + public: RunConfig(); ~RunConfig(); @@ -82,7 +84,9 @@ class RunConfig { int parse(Network* network, const char* file = NULL); int parseExpression(Network* network, const char* expr); void setParameter(const std::string& param, double value); - + void setCustomPopOutputExpression(Expression* expr) { custom_pop_output_expression = (PopExpression*) expr; } + bool hasCustomPopOutput() const { return custom_pop_output_expression != NULL; } + PopExpression* getCustomPopOutputExpression() { return custom_pop_output_expression; } RandomGeneratorFactory* getRandomGeneratorFactory() const; unsigned int getInitPop() const {return init_pop;} double getPopBase() const {return pop_base;} diff --git a/engine/src/RunConfigGrammar.l b/engine/src/RunConfigGrammar.l index 2715ed815..5f8083134 100644 --- a/engine/src/RunConfigGrammar.l +++ b/engine/src/RunConfigGrammar.l @@ -69,6 +69,9 @@ exppart [eE](\-|\+)?[0-9]+ "//".* { } [Tt][Rr][Uu][Ee] { yylval.l = 1; return INTEGER; } [Ff][Aa][Ll][Ss][Ee] { yylval.l = 0; return INTEGER; } +"#cell" { return NUMBER_CELL; } +"custom_pop_output" { return CUSTOM_POP_OUTPUT; } +" -- " { return NODE_SEP; } "$"[a-zA-Z_][a-zA-Z_0-9]* { yylval.str = strdup(yytext); return VARIABLE; diff --git a/engine/src/RunConfigGrammar.y b/engine/src/RunConfigGrammar.y index eeb7bced5..070f02c09 100644 --- a/engine/src/RunConfigGrammar.y +++ b/engine/src/RunConfigGrammar.y @@ -68,6 +68,7 @@ extern std::string yy_error_head(); double d; long long l; std::vector* node_list; + const Node* node; std::vector* expr_list; IStateGroup::ProbaIState* istate_expr; PopIStateGroup::PopProbaIState* pop_istate_expr; @@ -93,6 +94,7 @@ extern std::string yy_error_head(); %type expression %type symbol_list %type symbol_istate_list +%type state %type expression_list %type istate_expression %type istate_expression_list @@ -107,7 +109,7 @@ extern std::string yy_error_head(); %token DOUBLE %token INTEGER -%token LOGAND LOGOR LOGXOR LOGNOT EQUAL NOT_EQUAL NODE GTEQ LTEQ +%token LOGAND LOGOR LOGXOR LOGNOT EQUAL NOT_EQUAL GTEQ LTEQ CUSTOM_POP_OUTPUT NUMBER_CELL NODE_SEP %% @@ -122,6 +124,7 @@ translation_unit: decl decl: var_decl | runconfig_decl | node_attr_decl +| custom_pop_output_decl | ';' ; @@ -303,6 +306,12 @@ var_decl: VARIABLE '=' expression ';' } ; +custom_pop_output_decl: CUSTOM_POP_OUTPUT '=' expression ';' +{ + config->setCustomPopOutputExpression($3); +} +; + primary_expression: INTEGER { $$ = new ConstantExpression($1); @@ -320,8 +329,30 @@ primary_expression: INTEGER { $$ = new ParenthesisExpression($2); } +| NUMBER_CELL '(' state ')' +{ + NetworkState network_state; + for (auto* node : *$3) { + network_state.setNodeState(node, 1); + } + StateExpression* state_expr = new StateExpression(network_state, network); + + $$ = new PopExpression(state_expr); +} ; +state: SYMBOL +{ + $$ = new std::vector(); + $$->push_back(network->getNode($1)); +} +| state NODE_SEP SYMBOL +{ + $$ = $1; + $$->push_back(network->getNode($3)); +} +; + argument_expression_list: expression { $$ = new ArgumentList(); diff --git a/engine/src/plot_dist.py b/engine/src/plot_dist.py new file mode 100644 index 000000000..916962f7b --- /dev/null +++ b/engine/src/plot_dist.py @@ -0,0 +1,25 @@ +import sys, os +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + + +if len(sys.argv) < 2: + print("Usage: plot_dist.py ") + exit(1) + +dist_file = sys.argv[1] +first_index = -1 +with open(dist_file, 'r') as f: + for i, line in enumerate(f): + if i == 0: + first_index = next(i for i, col in enumerate(line.strip("\n").split("\t")) if col == "State") + if line.startswith(sys.argv[2]): + pop_sizes = [int(x) for x in line.strip("\n").split("\t")[first_index::3]] + pop_probas = line.strip("\n").split("\t")[first_index+1::3] + res = pd.Series(pop_probas, index=pop_sizes).astype(float) + res.sort_index(inplace=True) + res.plot() + plt.show() + break +