diff --git a/CMakeLists.txt b/CMakeLists.txt index 8e3c745..61bb9c3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -615,11 +615,16 @@ if (WCS_HAS_ROSS) $ $ $) - target_link_libraries(des-bin PRIVATE ROSS::ROSS MPI::MPI_CXX wcs ${LIB_FILESYSTEM}) + target_link_libraries(des-bin PRIVATE wcs ROSS::ROSS MPI::MPI_CXX ${LIB_FILESYSTEM}) set_target_properties(des-bin PROPERTIES OUTPUT_NAME des) list(APPEND WCS_EXEC_TARGETS des-bin) endif (WCS_HAS_ROSS) +# add executable hybrid +if (WCS_HAS_ROSS) + list(APPEND WCS_EXEC_TARGETS hybrid-bin) +endif (WCS_HAS_ROSS) + # Build tests add_executable( t_state_rngen-bin src/utils/unit_tests/t_state_rngen.cpp ) target_include_directories(t_state_rngen-bin PUBLIC @@ -713,6 +718,7 @@ set(INCLUDE_INSTALL_DIRS "${CMAKE_SOURCE_DIR}/src/proto" "${CMAKE_SOURCE_DIR}/src/reaction_network" "${CMAKE_SOURCE_DIR}/src/sim_methods" + "${CMAKE_SOURCE_DIR}/src/hybrid" "${CMAKE_SOURCE_DIR}/src/utils") set(LIB_INSTALL_DIR "${CMAKE_BINARY_DIR}") set(EXTRA_CMAKE_MODULE_DIR "${CMAKE_SOURCE_DIR}/cmake/modules") @@ -787,11 +793,16 @@ install( DIRECTORY "${PROJECT_SOURCE_DIR}/src/utils" DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}" FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp") +install( + DIRECTORY "${PROJECT_SOURCE_DIR}/src/hybrid" + DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}" + FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp" PATTERN "*.hh") install( FILES "${PROJECT_BINARY_DIR}/wcs_config.hpp" "${PROJECT_SOURCE_DIR}/src/wcs_types.hpp" "${PROJECT_SOURCE_DIR}/src/bgl.hpp" "${PROJECT_SOURCE_DIR}/src/des.hpp" + "${PROJECT_SOURCE_DIR}/src/hybrid/hybrid.hpp" "${PROJECT_SOURCE_DIR}/src/wcs-ross-bf.hpp" DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}") diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c55fc58..e55c471 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -16,6 +16,7 @@ add_subdirectory(proto) add_subdirectory(reaction_network) add_subdirectory(sim_methods) add_subdirectory(utils) +add_subdirectory(hybrid) # Propagate the files up the tree set(WCS_HEADERS "${WCS_HEADERS}" "${THIS_DIR_HEADERS}" PARENT_SCOPE) diff --git a/src/hybrid/CMakeLists.txt b/src/hybrid/CMakeLists.txt new file mode 100644 index 0000000..51c3837 --- /dev/null +++ b/src/hybrid/CMakeLists.txt @@ -0,0 +1,58 @@ +# blt_add_executable( +# NAME hybrid +# SOURCES actor.cc +# DEPENDS_ON ROSS mpi +# ) +# find_package(Python COMPONENTS Interpreter) + +set(actor_srcs ${CMAKE_CURRENT_SOURCE_DIR}/Funnel.hh) + +add_custom_command(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc + COMMAND ${Python_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/generateInfrastructure.py ${actor_srcs} > ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc + MAIN_DEPENDENCY generateInfrastructure.py + DEPENDS generateInfrastructure.py ${actor_srcs} + VERBATIM) + +set_source_files_properties( + ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc + PROPERTIES GENERATED TRUE +) + +if (WCS_HAS_ROSS) + add_executable(hybrid-bin + hybrid.cpp + hybrid.hpp + actor.hh + Checkpoint.hh + simtime.hh + actor.cc + ${actor_srcs} + ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc + ) + target_link_libraries(hybrid-bin PUBLIC ROSS::ROSS MPI::MPI_CXX wcs ${LIB_FILESYSTEM}) + target_include_directories(hybrid-bin SYSTEM PUBLIC ${ROSS_INCLUDE_DIRS}) + target_include_directories(hybrid-bin PUBLIC + $ + $ + $ + $) + + set_target_properties(hybrid-bin PROPERTIES OUTPUT_NAME hybrid) +endif() + +set_full_path(THIS_DIR_HEADERS + ) + +set_full_path(THIS_DIR_SOURCES + ) + +# # Add the subdirectories + + +# set(THIS_DIR_SOURCES "${THIS_DIR_SOURCES}" "${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc") + +# # Propagate the files up the tree +set(WCS_HEADERS "${WCS_HEADERS}" "${THIS_DIR_HEADERS}" PARENT_SCOPE) +set(WCS_SOURCES "${WCS_SOURCES}" "${THIS_DIR_SOURCES}" PARENT_SCOPE) + + diff --git a/src/hybrid/Checkpoint.hh b/src/hybrid/Checkpoint.hh new file mode 100644 index 0000000..cee1f6f --- /dev/null +++ b/src/hybrid/Checkpoint.hh @@ -0,0 +1,46 @@ +#pragma once + +#include +#include +#include + +// This could be done any number of ways faster/better than a list, +// but keep it for now. We can optimize this later, so long as we +// have the following interface. +class CHECKPOINTQueue { + public: + void push_back(std::string newString) { queue.push_back(newString); } + void pop_back() { queue.pop_back(); } + std::string back() { return queue.back(); } + std::string front() { return queue.front(); } + void pop_front() { queue.pop_front(); } + + std::list queue; +}; + +class InputCHECKPOINT { + public: + std::string str() { return sss.str(); } + + template + InputCHECKPOINT& operator&(const TTT& ttt) { + sss.write(reinterpret_cast(&ttt), sizeof(TTT)); + return *this; + } + + std::stringstream sss; +}; + +class OutputCHECKPOINT { + public: + OutputCHECKPOINT(std::string newString) { + sss.str(newString); + } + template + OutputCHECKPOINT& operator&(TTT& ttt) { + sss.read(reinterpret_cast(&ttt), sizeof(TTT)); + return *this; + } + + std::stringstream sss; +}; diff --git a/src/hybrid/Funnel.hh b/src/hybrid/Funnel.hh new file mode 100644 index 0000000..bde220a --- /dev/null +++ b/src/hybrid/Funnel.hh @@ -0,0 +1,609 @@ +#pragma once + +#include "actor.hh" +#include "reaction_network/network.hpp" +// #include "utils/graph_factory.hpp" +// #include +#include +#include +#include +#include "wcs_types.hpp" + +#define MAX_TIME 100 //FIXME! + +typedef wcs::Network::v_desc_t SpeciesTag; +typedef std::size_t ReactionTag; +typedef wcs::reaction_rate_t Real; +typedef std::size_t SpeciesValue; + + +typedef std::function&)> RateFunction; + + +struct NoopMsg { +}; +struct RateMsg { + Real rate; +}; +struct ReactionMsg { + ReactionTag reaction; +}; +struct FunnelRateExchangeMsg { + ReactionTag reaction; + bool nrmMode; + Real rate; + Real normalContrib; +}; +struct SwitchModeMsg { + ReactionTag reaction; + bool nrmMode; +}; + +class NextReactionMethodCrucible : public LP { + public: + + ACTOR_DEF(); + + SLOT(firedMyself, NoopMsg); + void FORWARD(firedMyself)(Time tnow, NoopMsg& msg, tw_bf*) { + ReactionMsg r; + r.reaction = _rxntag; + fire(0, r); + _countdown = newCountdown(); + fireNextReaction(tnow); + } + void COMMIT(firedMyself)(Time tnow, NoopMsg& msg, tw_bf*) { + //std::cout << tnow << ": Fired Reaction " << _rxntag << "!\n"; + } + template + void CHECKPOINT(firedMyself)(Checkpointer& ar, Time tnow, NoopMsg& msg, tw_bf*) { + ar & _countdown & _messageOutstanding & _lastUpdatedTime & _simpleRand; + } + + SLOT(updatedRate, RateMsg); + void FORWARD(updatedRate)(Time tnow, RateMsg& msg, tw_bf*) { + if (_rate != 0) { + EVENT_CANCEL(_messageOutstanding); + _countdown -= countdownChange(tnow); + } + _rate = msg.rate; + if (_rate != 0) { + fireNextReaction(tnow); + } + } + template + void CHECKPOINT(updatedRate)(Checkpointer& ar, Time tnow, RateMsg& msg, tw_bf*) { + ar & _messageOutstanding & _countdown & _rate & _lastUpdatedTime & _simpleRand; + } + + NextReactionMethodCrucible() { + _rate = 0; + _countdown = newCountdown(); + _messageOutstanding = PDES_NULL_EVENT; + } + + SIGNAL(fire,ReactionMsg); + + void setTag(const ReactionTag mytag) { + _rxntag = mytag; + } + + void setSeed(const int seed) { + _simpleRand.seed(seed); + } + + private: + ReactionTag _rxntag; + Real _countdown; + Real _rate; + std:: mt19937 _simpleRand; + PdesEvent _messageOutstanding; + Real _lastUpdatedTime; + // PdesLpid _self_lpid; + + void fireNextReaction(Time tnow) { + _lastUpdatedTime = tnow; + Time timeToFire = firingTime(_countdown,_rate); + NoopMsg msg; + _messageOutstanding = EVENT_SEND(twlp->gid, firedMyself, timeToFire, msg); + } + + Real newCountdown() { + //From Anderson, Algo 3 + //choose a uniform random number 0-1 + // std::random_device rd; //Will be used to obtain a seed for the random number engine + // std::mt19937 gen(rd()); //Standard mersenne_twister_engine seeded with rd() + std::uniform_real_distribution<> dis(0.0, 1.0); + Real uniform_random_number; + uniform_random_number = dis(_simpleRand); + + // std::cout << "Uniform: " << uniform_random_number << '\n'; + return log(1/uniform_random_number); + } + + Real countdownChange(Time tnow) { + return _rate*(tnow-_lastUpdatedTime); + } + + Real firingTime(Real countdown, Real rate) { + return countdown/rate; + } +}; + +extern +std::vector> stoic; + +class Funnel : public LP { + public: + + ACTOR_DEF(); + + Funnel() { + //Things needed to initiate the model + _rateNRM = 0; + _messageOutstanding = PDES_NULL_EVENT; + _numReactionsInBox = 0; + + + //tolerances the users can adjust + _tolerance = 0.003; + _sigma_max = 6; + _numReactionCutoff = 100; + _maxTimestep = MAX_TIME; + } + + + SLOT(firedReaction, ReactionMsg); + void FORWARD(firedReaction)(Time tnow, ReactionMsg& msg, tw_bf* twbf) { + //update species + for (const auto& kv : stoic[msg.reaction]) { + auto iter = _indexFromSpecTag.find(kv.first); + if (iter != _indexFromSpecTag.end()) { + _species[iter->second] += kv.second; + } + } + updateSpeciesBasedOnRates(tnow); + bool allInBounds = checkRateBounds(tnow); + if (!allInBounds) { + interruptNextWakeup(tnow); + setupFutureEvents(tnow, true); + } + } + void BACKWARD(firedReaction)(Time tnow, ReactionMsg& msg, tw_bf*) { + } + void COMMIT(firedReaction)(Time tnow, ReactionMsg& msg, tw_bf*) { + //std::cout << tnow << ": NRM " << msg.reaction << "->" << _rxntag << std::endl; + //printSpecies(tnow); + } + template + void CHECKPOINT(firedReaction)(Checkpointer& ar, Time tnow, ReactionMsg& msg, tw_bf*) { + } + + + SLOT(changedRate, FunnelRateExchangeMsg); + void FORWARD(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf* twbf) { + updateSpeciesBasedOnRates(tnow); + updateSavedRates(msg.reaction, msg.nrmMode, msg.rate, msg.normalContrib); + interruptNextWakeup(tnow); + //std::cout << tnow << ": %%% " << ((tnow-_lastRecomputeTime)>=_maxTimestep) << " " << checkRateBounds(tnow) << std::endl; + setupFutureEvents(tnow, (tnow-_lastRecomputeTime)>=_maxTimestep || !checkRateBounds(tnow)); + } + void BACKWARD(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { + } + void COMMIT(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { + // std::cout << tnow << ": Changed rate for " << msg.reaction << "->" + // << _rxntag << ": " + // << (msg.nrmMode ? "NRM " : "SDE ") + // << msg.rate << " " + // << msg.normalContrib << std::endl; + printSpecies(tnow); + } + template + void CHECKPOINT(changedRate)(Checkpointer& ar, Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { + } + + SLOT(wakeup, NoopMsg); + void FORWARD(wakeup)(Time tnow, NoopMsg& msg, tw_bf* twbf) { + // This gets called when we're sure that we're going to go + // out-of-bounds due to our tau condition, based on what we + // currently know about the SDE rates. + updateSpeciesBasedOnRates(tnow); + setupFutureEvents(tnow, (tnow-_lastRecomputeTime)>=_maxTimestep || !checkRateBounds(tnow)); + } + void COMMIT(wakeup)(Time tnow, NoopMsg& msg, tw_bf*) { + //std::cout << tnow << ": Woke Reaction " << _rxntag << "!\n"; + //printSpecies(tnow); + } + template + void CHECKPOINT(wakeup)(Checkpointer& ar, Time tnow, NoopMsg& msg, tw_bf*) { + } + + SLOT(changedMode, SwitchModeMsg); + void FORWARD(changedMode)(Time tnow, SwitchModeMsg& msg, tw_bf*) { + updateSpeciesBasedOnRates(tnow); + int irxn = _indexFromRxnTag[msg.reaction]; + updateSavedRates(msg.reaction, msg.nrmMode, _rxnRate[irxn], _normalContrib[irxn]); + bool allInBounds = checkRateBounds(tnow); + if (!allInBounds) { + interruptNextWakeup(tnow); + setupFutureEvents(tnow, true); + } + } + void COMMIT(changedMode)(Time tnow, SwitchModeMsg& msg, tw_bf*) { + //std::cout << tnow << ": Changed mode " << msg.reaction << "->" << _rxntag << ": " << msg.nrmMode << std::endl; + } + template + void CHECKPOINT(changedMode)(Checkpointer& ar, Time tnow, SwitchModeMsg& msg, tw_bf*) { + } + + SLOT(printData, NoopMsg); + void FORWARD(printData)(Time tnow, NoopMsg& msg, tw_bf*) { + } + void COMMIT(printData)(Time tnow, NoopMsg& msg, tw_bf*) { + printSpecies(tnow); + } + + SIGNAL(updateRateNRM, RateMsg); + SIGNAL(updateRateFunnel, FunnelRateExchangeMsg); + SIGNAL(updateModeFunnel, SwitchModeMsg); + + void addSpecies(SpeciesTag tag, SpeciesValue initialValue, SpeciesValue trustRegion) { + int oldSize = _indexFromSpecTag.size(); + _indexFromSpecTag[tag] = oldSize; + _tagFromSpecIndex.push_back(tag); + _species.push_back(initialValue); + _lowerBound.push_back(0); + _upperBound.push_back(0); + _speciesRate.push_back(0); + _trustRegion.push_back(trustRegion); + } + + void setTag(ReactionTag tag) { + _rxntag = tag; + } + + void setSeed(const int seed) { + _simpleRand.seed(seed); + } + + void setupSendToReaction(ReactionTag otherRxnTag, tw_lpid destID) { + if (otherRxnTag != _rxntag) { + CONNECT(*this,Funnel::updateRateFunnel,destID,Funnel::changedRate); + CONNECT(*this,Funnel::updateModeFunnel,destID,Funnel::changedMode); + } + } + + void setupRecvFromReaction(ReactionTag otherRxnTag) { + int index = _tagFromRxnIndex.size(); + _tagFromRxnIndex.push_back(otherRxnTag); + _indexFromRxnTag[otherRxnTag] = index; + _sqrtTimestep.push_back(0); + _normalContrib.push_back(0); + _rxnRate.push_back(0); + _rxnOverflow.push_back(0); + _rxnNrmMode.push_back(0); + if (otherRxnTag == _rxntag) { + _thisrxnindex = index; + } + } + + bool dependsOnSpecies(int specTag) { + return _indexFromSpecTag.find(specTag) != _indexFromSpecTag.end(); + } + + void beginReactions(Time tnow) { + _lastRecomputeTime = tnow; + _lastUpdatedTime = tnow; + + //recalculate rate + setupFutureEvents(tnow, false); + printSpecies(tnow); + } + + FunnelRateExchangeMsg initialRate(Time tnow, bool recompute) { + + FunnelRateExchangeMsg msg; + msg.reaction = _rxntag; + + if (recompute) { + msg.rate = recalculateRate(tnow); + } else { + msg.rate = _rxnRate[_thisrxnindex]; + } + Real delta_time; + getNextMode(msg.rate, delta_time, msg.nrmMode); + if (recompute) { + std::normal_distribution<> dis(0,1); + Real unit_normal_random = dis(_simpleRand); + msg.normalContrib = unit_normal_random * std::sqrt(msg.rate); + } + return msg; + } + + void setInitialRate(FunnelRateExchangeMsg& msg) { + updateSavedRates(msg.reaction, msg.nrmMode, msg.rate, msg.normalContrib); + } + + void printSpecies(Time tnow) { + std::cout << tnow << ": " + << "Reaction " << _rxntag << ": (wakeup " << _wakeupTime << ") (nrm=" << _rxnNrmMode[_thisrxnindex] << " rate=" << _rxnRate[_thisrxnindex] << ")" ; + for (std::size_t ii=0; ii<_tagFromSpecIndex.size(); ii++) { + std::cout << " {" << _tagFromSpecIndex[ii] << ": " + << "(" << _lowerBound[ii] << "<" + << _species[ii] + << "<" << _upperBound[ii] << ") " + << _speciesRate[ii] + << "}"; + } + std::cout << std::endl; + // int both_species = 0; + // for (int ii=0; ii<_tagFromSpecIndex.size(); ii++) { + // if (_tagFromSpecIndex[ii] == 2 || _tagFromSpecIndex[ii] == 4){ + // both_species = both_species +1; + // } + // } + // if (both_species == 2 ) { + // std::cout << tnow << " "; + // for (int ii=0; ii<_tagFromSpecIndex.size(); ii++) { + // if (_tagFromSpecIndex[ii] == 2 || _tagFromSpecIndex[ii] == 4){ + // std::cout << _tagFromSpecIndex[ii] << ": " << _species[ii] << " "; + // } + // } + // std::cout << std::endl; + // } + + } + + public: + + RateFunction _rateFunction; //function to use for rate calculations + Real _tolerance; // tolerance for the tau independence condition + Real _sigma_max; // tolerance for switching between SDE and NRM + int _numReactionCutoff; // min number of reactions before we consider SDE + Real _maxTimestep; // maximum timestep allowed in SDE mode. + + + private: + std::unordered_map _indexFromSpecTag; // internal index from tag + std::vector _tagFromSpecIndex; // tag from internal index + std::vector _species; //species counts + std::vector _speciesRate; //first derivative of the species + std::vector _lowerBound; // lb of the tau independence region + std::vector _upperBound; // ub of the tau independence region + std::vector _trustRegion; //how far should we trust this species for a trust region? Think of this as a max for the ub/lb, independent of the derivative + + std::unordered_map _indexFromRxnTag; //internal index from tag + std::vector _tagFromRxnIndex; // tag from internal index + std::vector _sqrtTimestep; // sqrt of the timestep, used for SDE + std::vector _normalContrib; // sqrt(rate)*normal, used for SDE + std::vector _rxnOverflow; // Accumulator for SDE updates + std::vector _rxnRate; // rate for each reaction + std::vector _rxnNrmMode; // boolean, is the rxn in nrmMode + + Real _lastUpdatedTime; + Real _lastRecomputeTime; + Real _wakeupTime; + + ReactionTag _rxntag; + int _thisrxnindex; + + Real _rateNRM; + PdesEvent _messageOutstanding; + Real _numReactionsInBox; + + std:: mt19937 _simpleRand; + + Real recalculateRate(Time tnow) { + //recalculate rate + Real newRate = _rateFunction(_species); + _numReactionsInBox = _numReactionCutoff+1; + //compute gradient + std::vector gradient(_species.size()); + for (unsigned int ii=0; ii<_species.size(); ii++) { + std::vector modSpecies(_species); + modSpecies[ii]++; + Real gradRate = _rateFunction(modSpecies); + gradient[ii] = gradRate-newRate; + } + + int nonZeroGradCount = 0; + for (auto grad : gradient) { + if (grad!=0) { + nonZeroGradCount++; + } + } + //compute bounding box + for (unsigned int ii=0; ii<_species.size(); ii++) { + Real maxDelta = _trustRegion[ii]; + Real delta; + if (gradient[ii] == 0) { + delta = maxDelta; + } else { + delta = _tolerance*newRate/fabs(gradient[ii])/sqrt(nonZeroGradCount)/2; + if (delta > maxDelta) { + delta = maxDelta; + } + } + _lowerBound[ii] = _species[ii]-delta; + _upperBound[ii] = _species[ii]+delta; + + SpeciesTag specTag = _tagFromSpecIndex[ii]; + auto iter = stoic[_rxntag].find(specTag); + if (iter != stoic[_rxntag].end()) { + int sss = ((iter->second >= 0) ? iter->second : -iter->second); + if (sss) { + _numReactionsInBox = std::min(_numReactionsInBox, delta/sss); + } + } + } + _lastRecomputeTime = tnow; + return newRate; + } + + bool checkRateBounds(Time tnow) const { + //are we within our bounds + bool allWithinBounds=true; + for (unsigned int ii=0; ii<_species.size(); ii++) { + if (!(_lowerBound[ii] <= _species[ii] && _species[ii] <= _upperBound[ii])) { + allWithinBounds = false; + break; + } + } + return allWithinBounds; + } + + void updateSpeciesBasedOnRates(Time tnow) { + + // Use SDE euler method to update the rates for things in SDE + // mode. + if (tnow <= _lastUpdatedTime) { + return; + } + + for (unsigned int irxn=0; irxn<_tagFromRxnIndex.size(); irxn++) { + if (_rxnNrmMode[irxn]==false && _rxnRate[irxn] != 0) { + _rxnOverflow[irxn] += _rxnRate[irxn]*(tnow-_lastUpdatedTime); + if (_normalContrib[irxn] != 0) { + Real c_old = _sqrtTimestep[irxn]; + Real c_now = std::sqrt(tnow - _lastUpdatedTime + c_old * c_old); + _rxnOverflow[irxn] += _normalContrib[irxn]*(c_now-c_old); + _sqrtTimestep[irxn] = c_now; + } + if (_rxnOverflow[irxn] >= 1) { + int rxnCount = int(_rxnOverflow[irxn]); + for (const auto& kv : stoic[_tagFromRxnIndex[irxn]]) { + auto iter = _indexFromSpecTag.find(kv.first); + if (iter != _indexFromSpecTag.end()) { + _species[iter->second] += rxnCount*kv.second; + } + } + _rxnOverflow[irxn] -= rxnCount; + } + } + } + _lastUpdatedTime = tnow; + } + + void updateSavedRates(ReactionTag otherRxnTag, bool nrmMode, Real rate, Real normalContrib) { + int irxn = _indexFromRxnTag[otherRxnTag]; + _normalContrib[irxn] = normalContrib; + _sqrtTimestep[irxn] = 0; + Real delta_rate = rate-_rxnRate[irxn]; + _rxnRate[irxn] = rate; + _rxnNrmMode[irxn] = nrmMode; + if (delta_rate!=0) { + for (const auto& kv : stoic[otherRxnTag]) { + auto iter = _indexFromSpecTag.find(kv.first); + if (iter != _indexFromSpecTag.end()) { + _speciesRate[iter->second] += kv.second*delta_rate; + } + } + } + } + + void interruptNextWakeup(Time tnow) { + EVENT_CANCEL(_messageOutstanding); + } + + void scheduleNextWakeup(Time tnow, Real delta_time) { + NoopMsg msg; + _messageOutstanding = EVENT_SEND(twlp->gid, wakeup, delta_time, msg); + _wakeupTime = tnow+delta_time; + } + + void setupFutureEvents(Time tnow, bool recompute) { + + bool nrmMode; + Real delta_time; + Real nextRate; + if (recompute) { + nextRate = recalculateRate(tnow); + } else { + nextRate = _rxnRate[_thisrxnindex]; + } + getNextMode(nextRate, delta_time, nrmMode); + if (recompute) { + FunnelRateExchangeMsg newRateMsg; + newRateMsg.reaction = _rxntag; + newRateMsg.rate = nextRate; + newRateMsg.nrmMode = nrmMode; + + std::normal_distribution<> dis(0,1); + Real unit_normal_random = dis(_simpleRand); + newRateMsg.normalContrib = unit_normal_random * std::sqrt(nextRate); + + updateSavedRates(newRateMsg.reaction, newRateMsg.nrmMode, newRateMsg.rate, newRateMsg.normalContrib); + updateRateFunnel(0, newRateMsg); + } + Real newRateNRM = (nrmMode==true ? nextRate : 0); + if (newRateNRM != _rateNRM) { + RateMsg newRateMsg; + newRateMsg.rate = newRateNRM; + updateRateNRM(0, newRateMsg); + _rateNRM = newRateNRM; + + if (!recompute) { + SwitchModeMsg newSwitchMsg; + newSwitchMsg.reaction = _rxntag; + newSwitchMsg.nrmMode = nrmMode; + updateSavedRates(_rxntag, nrmMode, nextRate, _normalContrib[_thisrxnindex]); + updateModeFunnel(0, newSwitchMsg); + } + } + scheduleNextWakeup(tnow, delta_time); + } + + void getNextMode(Real nextRate, Real& delta_time, bool& nextNrmMode) const { + + delta_time = findWakeupInterval(nextRate); + if (_numReactionsInBox > _numReactionCutoff && nextRate*(delta_time+_sqrtTimestep[_thisrxnindex]*_sqrtTimestep[_thisrxnindex]) > _sigma_max*_sigma_max) { + nextNrmMode = false; + } else { + nextNrmMode = true; + } + } + + Real findWakeupInterval(Real nextRate) const { + // find when species will cross boundary + Real minTime = _maxTimestep; + Real delta_rate = nextRate-(_rxnRate[_thisrxnindex]); + + for (unsigned int ispec=0; ispec<_species.size(); ispec++) { + Real thisRate = _speciesRate[ispec]; + if (delta_rate != 0) { + auto iter=stoic[_rxntag].find(_tagFromSpecIndex[ispec]); + if (iter != stoic[_rxntag].end()) { + thisRate += delta_rate*iter->second; + } + } + if (thisRate > 0) { + minTime = std::min(minTime,(_upperBound[ispec]+1 - _species[ispec]) / thisRate); + } else if (thisRate < 0) { + minTime = std::min(minTime,(_lowerBound[ispec]-1 - _species[ispec]) / thisRate); + } else { + //do nothing + } + //std::cout << "%%% " << thisRate << " " << minTime << std::endl; + } + return minTime; + } +}; + + +class Heartbeat : public LP { + public: + + ACTOR_DEF(); + + void setInterval(double interval) { _interval = interval; } + + SLOT(activated, NoopMsg); + void FORWARD(activated)(Time tnow, NoopMsg& msg, tw_bf*) { + activate(_interval, msg); + } + + SIGNAL(activate, NoopMsg); + + private: + double _interval; +}; diff --git a/src/hybrid/actor.cc b/src/hybrid/actor.cc new file mode 100644 index 0000000..1d97998 --- /dev/null +++ b/src/hybrid/actor.cc @@ -0,0 +1,13 @@ + +#include "actor.hh" +#include + +void my_event_cancel(PdesEvent evtptr) { + if (evtptr) { + tw_event_rescind(evtptr); + } +} + +tw_lpid BaseLP::getID() const { + return twlp->gid; +} diff --git a/src/hybrid/actor.hh b/src/hybrid/actor.hh new file mode 100644 index 0000000..082976b --- /dev/null +++ b/src/hybrid/actor.hh @@ -0,0 +1,135 @@ +#pragma once + +#include +#include +#include +#include +// #include "wcs-ross-bf.hpp" +#include "simtime.hh" +#include "Checkpoint.hh" + +//// Classes that make the whole API work + +struct EventData; //will be defined by a perl script. +typedef std::size_t Tag; +typedef tw_lpid PdesLpid; +typedef tw_event* PdesEvent; +#define PDES_NULL_EVENT NULL + +class BaseLP { + public: + virtual ~BaseLP() {} + + struct ObserverData { + tw_lpid dest; + Tag slot; + }; + tw_lp* twlp; //this has to be initialized properly! TODO + + tw_lpid getID() const; +}; + +extern CHECKPOINTQueue checkpointQueue; + +template +class LP : public BaseLP { + public: + typedef void (Dependent::*MethodPtr)(Time tnow, void*, tw_bf*); + static void init_lp(void *state, tw_lp *me); + static void FORWARD_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); + static void BACKWARD_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); + static void COMMIT_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); + static void final_lp(void *state, tw_lp *me); + static MethodPtr FORWARDDispatch[]; + static MethodPtr BACKWARDDispatch[]; + static MethodPtr COMMITDispatch[]; + + template + void FORWARDWithCHECKPOINT(Time tnow, typename SlotType::MsgType& msg, tw_bf* twbf) { + InputCHECKPOINT ic; + static_cast(this)->Dependent::template CHECKPOINTFull(ic, tnow, msg, twbf); + static_cast(this)->Dependent::template FORWARDFull(tnow, msg, twbf); + checkpointQueue.push_back(ic.str()); + } + template + void BACKWARDWithCHECKPOINT(Time tnow, typename SlotType::MsgType& msg, tw_bf* twbf) { + static_cast(this)->Dependent::template BACKWARDFull(tnow, msg, twbf); + OutputCHECKPOINT oc(checkpointQueue.back()); + static_cast(this)->Dependent::template CHECKPOINTFull(oc, tnow, msg, twbf); + checkpointQueue.pop_back(); + } + template + void COMMITWithCHECKPOINT(Time tnow, typename SlotType::MsgType& msg, tw_bf* twbf) { + static_cast(this)->Dependent::template COMMITFull(tnow, msg, twbf); + checkpointQueue.pop_front(); + } +}; + +#define ACTOR_DEF() \ + template \ + void FORWARDFull(Time tnow, typename SlotType::MsgType& msg, tw_bf*); \ + template \ + void BACKWARDFull(Time tnow, typename SlotType::MsgType& msg, tw_bf*) {} \ + template \ + void COMMITFull(Time tnow, typename SlotType::MsgType& msg, tw_bf*) {} \ + template \ + inline void CHECKPOINTFull(Archiver& ar, Time tnow, typename SlotType::MsgType& msg, tw_bf*); + + +#define TAGIFY(SlotType) SlotType::tag + +#define FORWARD(SlotType) FORWARD_##SlotType +#define BACKWARD(SlotType) BACKWARD_##SlotType +#define COMMIT(SlotType) COMMIT_##SlotType +#define CHECKPOINT(SlotType) CHECKPOINT_##SlotType + +#define SLOT(name, msgType) \ +class name { \ + public: \ + typedef msgType MsgType; \ + static Tag tag; \ +} + +#define SIGNAL(name, type) \ + class name##_signal_class { public: typedef type MsgType; }; \ + void name(Time trecv, type& msg) { \ + /*for all LPs attached to this signal*/ \ + for (auto&& observerData : name##_listeners) { \ + my_event_send(observerData.dest, trecv, twlp, observerData.slot, msg); \ + } \ + } \ + std::list name##_listeners + +#define CONNECT(sendingLP, signal, dest, slot) signal_add_listener((sendingLP).signal##_listeners,dest) + +template +inline +typename std::enable_if::value>::type +signal_add_listener(std::list& listeners, const BaseLP& dest) { + tw_lpid destID = dest.getID(); + signal_add_listener(listeners, destID); +} + +template +inline +typename std::enable_if::value>::type +signal_add_listener(std::list& listeners, const tw_lpid destID) { + listeners.push_back({destID,TAGIFY(SlotType)}); +} + +#define EVENT_SEND(lpid, slotClass, msgTime, msg) my_event_send(lpid, msgTime, twlp, msg) + +#define EVENT_CANCEL(event) my_event_cancel(event) + +template +PdesEvent my_event_send(tw_lpid dest, Time trecv, tw_lp* twlp, typename SlotType::MsgType& msg) +{ + return my_event_send(dest, trecv, twlp, TAGIFY(SlotType), msg); +} + +template +PdesEvent my_event_send(tw_lpid dest, Time trecv, tw_lp* twlp, Tag tag, TTT& msg); + +void my_event_cancel(PdesEvent evtptr); + +std::size_t messageSize(); diff --git a/src/hybrid/generateInfrastructure.py b/src/hybrid/generateInfrastructure.py new file mode 100755 index 0000000..e5cbc1c --- /dev/null +++ b/src/hybrid/generateInfrastructure.py @@ -0,0 +1,188 @@ +#!/usr/bin/env python + +import re +from collections import OrderedDict + +class Slot: + def __init__(self, name, msg): + self.name = name + self.msg = msg + self.modes = set() + def hasMode(self, mode): + return mode in self.modes + def addMode(self, mode): + self.modes.add(mode) + +class Actor: + def __init__(self, name): + self.name = name + self.slots = {} + def addSlot(self,slot): + self.slots[slot.name] = slot + def getSlot(self,slotName): + return self.slots[slotName] + def getSlots(self): + return self.slots.values() + +def main(): + import sys + + filenames = sys.argv[1:] + + actors = OrderedDict() + messages = set() + for filename in filenames: + for line in open(filename, "r"): + mmm = re.search(r'public LP<([A-Za-z0-9_]+)>',line) + if mmm: + actorName = mmm.group(1) + actors[actorName] = Actor(actorName) + mmm = re.search(r'SLOT\(\s*([A-Za-z0-9_]+)\s*,\s*([A-Za-z0-9_]+)\s*\)',line) + if mmm: + slotName = mmm.group(1) + msgName = mmm.group(2) + actors[actorName].addSlot(Slot(slotName,msgName)) + for eventType in ['FORWARD','BACKWARD','COMMIT','CHECKPOINT']: + mmm = re.search("void "+eventType+r'\(([A-Za-z0-9_]+)\)',line) + if mmm: + slotName = mmm.group(1) + actors[actorName].getSlot(slotName).addMode(eventType) + + print('#include "actor.hh"') + print('#include "Checkpoint.hh"') + + for filename in filenames: + print('#include "%s"' % filename) + + print(''' +struct EventData { + Tag tag; + union { + char msg;''') + allMessages = set() + for actor in actors.values(): + for slot in actor.getSlots(): + allMessages.add(slot.msg) + ii=0 + for msg in allMessages: + print(" %s msg%d;" % (msg, ii)) + ii += 1 + print(''' + }; +}; + +std::size_t messageSize() { + return sizeof(EventData); +} + +CHECKPOINTQueue checkpointQueue; + +template +PdesEvent my_event_send(tw_lpid dest, Time trecv, tw_lp* twlp, Tag tag, TTT& msg) +{ + tw_stime recv_ts = TW_STIME_ADD(tw_now(twlp), trecv); + if (recv_ts >= g_tw_ts_end) { + return NULL; + } + tw_event* evtptr = tw_event_new(dest, trecv, twlp); + EventData* data = static_cast(tw_event_data(evtptr)); + data->tag = tag; + memcpy(&data->msg, &msg, sizeof(msg)); + tw_event_send(evtptr); + return evtptr; +}''') + for msg in allMessages: + print("template PdesEvent my_event_send(tw_lpid dest, Time trecv, tw_lp* twlp, Tag tag, %s& msg);" % msg) + for actor in actors.values(): + for slot in actor.getSlots(): + ttt = { + "actor" : actor.name, + "slot" : slot.name, + "msg" : slot.msg, + } + for eventType in "FORWARD","BACKWARD","COMMIT": + ttt['event'] = eventType + if slot.hasMode(eventType): + print("""template <> +inline void %(actor)s::%(event)sFull<%(actor)s::%(slot)s>(Time tnow, %(msg)s& msg, tw_bf* twbf) { %(event)s(%(slot)s)(tnow,msg,twbf); }""" % ttt) + if slot.hasMode('CHECKPOINT'): + print("""template <> +inline void %(actor)s::CHECKPOINTFull<%(actor)s::%(slot)s>(InputCHECKPOINT& ar, Time tnow, %(msg)s& msg, tw_bf* twbf) { CHECKPOINT(%(slot)s)(ar,tnow,msg,twbf); } +template <> +inline void %(actor)s::CHECKPOINTFull<%(actor)s::%(slot)s>(OutputCHECKPOINT& ar, Time tnow, %(msg)s& msg, tw_bf* twbf) { CHECKPOINT(%(slot)s)(ar,tnow,msg,twbf); } +""" % ttt) + for actor in actors.values(): + for (ii,slot) in enumerate(actor.getSlots()): + print("Tag %s::%s::tag = %d;" % (actor.name,slot.name,ii)) + + for actor in actors.values(): + for eventType in ["FORWARD","BACKWARD","COMMIT"]: + print("template<> %s::MethodPtr LP<%s>::%sDispatch[] = {" % (actor.name, actor.name,eventType)) + for slot in actor.getSlots(): + ttt = { + "actor" : actor.name, + "slot" : slot.name, + "event" : eventType, + } + if slot.hasMode('CHECKPOINT'): + method = '%(actor)s::%(event)sWithCHECKPOINT<%(actor)s::%(slot)s>' % ttt + else: + method = '%(actor)s::%(event)sFull<%(actor)s::%(slot)s>' % ttt + ttt['method'] = method + print(" reinterpret_cast<%(actor)s::MethodPtr>(&%(method)s)," % ttt) + print("0};") + + print(""" +///This stuff needs EventData to be defined, but doesn't need specific generation + +template +void LP::init_lp(void *state, tw_lp *thislp) { + Dependent* derivedObj = new (reinterpret_cast(state)) Dependent(); + BaseLP *lpptr = static_cast(derivedObj); + lpptr->twlp = thislp; +} + +template +void LP::FORWARD_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp) { + Time tnow = tw_now(thislp); + EventData *e = reinterpret_cast(evt); + Dependent* derivedObj = reinterpret_cast(state); + (derivedObj->*FORWARDDispatch[e->tag])(tnow, &e->msg, bf); +} + +template +void LP::BACKWARD_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp) { + Time tnow = tw_now(thislp); + EventData *e = reinterpret_cast(evt); + Dependent* derivedObj = reinterpret_cast(state); + (derivedObj->*BACKWARDDispatch[e->tag])(tnow, &e->msg, bf); +} + +template +void LP::COMMIT_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp) { + Time tnow = tw_now(thislp); + EventData *e = reinterpret_cast(evt); + Dependent* derivedObj = reinterpret_cast(state); + (derivedObj->*COMMITDispatch[e->tag])(tnow, &e->msg, bf); +} + +template +void LP::final_lp(void *state, tw_lp *) { + Dependent* derivedObj = reinterpret_cast(state); + derivedObj->~Dependent(); +} + +""") + for actor in actors.values(): + print(""" +template void LP<%(actor)s>::FORWARD_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); +template void LP<%(actor)s>::BACKWARD_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); +template void LP<%(actor)s>::COMMIT_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); +template void LP<%(actor)s>::final_lp(void *state, tw_lp *me); +template void LP<%(actor)s>::init_lp(void *state, tw_lp *me); + +""" % {"actor" : actor.name}) + + +if __name__=="__main__": + main() diff --git a/src/hybrid/hybrid.cpp b/src/hybrid/hybrid.cpp new file mode 100644 index 0000000..ae2ced1 --- /dev/null +++ b/src/hybrid/hybrid.cpp @@ -0,0 +1,455 @@ +/****************************************************************************** + * * + * Copyright 2020 Lawrence Livermore National Security, LLC and other * + * Whole Cell Simulator Project Developers. See the top-level COPYRIGHT * + * file for details. * + * * + * SPDX-License-Identifier: MIT * + * * + ******************************************************************************/ + +#if defined(WCS_HAS_CONFIG) +#include "wcs_config.hpp" +#else +#error "no config" +#endif + +#if defined(WCS_HAS_ROSS) + +#include +#include // memset +#include +#include +#include +#include "utils/write_graphviz.hpp" +#include "utils/timer.hpp" +#include "utils/to_string.hpp" +#include "reaction_network/network.hpp" +#include "reaction_network/species.hpp" +#include "reaction_network/reaction.hpp" +#include "reaction_network/vertex.hpp" +#include "hybrid.hpp" +#include "hybrid/simtime.hh" +#include "hybrid/Funnel.hh" + +#define OPTIONS "hi:s:r:t:" +static const struct option longopts[] = { + {"help", no_argument, 0, 'h'}, + {"interval", required_argument, 0, 'i'}, + {"seedbase", required_argument, 0, 's'}, + {"trust_region", required_argument, 0, 'r'}, + {"simulation_time", required_argument, 0, 't'}, + { 0, 0, 0, 0}, +}; + +using r_prop_t = wcs::Network::r_prop_t; +// /// The type of BGL vertex descriptor for graph_t +using v_desc_t = boost::graph_traits::vertex_descriptor; +double interval = 0.2; +int seedBase = 137; +int trust_region = 20000; +int simulation_time = 100; +int numReactions = 0; +int numSpecies = 0; +const wcs::Network::reaction_list_t* reaction_list; +const wcs::Network::species_list_t* species_list; +const LIBSBML_CPP_NAMESPACE::Model* model; +std::string filename = ""; +wcs::Network rnet; + +void do_nothing(void *,void *,void *) {} + +tw_peid ctr_map(tw_lpid gid) { + return (tw_peid) gid / g_tw_nlp; +} + +template +class RegisterBufferAndOffset { + public: + static std::vector buffer; + static int offset; + + static void init_f(void *state, tw_lp* me) + { + LP::init_lp(state, me); + buffer[me->gid-offset] = reinterpret_cast(state); + } + + static void init(const int bufferSize, const int _offset) + { + buffer.resize(bufferSize); + offset = _offset; + } + static TTT** getBuffer() { + return &buffer[0]; + } +}; + +tw_lptype Funnellps[] = { + { (init_f) RegisterBufferAndOffset::init_f, + (pre_run_f) do_nothing, + (event_f) Funnel::FORWARD_event, + (revent_f) Funnel::BACKWARD_event, + (commit_f) Funnel::COMMIT_event, + (final_f) Funnel::final_lp, + (map_f) ctr_map, + //NULL, + sizeof(Funnel) + }, + {0} +}; +template<> +std::vector RegisterBufferAndOffset::buffer = std::vector(); +template<> +int RegisterBufferAndOffset::offset = 0; + +tw_lptype NRMlps[] = { + { (init_f) RegisterBufferAndOffset::init_f, + (pre_run_f) do_nothing, + (event_f) NextReactionMethodCrucible::FORWARD_event, + (revent_f) NextReactionMethodCrucible::BACKWARD_event, + (commit_f) NextReactionMethodCrucible::COMMIT_event, + (final_f) NextReactionMethodCrucible::final_lp, + (map_f) ctr_map, + //NULL, + sizeof(NextReactionMethodCrucible) + }, + {0} +}; + +template<> +std::vector RegisterBufferAndOffset::buffer = std::vector(); +template<> +int RegisterBufferAndOffset::offset = 0; + + +tw_lptype Heartbeatlps[] = { + { (init_f) RegisterBufferAndOffset::init_f, + (pre_run_f) do_nothing, + (event_f) Heartbeat::FORWARD_event, + (revent_f) Heartbeat::BACKWARD_event, + (commit_f) Heartbeat::COMMIT_event, + (final_f) Heartbeat::final_lp, + (map_f) ctr_map, + //NULL, + sizeof(Heartbeat) + }, + {0} +}; + + +template<> +std::vector RegisterBufferAndOffset::buffer = std::vector(); +template<> +int RegisterBufferAndOffset::offset = 0; + +// void post_lpinit_setup(tw_pe* pe); + +tw_petype MyPEType = { + NULL, + post_lpinit_setup, + NULL, + NULL +}; + + +std::vector> stoic; + + +void print_usage(const std::string exec, int code) +{ + std::cerr << + "Usage: " << exec << " .xml\n" + " Run the simulation using the hybrid simulation mode which combines \n" + " Stochastic differential equations (SDEs) with stochastic simulation \n" + " algorithm (SSA) on a reaction network graph (.xml).\n" + "\n" + " OPTIONS:\n" + " -h, --help\n" + " Display this usage information\n" + "\n" + " -i, --interval\n" + " Specify the interval\n" + "\n" + " -s, --seedbase\n" + " Specify the seed for random number generator\n" + "\n" + " -r, --trust_region\n" + " Specify the trust region\n" + "\n" + " -t, --simulation_time\n" + " Specify the simulation time\n" + "\n"; + exit(code); +} + + +int main(int argc, char **argv) +{ + // get rid of error if compiled w/ MEMORY queues + //g_tw_memory_nqueues=1; + //g_tw_memory_nqueues = 16; // give at least 16 memory queue event + + //printf("Calling tw_opt_add\n"); + //tw_opt_add(app_opt); + + + int c; + int rc = EXIT_SUCCESS; + std::string outfile; + while ((c = getopt_long(argc, argv, OPTIONS, longopts, NULL)) != -1) { + switch (c) { + case 'h': /* --help */ + print_usage(argv[0], 0); + break; + case 'i': /* --interval */ + interval = atof(optarg); + break; + case 's': /* --seedbase */ + seedBase = static_cast(atoi(optarg)); + break; + case 'r': /* --trust_region */ + trust_region = static_cast(atoi(optarg)); + break; + case 't': /* --simulation_time */ + simulation_time = static_cast(atoi(optarg)); + break; + default: + print_usage(argv[0], 1); + break; + } + } + if (optind != (argc - 1)) { + print_usage (argv[0], 1); + } + std::string fn(argv[optind]); + filename = fn; + rnet.load(fn); + rnet.init(); + const wcs::Network::graph_t& g = rnet.graph(); + + LIBSBML_CPP_NAMESPACE::SBMLReader reader; + const LIBSBML_CPP_NAMESPACE::SBMLDocument* document + = reader.readSBML(fn); + + model = document->getModel(); + + + //FIXME, needs to come from graph => DONE + + reaction_list = &rnet.reaction_list(); + species_list = &rnet.species_list(); + numReactions = reaction_list->size(); + numSpecies = species_list->size(); + + //printf("Calling tw_init\n"); + tw_init(&argc, &argv); + //printf("tw_init returned.\n"); + // g_tw_nlp = numReactions+numReactions+numReactions; + g_tw_nlp = numReactions+numReactions+1; + g_tw_events_per_pe = g_tw_nlp * 100; + + //FIXME, yeom2 . I need to fill out these arrays with pointers to the actual objects. + RegisterBufferAndOffset::init(numReactions,0); + RegisterBufferAndOffset::init(numReactions,numReactions); + RegisterBufferAndOffset::init(1,2*numReactions); + + tw_define_lps(g_tw_nlp, messageSize()); + + for(int ii = 0; ii < numReactions ; ii++) { + tw_lp_settype(ii, &Funnellps[0]); + } + for(int ii = numReactions; ii < numReactions+numReactions; ii++) { + tw_lp_settype(ii, &NRMlps[0]); + } + tw_lp_settype(2*numReactions, &Heartbeatlps[0]); + + tw_pe_settype(&MyPEType); + g_st_sampling_end = simulation_time; + if (g_st_sampling_end > 0) { + g_tw_ts_end = g_st_sampling_end; // default 100000 + } else { + g_tw_ts_end = 1; // default 100000 + } + + // g_tw_events_per_pe =10; //default 800 + if( g_tw_mynode == 0 ) + { + std::cout << "=========================================" << std::endl; + std::cout << "WCS ROSS Configuration.............." << std::endl; + // std::cout << " run_id:\t" + std::string(run_id) << std::endl; + std::cout << " nlp_per_pe:\t" << g_tw_nlp << std::endl; + std::cout << " g_tw_ts_end:\t" << g_tw_ts_end << std::endl; + std::cout << " g_tw_events_per_pe:\t" << g_tw_events_per_pe << std::endl; + + std::cout << " gvt-interval:\t" << g_tw_gvt_interval << std::endl;; + std::cout << " extramem:\t" << g_tw_events_per_pe_extra << std::endl; + std::cout << " ......................................" << std::endl; + std::cout << " Num nodes:\t" << tw_nnodes() << std::endl; + // std::cout << " Message size:\t" << sizeof(WCS_Message) << std::endl; + std::cout << "=========================================" + << std::endl << std::endl; + } + + + tw_run(); + + tw_end(); + +// return 0; +return EXIT_SUCCESS; +} + +void post_lpinit_setup(tw_pe* pe) { + + Funnel** funnels = RegisterBufferAndOffset::getBuffer(); + NextReactionMethodCrucible** crucibles = RegisterBufferAndOffset::getBuffer(); + Heartbeat** heartbeats = RegisterBufferAndOffset::getBuffer(); + + + heartbeats[0]->setInterval(interval); + + for (int ireaction=0; ireactionsetTag(ireaction); + crucibles[ireaction]->setSeed(seedBase+ireaction); + + funnels[ireaction]->setTag(ireaction); + funnels[ireaction]->setSeed(seedBase+numReactions+ireaction); + } + + + const wcs::Network::graph_t& g = rnet.graph(); + + const wcs::Network::species_list_t& species_list = rnet.species_list(); + const wcs::Network::reaction_list_t& reaction_list = rnet.reaction_list(); + const LIBSBML_CPP_NAMESPACE::ListOfReactions* model_reactions + = model->getListOfReactions(); + const LIBSBML_CPP_NAMESPACE::ListOfSpecies* model_species + = model->getListOfSpecies(); + + + + for (size_t i = 0u; i < reaction_list.size(); ++i) { + // std::cout << "reaction" << i << std::endl; + const auto& vd = reaction_list[i]; + + const auto &reaction = *(model_reactions->get(g[vd].get_label())); + const LIBSBML_CPP_NAMESPACE::ListOfSpeciesReferences* reaction_reactants + = reaction.getListOfReactants(); + const LIBSBML_CPP_NAMESPACE::ListOfSpeciesReferences* reaction_modifiers + = reaction.getListOfModifiers(); + + // std::cout << "Reaction label " << g[vd].get_label() << std::endl; + const auto& rp = g[vd].property(); + funnels[i]->_rateFunction = rp.get_calc_rate_fn(); + // reactant species + // std::cout << "reactants " << std::endl; + for (const auto ei_in : boost::make_iterator_range(boost::in_edges(vd, g))) { + const auto vd_reactant = boost::source(ei_in, g); + const auto& sv_reactant = g[vd_reactant]; + + const auto& sp_reactant = sv_reactant.property(); + const auto stoichio = g[ei_in].get_stoichiometry_ratio(); + // Include only reactants without the constant reactants/modifiers which have been converted into modifiers + if (reaction_reactants->get(sv_reactant.get_label()) != NULL || reaction_modifiers->get(sv_reactant.get_label()) != NULL) { + funnels[i]->addSpecies(vd_reactant, sp_reactant.get_count(), trust_region); + } + } + } + + //Connect up the funnels to the crucibles they own. + for( int ii=0; ii> reactionConnectivity; + + for (size_t ireaction = 0u; ireaction < reaction_list.size(); ++ireaction) { + funnels[ireaction]->setupRecvFromReaction(ireaction); + std::unordered_map dependentSpecies; + std::unordered_map::iterator speciesit; + + const auto& vd = reaction_list[ireaction]; + + // reactant species + for (const auto ei_in : boost::make_iterator_range(boost::in_edges(vd, g))) { + const auto vd_reactant = boost::source(ei_in, g); + const auto& sv_reactant = g[vd_reactant]; + + const auto stoichio = g[ei_in].get_stoichiometry_ratio(); + // Include only reactants without modifiers + if (stoichio > 0) { + dependentSpecies.insert(std::make_pair(vd_reactant, stoichio * -1)); + } + } + // product species + for (const auto ei_out : boost::make_iterator_range(boost::out_edges(vd, g))) { + const auto vd_product = boost::target(ei_out, g); + const auto& sv_product = g[vd_product]; + + const auto stoichio = g[ei_out].get_stoichiometry_ratio(); + speciesit = dependentSpecies.find(vd_product); + if (speciesit == dependentSpecies.cend()) { + dependentSpecies.insert(std::make_pair(vd_product, stoichio)); + } else { + speciesit->second += stoichio; + } + } + stoic.push_back(dependentSpecies); + std::set rxnConnectSet; + for (size_t jreaction = 0u; jreaction < reaction_list.size(); ++jreaction) { + // const auto& vdj = reaction_list[jreaction]; + for (auto speciesTag : dependentSpecies) { + if (funnels[jreaction]->dependsOnSpecies(speciesTag.first)) { + CONNECT(*(crucibles[ireaction]),NextReactionMethodCrucible::fire, + *(funnels[jreaction]),Funnel::firedReaction); + rxnConnectSet.insert(jreaction); + if (ireaction != jreaction) { + funnels[ireaction]->setupSendToReaction(jreaction, funnels[jreaction]->getID()); + funnels[jreaction]->setupRecvFromReaction(ireaction); + } + break; + } + } + } + reactionConnectivity.emplace_back(std::move(rxnConnectSet)); + } + + //connect the printing to the heartbeat for output + CONNECT((*heartbeats[0]),Heartbeat::activate,(*heartbeats[0]),Heartbeat::activated); + for (int ireaction=0; ireactioninitialRate(tnow, true); + for (auto& jreaction : reactionConnectivity[ireaction]) { + funnels[jreaction]->setInitialRate(rate); + } + } + + //Setup the initial modes: + for (size_t ireaction = 0u; ireaction < reaction_list.size(); ++ireaction) { + auto rate = funnels[ireaction]->initialRate(tnow, false); + for (auto& jreaction : reactionConnectivity[ireaction]) { + funnels[jreaction]->setInitialRate(rate); + } + } + + //setup the initial messages + for (int ireaction=0; ireactionbeginReactions(tnow); + } + NoopMsg msg; + //heartbeats[0]->activate(tnow, msg); +} + + + + +#endif // defined(WCS_HAS_ROSS) diff --git a/src/hybrid/hybrid.hpp b/src/hybrid/hybrid.hpp new file mode 100644 index 0000000..ec12ac3 --- /dev/null +++ b/src/hybrid/hybrid.hpp @@ -0,0 +1,53 @@ +/****************************************************************************** + * * + * Copyright 2020 Lawrence Livermore National Security, LLC and other * + * Whole Cell Simulator Project Developers. See the top-level COPYRIGHT * + * file for details. * + * * + * SPDX-License-Identifier: MIT * + * * + ******************************************************************************/ + +#ifndef WCS_HYBRID_HPP +#define WCS_HYBRID_HPP + +#if defined(WCS_HAS_CONFIG) +#include "wcs_config.hpp" +#else +#error "no config" +#endif + +#if defined(WCS_HAS_ROSS) + +#include +#include +#include +#include "reaction_network/network.hpp" +#include "sim_methods/ssa_nrm.hpp" +#include "params/ssa_params.hpp" + + + +void do_nothing(void *,void *,void *); +tw_peid ctr_map(tw_lpid gid); +void post_lpinit_setup(tw_pe* pe); + + + + +// //----------------- +// // ROSS options +// //----------------- +// static unsigned int nlp_per_pe = 1u; +// static char run_id[1024] = "wcs-ross"; + +// const tw_optdef app_opt[] = +// { +// TWOPT_GROUP("WCS ROSS"), +// TWOPT_UINT("nlp", nlp_per_pe, "Number of LPs per processor"), +// TWOPT_CHAR("run", run_id, "User supplied run name"), +// TWOPT_END() +// }; + +#endif // defined(WCS_HAS_ROSS) +#endif // WCS_HYBRID_HPP diff --git a/src/hybrid/simtime.hh b/src/hybrid/simtime.hh new file mode 100644 index 0000000..6a5d415 --- /dev/null +++ b/src/hybrid/simtime.hh @@ -0,0 +1,10 @@ +#pragma once + +#ifndef KMCSIM_TYPES_ONLY +#include +#include +// #include "wcs-ross-bf.hpp" +#endif + +typedef tw_stime Time; + diff --git a/src/reaction_network/network.cpp b/src/reaction_network/network.cpp index 85a91a4..b6d504a 100644 --- a/src/reaction_network/network.cpp +++ b/src/reaction_network/network.cpp @@ -154,9 +154,9 @@ void Network::loadSBML(const std::string sbml_filename, const bool reuse) generate_cxx_code code_generator(lib_filename, !reuse); typename params_map_t::const_iterator pit; + struct timespec t_1, t_2; code_generator.generate_code(*model, m_dep_params_f, m_dep_params_nf, m_rate_rules_dep_map); - const std::string library_file = code_generator.compile_code(); using std::operator<<; @@ -192,6 +192,7 @@ void Network::init() m_species.reserve(num_vertices); v_iter_t vi, vi_end; + for (boost::tie(vi, vi_end) = boost::vertices(m_graph); vi != vi_end; ++vi) { const v_prop_t& v = m_graph[*vi]; const auto vt = static_cast(v.get_typeid()); @@ -238,6 +239,7 @@ void Network::init() const auto st = m_graph[ei_in].get_stoichiometry_ratio(); involved_species.insert(std::make_pair(m_graph[reactant].get_label(), std::make_pair(reactant, st))); + } } diff --git a/src/reaction_network/reaction_base.cpp b/src/reaction_network/reaction_base.cpp index 91b1ab6..57d2f0a 100644 --- a/src/reaction_network/reaction_base.cpp +++ b/src/reaction_network/reaction_base.cpp @@ -125,6 +125,14 @@ void ReactionBase::set_calc_rate_fn( m_calc_rate = calc_rate; } +const std::function< + reaction_rate_t (const std::vector&) + >& ReactionBase::get_calc_rate_fn() const +{ + return m_calc_rate; +} + + reaction_rate_t ReactionBase::calc_rate(std::vector&& params) { params.push_back(m_rate_const); diff --git a/src/reaction_network/reaction_base.hpp b/src/reaction_network/reaction_base.hpp index be799dd..37b475b 100644 --- a/src/reaction_network/reaction_base.hpp +++ b/src/reaction_network/reaction_base.hpp @@ -49,6 +49,10 @@ class ReactionBase : public VertexPropertyBase { void set_calc_rate_fn(const std::function< reaction_rate_t (const std::vector&) >& calc_rate); + const std::function< + reaction_rate_t (const std::vector&) + >& get_calc_rate_fn() const; + /// Overwrite the reaction rate to a given value void set_rate(const reaction_rate_t rate); reaction_rate_t get_rate() const; diff --git a/src/reaction_network/vertex.cpp b/src/reaction_network/vertex.cpp index 859c83f..df62d22 100644 --- a/src/reaction_network/vertex.cpp +++ b/src/reaction_network/vertex.cpp @@ -30,7 +30,7 @@ std::map Vertex::str_vt { }; Vertex::Vertex() -: m_type(_undefined_), m_label(""), +: m_type(_undefined_), m_label(""), m_annotation(""), m_pid(unassigned_partition), m_p(nullptr) { @@ -41,6 +41,7 @@ Vertex::Vertex(const Vertex& rhs) : m_type(rhs.m_type), m_typeid(rhs.m_typeid), m_label(rhs.m_label), + m_annotation(rhs.m_annotation), m_pid(rhs.m_pid), m_p((!rhs.m_p)? nullptr : rhs.m_p.get()->clone()) { @@ -53,6 +54,7 @@ Vertex::Vertex(Vertex&& rhs) noexcept { if (this != &rhs) { m_label = std::move(rhs.m_label); + m_annotation = std::move(rhs.m_annotation); m_p = std::move(rhs.m_p); reset(rhs); @@ -65,6 +67,7 @@ Vertex& Vertex::operator=(const Vertex& rhs) m_type = rhs.m_type; m_typeid = rhs.m_typeid; m_label = rhs.m_label; + m_annotation = rhs.m_annotation; m_pid = rhs.m_pid; m_p = (!rhs.m_p)? nullptr : rhs.m_p.get()->clone(); } @@ -77,6 +80,7 @@ Vertex& Vertex::operator=(Vertex&& rhs) noexcept m_type = rhs.m_type; m_typeid = rhs.m_typeid; m_label = std::move(rhs.m_label); + m_annotation = std::move(rhs.m_annotation); m_pid = rhs.m_pid; m_p = std::move(rhs.m_p); @@ -150,6 +154,15 @@ partition_id_t Vertex::get_partition() const { return m_pid; } +void Vertex::set_annotation(const std::string& f) +{ + m_annotation = f; +} + +std::string Vertex::get_annotation() const +{ + return m_annotation; +} /**@}*/ } // end of namespace wcs diff --git a/src/reaction_network/vertex.hpp b/src/reaction_network/vertex.hpp index d394bc4..fda8de4 100644 --- a/src/reaction_network/vertex.hpp +++ b/src/reaction_network/vertex.hpp @@ -70,7 +70,7 @@ class Vertex { #if defined(WCS_HAS_SBML) template - Vertex(const LIBSBML_CPP_NAMESPACE::Species& species, const G& g); + Vertex(const LIBSBML_CPP_NAMESPACE::Model& model, const LIBSBML_CPP_NAMESPACE::Species& species, const G& g); template Vertex(const LIBSBML_CPP_NAMESPACE::Model& model, const @@ -89,6 +89,8 @@ class Vertex { std::string get_label() const; void set_partition(const partition_id_t pid); partition_id_t get_partition() const; + void set_annotation(const std::string& f); + std::string get_annotation() const; template P& property() const; template P& checked_property() const; @@ -106,6 +108,7 @@ class Vertex { int m_typeid; ///< The vertex type in integer form used for GraphML parsing std::string m_label; ///< The vertex label partition_id_t m_pid; ///< The id of the partition to which this edge belongs + std::string m_annotation; ///< vertex annotation /** * The pointer to the detailed property object, which is polymorphic. @@ -153,7 +156,8 @@ Vertex::Vertex(const VertexFlat& flat, const G& g) #if defined(WCS_HAS_SBML) template -Vertex::Vertex(const LIBSBML_CPP_NAMESPACE::Species& species, const G& g) +Vertex::Vertex(const LIBSBML_CPP_NAMESPACE::Model& model, +const LIBSBML_CPP_NAMESPACE::Species& species, const G& g) : m_type(_species_), m_typeid(static_cast(_species_)), m_label(species.getIdAttribute()), @@ -162,15 +166,59 @@ Vertex::Vertex(const LIBSBML_CPP_NAMESPACE::Species& species, const G& g) { m_p = std::unique_ptr(new Species); - if (!isnan(species.getInitialAmount())) { + if (species.isSetInitialAmount()) { dynamic_cast(m_p.get())-> set_count(static_cast(species.getInitialAmount())); + // TODO: if units of species are mole then species Initial Amount has + // to be multiplied by Avogadro number + // TODO: species.getInitialConcentration() should be multiplied by // compartment volume and molarity (depending on the unit of // concentration) should be converted to count - } else if (!isnan(species.getInitialConcentration())) { + } else if (species.isSetInitialConcentration()) { + const LIBSBML_CPP_NAMESPACE::ListOfCompartments* compartment_list + = model.getListOfCompartments(); + const LIBSBML_CPP_NAMESPACE::ListOfUnitDefinitions* unit_definition_list + = model.getListOfUnitDefinitions(); + double compartment_size = compartment_list->get(species.getCompartment())->getSize(); + // Find compartment units + std::string compartment_unit ; + if (compartment_list->get(species.getCompartment())->isSetUnits()){ + compartment_unit = compartment_list->get(species.getCompartment())->getUnits(); + } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 3) { + compartment_unit = model.getVolumeUnits(); + } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 2) { + compartment_unit = model.getAreaUnits(); + } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 1) { + compartment_unit = model.getLengthUnits(); + } + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition_comp = unit_definition_list->get(compartment_unit); + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list_comp + = unit_definition_comp->getListOfUnits(); + unsigned int unitsSize_comp = unit_list_comp->size(); + double comp_unit = 1.0; + for (unsigned int iu = 0u; iu < unitsSize_comp; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list_comp->get(iu); + comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); + } + // Find substance units + std::string substance_unit ; + substance_unit = model.getSubstanceUnits(); + + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition_substance = unit_definition_list->get(substance_unit); + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list_substance + = unit_definition_substance->getListOfUnits(); + unsigned int unitsSize_substance = unit_list_substance->size(); + double sub_unit = 1.0; + for (unsigned int iu = 0u; iu < unitsSize_substance; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list_substance->get(iu); + sub_unit = sub_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); + } + + double avog_num = 6.02214e+23; dynamic_cast(m_p.get())-> - set_count(static_cast(species.getInitialConcentration())); + set_count(static_cast(species.getInitialConcentration() * (6.02214E23 * + compartment_size) * comp_unit * sub_unit)); } } @@ -209,12 +257,31 @@ Vertex::Vertex( //Add parameters const LIBSBML_CPP_NAMESPACE::ListOfParameters* parameter_list = model.getListOfParameters(); - unsigned int parametersSize = parameter_list->size(); - + unsigned int parametersSize = parameter_list->size(); //Add local parameters - const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list - = reaction.getKineticLaw()->getListOfLocalParameters(); - unsigned int num_local_parameters = local_parameter_list->size(); + if (model.getLevel() > 2) { + const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfLocalParameters(); + unsigned int num_local_parameters = reaction.getKineticLaw()->getNumLocalParameters(); + + // Create an unordered_set for all local parameters of reaction + for (unsigned int pi = 0u; pi < num_local_parameters; pi++) { + const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter + = local_parameter_list->get(pi); + lpset.insert(localparameter->getIdAttribute()); + } + } else { + const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfParameters(); + unsigned int num_local_parameters = reaction.getKineticLaw()->getNumParameters(); + // Create an unordered_set for all local parameters of reaction + for (unsigned int pi = 0u; pi < num_local_parameters; pi++) { + const LIBSBML_CPP_NAMESPACE::Parameter* localparameter + = local_parameter_list->get(pi); + lpset.insert(localparameter->getIdAttribute()); + } + } + sbml_utils sbml_o; pset=sbml_o.get_reaction_parameters(model, reaction); @@ -225,12 +292,6 @@ Vertex::Vertex( mpset.insert(parameter->getIdAttribute()); } - // Create an unordered_set for all local parameters of reaction - for (unsigned int pi = 0u; pi < num_local_parameters; pi++) { - const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter - = local_parameter_list->get(pi); - lpset.insert(localparameter->getIdAttribute()); - } // Put initial assignments in map using all_assignments_type @@ -284,16 +345,30 @@ Vertex::Vertex( } // reaction local parameters lpit = lpset.find(s_label) ; - if (lpit != lpset.end()) { - std::stringstream ss; - ss << local_parameter_list->get(s_label)->getValue(); - std::string parametervalue = ss.str(); - wholeformula = wholeformula + "var " + s_label - + " := " + parametervalue + "; "; - } + if (model.getLevel() > 2) { + const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfLocalParameters(); + if (lpit != lpset.end()) { + std::stringstream ss; + ss << local_parameter_list->get(s_label)->getValue(); + std::string parametervalue = ss.str(); + wholeformula = wholeformula + "var " + s_label + + " := " + parametervalue + "; "; + } + } else { + const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfParameters(); + if (lpit != lpset.end()) { + std::stringstream ss; + ss << local_parameter_list->get(s_label)->getValue(); + std::string parametervalue = ss.str(); + wholeformula = wholeformula + "var " + s_label + + " := " + parametervalue + "; "; + } + } } - //Add compartnents + //Add compartments const LIBSBML_CPP_NAMESPACE::ListOfCompartments* compartment_list = model.getListOfCompartments(); unsigned int compartmentsSize = compartment_list->size(); @@ -317,6 +392,12 @@ Vertex::Vertex( wholeformula = wholeformula + "m_rate := " + formula + ";"; dynamic_cast*>(m_p.get())->set_rate_formula(wholeformula); + if (reaction.isSetAnnotation()) { + std::string annot= reaction.getAnnotationString() ; + m_annotation = annot; + } else { + m_annotation = ""; + } #if !defined(WCS_HAS_EXPRTK) dynamic_cast*>(m_p.get())->ReactionBase::set_calc_rate_fn(reaction_function_rate); diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt index 47ccdb9..2faac5a 100644 --- a/src/utils/CMakeLists.txt +++ b/src/utils/CMakeLists.txt @@ -12,6 +12,7 @@ set_full_path(THIS_DIR_HEADERS rngen_impl.hpp samples_ssa.hpp sbml_utils.hpp + FunctionInterfaceForJIT.hpp seed.hpp state_io.hpp state_io_impl.hpp @@ -39,6 +40,7 @@ set_full_path(THIS_DIR_SOURCES print_vertices.cpp samples_ssa.cpp sbml_utils.cpp + FunctionInterfaceForJIT.cpp trajectory.cpp trace_ssa.cpp trace_generic.cpp diff --git a/src/utils/FunctionInterfaceForJIT.cpp b/src/utils/FunctionInterfaceForJIT.cpp new file mode 100644 index 0000000..b85f2be --- /dev/null +++ b/src/utils/FunctionInterfaceForJIT.cpp @@ -0,0 +1,343 @@ +/****************************************************************************** + * * + * Copyright 2020 Lawrence Livermore National Security, LLC and other * + * Whole Cell Simulator Project Developers. See the top-level COPYRIGHT * + * file for details. * + * * + * SPDX-License-Identifier: MIT * + * * + ******************************************************************************/ + +#include "utils/FunctionInterfaceForJIT.hpp" + +#if defined(WCS_HAS_CONFIG) +#include "wcs_config.hpp" +#else +#error "no config" +#endif + + +#if defined(WCS_HAS_SBML) +#include +#include +#endif // defined(WCS_HAS_SBML) + +#include +#include +// #include +// #include + +namespace wcs { +/** \addtogroup wcs_utils + * @{ */ + +FunctionInterfaceForJIT::FunctionInterfaceForJIT() +{} + +#if defined(WCS_HAS_SBML) +using ast_function_ptr = std::function; +using node_structure = std::tuple; +using node_map = std::unordered_map; +node_map nodetypes = { + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ABS,{nullptr,"std::abs(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCCOS,{nullptr,"std::acos(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCCOSH,{nullptr,"std::acosh(", "", ")"}}, + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCCOT,{nullptr,"std::(", "", ")"}}, // arccotangent + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCCOTH,{nullptr,"std::(", "", ")"}}, // Hyperbolic arccotangent + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCCSC,{nullptr,"std::(", "", ")"}}, // arccosecant + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCCSCH,{nullptr,"std::(", "", ")"}}, // Hyperbolic arccosecant + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCSEC,{nullptr,"std::(", "", ")"}}, // arcsecant + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCSECH,{nullptr,"std::(", "", ")"}}, // Hyperbolic arcsecant + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCSIN,{nullptr,"std::asin(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCSINH,{nullptr,"std::asinh(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCTAN,{nullptr,"std::atan(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCTANH,{nullptr,"std::atanh(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_CEILING,{nullptr,"std::ceil(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_COS,{nullptr,"std::cos(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_COSH,{nullptr,"std::cosh(", "", ")"}}, + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_COT,{nullptr,"std::(", "", ")"}}, //cotangent + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_COTH,{nullptr,"std::(", "", ")"}}, // Hyperbolic cotangent + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_CSC,{nullptr,"std::(", "", ")"}}, // cosecant + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_CSCH,{nullptr,"std::(", "", ")"}}, // Hyperbolic cosecant + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_DELAY,{nullptr,"std::(", "", ")"}}, //Delay + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_EXP,{nullptr,"std::exp(", "", ")"}}, + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_FACTORIAL,{nullptr,"std::(", "", ")"}}, // Factorial + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_FLOOR,{nullptr,"std::floor(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_LN,{nullptr,"std::log(", "", ")"}}, // natural log + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_LOG,{nullptr,"", "", ""}}, // function pointer------- + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_PIECEWISE,{nullptr,"std::(", "", ")"}}, // Piecewise + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_POWER,{nullptr,"std::pow(", ", ", ")"}}, // NOT SURE + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_SEC,{nullptr,"std::(", "", ")"}}, // Secant + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_SECH,{nullptr,"std::(", "", ")"}}, // Hyperbolic Secant + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_SIN,{nullptr,"std::sin(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_SINH,{nullptr,"std::sinh(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_TAN,{nullptr,"std::tan(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_TANH,{nullptr,"std::tanh(", "", ")"}}, + + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_MIN,{nullptr,"std::min({", ", ", "})"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_MAX,{nullptr,"std::max({", ", ", "})"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_QUOTIENT,{nullptr,"", "", ""}}, // function pointer div ------- + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_RATE_OF,{nullptr,"", "", ""}}, // function pointer rate_of ------- + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_REM,{nullptr,"std::remainder(", ", ", ")"}}, + + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ROOT,{nullptr,"", "", ""}}, //function pointer + }; +/** + * Visits the given ASTNode as a function. For this node only the + * traversal is preorder. + */ +/* WCS addition: Add "std::" before min and "{...}" between parenthesis in order to support min */ +void FunctionInterfaceForJIT::FormulaFormatter_visitFunction_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ) +{ + unsigned int numChildren = ASTNode_getNumChildren(node); + unsigned int n; + + // int node_type = ASTNode_getType(node); + LIBSBML_CPP_NAMESPACE::ASTNodeType_t node_type = ASTNode_getType(node); + constexpr LIBSBML_CPP_NAMESPACE::ASTNodeType_t ast_func_min = LIBSBML_CPP_NAMESPACE::AST_FUNCTION_MIN; + constexpr LIBSBML_CPP_NAMESPACE::ASTNodeType_t ast_func_max = LIBSBML_CPP_NAMESPACE::AST_FUNCTION_MAX; + if (node_type == ast_func_max || node_type == ast_func_min) { + StringBuffer_append(sb, "std::"); + } + FormulaFormatter_format(sb, node); + + // support min in sbml + StringBuffer_appendChar(sb, '('); + if (node_type == 321 || node_type == 320) { + StringBuffer_appendChar(sb, '{'); + } + + if (numChildren > 0) + { + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, LIBSBML_CPP_NAMESPACE::ASTNode_getChild(node, 0), sb ); + } + + for (n = 1; n < numChildren; n++) + { + StringBuffer_appendChar(sb, ','); + StringBuffer_appendChar(sb, ' '); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, LIBSBML_CPP_NAMESPACE::ASTNode_getChild(node, n), sb ); + } + + if (node_type == 321 || node_type == 320) { + StringBuffer_appendChar(sb, '}'); + } + StringBuffer_appendChar(sb, ')'); +} + + +/** + * Visits the given ASTNode as the function "log(10, x)" and in doing so, + * formats it as "log10(x)" (where x is any subexpression). + */ +/* WCS addition: Change function call from FormulaFormatter_visit to FormulaFormatter_visit_wcs */ +void FunctionInterfaceForJIT::FormulaFormatter_visitLog10_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ) +{ + StringBuffer_append(sb, "log10("); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs(node, ASTNode_getChild(node, 1), sb); + StringBuffer_appendChar(sb, ')'); +} + + +/** + * Visits the given ASTNode as the function "root(2, x)" and in doing so, + * formats it as "sqrt(x)" (where x is any subexpression). + */ +/* WCS addition: Change function call from FormulaFormatter_visit to FormulaFormatter_visit_wcs */ +void FunctionInterfaceForJIT::FormulaFormatter_visitSqrt_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ) +{ + StringBuffer_append(sb, "sqrt("); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs(node, ASTNode_getChild(node, 1), sb); + StringBuffer_appendChar(sb, ')'); +} + + +/** + * Visits the given ASTNode as a unary minus. For this node only the + * traversal is preorder. + */ +/* WCS addition: Change function call from FormulaFormatter_visit to FormulaFormatter_visit_wcs */ +void FunctionInterfaceForJIT::FormulaFormatter_visitUMinus_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ) +{ + StringBuffer_appendChar(sb, '-'); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs ( node, ASTNode_getLeftChild(node), sb ); +} + + +/** + * Visits the given ASTNode and continues the inorder traversal. + */ +/* WCS addition: Change function call from FormulaFormatter_visit to FormulaFormatter_visit_wcs */ +void FunctionInterfaceForJIT::FormulaFormatter_visitOther_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ) +{ + unsigned int numChildren = ASTNode_getNumChildren(node); + int group = FormulaFormatter_isGrouped(parent, node); + unsigned int n; + + + if (group) + { + StringBuffer_appendChar(sb, '('); + } + + if (numChildren == 0) { + FormulaFormatter_format(sb, node); + } + + else if (numChildren == 1) + { + //I believe this would only be called for invalid ASTNode setups, + // but this could in theory occur. This is the safest + // behavior I can think of. + FormulaFormatter_format(sb, node); + StringBuffer_appendChar(sb, '('); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, ASTNode_getChild(node, 0), sb ); + StringBuffer_appendChar(sb, ')'); + } + + else { + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, ASTNode_getChild(node, 0), sb ); + + for (n = 1; n < numChildren; n++) + { + FormulaFormatter_format(sb, node); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, ASTNode_getChild(node, n), sb ); + } + } + + if (group) + { + StringBuffer_appendChar(sb, ')'); + } +} + + + +/** + * Visits the given ASTNode node. This function is really just a + * dispatcher to either SBML_formulaToString_visitFunction() or + * SBML_formulaToString_visitOther(). + */ +/* WCS addition: Change function calls to wcs modified fuctions */ +void FunctionInterfaceForJIT::FormulaFormatter_visit_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ) +{ + if (ASTNode_isLog10(node)) + { + StringBuffer_append(sb, "log10("); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs(node, ASTNode_getChild(node, 1), sb); + StringBuffer_appendChar(sb, ')'); + // FormulaFormatter_visitLog10_wcs(parent, node, sb); + } + else if (ASTNode_isSqrt(node)) + { + StringBuffer_append(sb, "sqrt("); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs(node, ASTNode_getChild(node, 1), sb); + StringBuffer_appendChar(sb, ')'); + // FormulaFormatter_visitSqrt_wcs(parent, node, sb); + } + else if (FormulaFormatter_isFunction(node)) + { + node_map::iterator nit; + LIBSBML_CPP_NAMESPACE::ASTNodeType_t node_type = ASTNode_getType(node); + nit = nodetypes.find(node_type); + if (nit != nodetypes.cend()) { + node_structure node_st = nit -> second; + ast_function_ptr node_f_ptr = std::get<0>(node_st); + if (node_f_ptr == nullptr) { + unsigned int numChildren = ASTNode_getNumChildren(node); + unsigned int n; + std::string fist_symbol = std::get<1>(node_st); + std::string middle_symbol = std::get<2>(node_st); + std::string end_symbol = std::get<3>(node_st); + StringBuffer_append(sb, fist_symbol.c_str()); + // FormulaFormatter_format(sb, node); + if (numChildren > 0) { + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, LIBSBML_CPP_NAMESPACE::ASTNode_getChild(node, 0), sb ); + } + for (n = 1; n < numChildren; n++) { + StringBuffer_append(sb, middle_symbol.c_str()); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, LIBSBML_CPP_NAMESPACE::ASTNode_getChild(node, n), sb ); + } + StringBuffer_append(sb, end_symbol.c_str()); + } else { // if there is a function pointer + + } + } else { // user functions + unsigned int numChildren = ASTNode_getNumChildren(node); + unsigned int n; + StringBuffer_appendChar(sb, '('); + StringBuffer_append(sb, ASTNode_getName(node)); + + if (numChildren > 0) + { + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, LIBSBML_CPP_NAMESPACE::ASTNode_getChild(node, 0), sb ); + } + + for (n = 1; n < numChildren; n++) + { + StringBuffer_appendChar(sb, ','); + StringBuffer_appendChar(sb, ' '); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, LIBSBML_CPP_NAMESPACE::ASTNode_getChild(node, n), sb ); + } + StringBuffer_appendChar(sb, ')'); + } + // FormulaFormatter_visitFunction_wcs(parent, node, sb); + } + else if (ASTNode_hasTypeAndNumChildren(node, LIBSBML_CPP_NAMESPACE::AST_MINUS, 1)) + { + StringBuffer_appendChar(sb, '-'); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs ( node, ASTNode_getLeftChild(node), sb ); + // FormulaFormatter_visitUMinus_wcs(parent, node, sb); + } + else if (ASTNode_hasTypeAndNumChildren(node, LIBSBML_CPP_NAMESPACE::AST_PLUS, 1) || ASTNode_hasTypeAndNumChildren(node, LIBSBML_CPP_NAMESPACE::AST_TIMES, 1)) + { + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs (node, ASTNode_getChild(node, 0), sb); + } + else if (ASTNode_hasTypeAndNumChildren(node, LIBSBML_CPP_NAMESPACE::AST_PLUS, 0)) + { + libsbml::StringBuffer_appendInt(sb, 0); + } + else if (ASTNode_hasTypeAndNumChildren(node, LIBSBML_CPP_NAMESPACE::AST_TIMES, 0)) + { + libsbml::StringBuffer_appendInt(sb, 1); + } + else + { + FormulaFormatter_visitOther_wcs(parent, node, sb); + } +} + +/* WCS addition: Change function call from FormulaFormatter_visit to FormulaFormatter_visit_wcs */ +std::string FunctionInterfaceForJIT::SBML_formulaToString_wcs (const ASTNode_t *tree) +{ + if (tree == nullptr) { + return nullptr; + } + + LIBSBML_CPP_NAMESPACE::StringBuffer_t *buf = LIBSBML_CPP_NAMESPACE::StringBuffer_create(1024); + + FormulaFormatter_visit_wcs(NULL, tree, buf); + std::string str(LIBSBML_CPP_NAMESPACE::StringBuffer_getBuffer(buf)); + safe_free(buf); + + return str; +} +#endif // defined(WCS_HAS_SBML) + +FunctionInterfaceForJIT::~FunctionInterfaceForJIT() +{} + +/**@}*/ +} // end of namespace wcs diff --git a/src/utils/FunctionInterfaceForJIT.hpp b/src/utils/FunctionInterfaceForJIT.hpp new file mode 100644 index 0000000..f8ca522 --- /dev/null +++ b/src/utils/FunctionInterfaceForJIT.hpp @@ -0,0 +1,74 @@ +/****************************************************************************** + * * + * Copyright 2020 Lawrence Livermore National Security, LLC and other * + * Whole Cell Simulator Project Developers. See the top-level COPYRIGHT * + * file for details. * + * * + * SPDX-License-Identifier: MIT * + * * + ******************************************************************************/ + +#ifndef __WCS_UTILS_FUNCTIONINTERFACEFORJIT__ +#define __WCS_UTILS_FUNCTIONINTERFACEFORJIT__ + +#if defined(WCS_HAS_CONFIG) +#include "wcs_config.hpp" +#else +#error "no config" +#endif + +#include +#include + +#if defined(WCS_HAS_SBML) +#include +#include +#endif // defined(WCS_HAS_SBML) + + +namespace wcs { +/** \addtogroup wcs_functions + * @{ */ + +#if defined(WCS_HAS_SBML) +class FunctionInterfaceForJIT { + public: + using ASTNode_t = LIBSBML_CPP_NAMESPACE::ASTNode_t; + using strbuff_t = LIBSBML_CPP_NAMESPACE::StringBuffer_t; + public: + FunctionInterfaceForJIT(); + ~FunctionInterfaceForJIT(); + + static std::string SBML_formulaToString_wcs (const ASTNode_t *tree); + + private: + static void FormulaFormatter_visit_wcs ( const ASTNode_t*parent, + const ASTNode_t *node, + strbuff_t *sb ); + + static void FormulaFormatter_visitFunction_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ); + + static void FormulaFormatter_visitLog10_wcs ( const ASTNode_t*parent, + const ASTNode_t *node, + strbuff_t *sb ); + + static void FormulaFormatter_visitSqrt_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ); + + static void FormulaFormatter_visitUMinus_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ); + + static void FormulaFormatter_visitOther_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ); + +}; +#endif // defined(WCS_HAS_SBML) + +/**@}*/ +} // end of namespace wcs +#endif // __WCS_UTILS_FUNCTIONINTERFACEFORJIT__ diff --git a/src/utils/generate_cxx_code.cpp b/src/utils/generate_cxx_code.cpp index 66f935e..f29fdea 100644 --- a/src/utils/generate_cxx_code.cpp +++ b/src/utils/generate_cxx_code.cpp @@ -31,6 +31,7 @@ #include #include #include "wcs_types.hpp" +#include "FunctionInterfaceForJIT.hpp" #include // std::thread::hardware_concurrency #if defined(WCS_HAS_SBML) @@ -170,6 +171,61 @@ generate_cxx_code::get_all_dependencies( return dependencies_no_dupl; } +static double count2Concentration_converter ( + const LIBSBML_CPP_NAMESPACE::Model& model, + const LIBSBML_CPP_NAMESPACE::Species& species, + std::string& compartment, + double& comp_unit, + double& sub_unit) +{ + const LIBSBML_CPP_NAMESPACE::ListOfCompartments* compartment_list + = model.getListOfCompartments(); + const LIBSBML_CPP_NAMESPACE::ListOfUnitDefinitions* unit_definition_list + = model.getListOfUnitDefinitions(); + compartment = species.getCompartment(); + double compartment_size = compartment_list->get(species.getCompartment())->getSize(); + double avog_num = 6.02214e+23; + + + std::string compartment_unit ; + if (compartment_list->get(species.getCompartment())->isSetUnits()){ + compartment_unit = compartment_list->get(species.getCompartment())->getUnits(); + } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 3) { + compartment_unit = model.getVolumeUnits(); + } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 2) { + compartment_unit = model.getAreaUnits(); + } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 1) { + compartment_unit = model.getLengthUnits(); + } + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition_comp = unit_definition_list->get(compartment_unit); + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list_comp + = unit_definition_comp->getListOfUnits(); + unsigned int unitsSize_comp = unit_list_comp->size(); + // double comp_unit = 1.0; + comp_unit = 1.0; + for (unsigned int iu = 0u; iu < unitsSize_comp; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list_comp->get(iu); + comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); + } + + // Find substance units + std::string substance_unit ; + substance_unit = model.getSubstanceUnits(); + + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition_substance = unit_definition_list->get(substance_unit); + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list_substance + = unit_definition_substance->getListOfUnits(); + unsigned int unitsSize_substance = unit_list_substance->size(); + sub_unit = 1.0; + for (unsigned int iu = 0u; iu < unitsSize_substance; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list_substance->get(iu); + sub_unit = sub_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); + } + + return avog_num * compartment_size * comp_unit * sub_unit; +} + + void generate_cxx_code::return_denominators( const LIBSBML_CPP_NAMESPACE::ASTNode& formula, @@ -188,7 +244,8 @@ generate_cxx_code::return_denominators( } else { if ( std::find(denominators.cbegin(), denominators.cend(), SBML_formulaToString(denominator)) == denominators.cend()) { - denominators.push_back(SBML_formulaToString(denominator)); + std::string new_denominator = FunctionInterfaceForJIT::SBML_formulaToString_wcs(denominator); + denominators.push_back(new_denominator); } } } @@ -365,15 +422,36 @@ void generate_cxx_code::find_used_params( // Find used parameters in the rates for (unsigned int ic = 0u; ic < num_reactions; ic++) { const LIBSBML_CPP_NAMESPACE::Reaction& reaction = *(reaction_list->get(ic)); - const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list - = reaction.getKineticLaw()->getListOfLocalParameters(); - unsigned int num_localparameters = local_parameter_list->size(); - //genfile << "reaction" << reaction.getIdAttribute() <<": "; + // const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + // = reaction.getKineticLaw()->getListOfLocalParameters(); + // unsigned int num_localparameters = local_parameter_list->size(); + // //genfile << "reaction" << reaction.getIdAttribute() <<": "; using reaction_parameters_t = std::unordered_set; reaction_parameters_t lpset; - for (unsigned int pi = 0u; pi < num_localparameters; pi++) { - const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); - lpset.insert(localparameter->getIdAttribute()); + // for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + // const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); + // lpset.insert(localparameter->getIdAttribute()); + // Check SBML level for local parameters + if (model.getLevel() > 2) { + const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfLocalParameters(); + unsigned int num_localparameters = local_parameter_list->size(); + //genfile << "reaction" << reaction.getIdAttribute() <<": "; + + for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); + lpset.insert(localparameter->getIdAttribute()); + } + } else { + const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfParameters(); + unsigned int num_localparameters = local_parameter_list->size(); + //genfile << "reaction" << reaction.getIdAttribute() <<": "; + + for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + const LIBSBML_CPP_NAMESPACE::Parameter* localparameter = local_parameter_list->get(pi); + lpset.insert(localparameter->getIdAttribute()); + } } std::vector dependencies_set = get_all_dependencies(*reaction.getKineticLaw()->getMath(), @@ -1068,21 +1146,58 @@ void generate_cxx_code::print_reaction_rates( genfile << "#include \"" + header + '"' + "\n\n"; for (unsigned int ic = rid_start; ic < rid_end; ic++) { const LIBSBML_CPP_NAMESPACE::Reaction& reaction = *(reaction_list->get(ic)); - const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list - = reaction.getKineticLaw()->getListOfLocalParameters(); - unsigned int num_localparameters = local_parameter_list->size(); + // const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + // = reaction.getKineticLaw()->getListOfLocalParameters(); + // unsigned int num_localparameters = local_parameter_list->size(); + // a map to keep the compartment and their units + using compartment_t = std::unordered_map ; + typename compartment_t::const_iterator rci; + compartment_t reaction_comp; + unsigned int num_localparameters = 0u; + double model_comp_unit = 1; + double model_sub_unit = 1; + int model_comp_unit_defined = 0; + if (model.getLevel() > 2) { + const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfLocalParameters(); + num_localparameters = local_parameter_list->size(); + } else { + const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfParameters(); + num_localparameters = local_parameter_list->size(); + } + genfile << "extern \"C\" " << Real << " wcs__rate_" << reaction.getIdAttribute() << "(const std::vector<" << Real <<">& __input) {\n"; - + + double sum_stoich = 0.0; //print reaction's local parameters using reaction_local_parameters_t = std::unordered_set; reaction_local_parameters_t localpset; - for (unsigned int pi = 0u; pi < num_localparameters; pi++) { - const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); - genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " - << localparameter->getValue() << ";\n"; - localpset.insert(localparameter->getIdAttribute()); + // for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + // const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); + // genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + // << localparameter->getValue() << ";\n"; + // localpset.insert(localparameter->getIdAttribute()); + if (model.getLevel() > 2) { + const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfLocalParameters(); + for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); + genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + << localparameter->getValue() << ";\n"; + localpset.insert(localparameter->getIdAttribute()); + } + } else { + const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfParameters(); + for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + const LIBSBML_CPP_NAMESPACE::Parameter* localparameter = local_parameter_list->get(pi); + genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + << localparameter->getValue() << ";\n"; + localpset.insert(localparameter->getIdAttribute()); + } } //genfile << " int __ii=0;\n"; //genfile << "printf(\"Number of input arguments of %s received : %lu \", \"" @@ -1116,12 +1231,15 @@ void generate_cxx_code::print_reaction_rates( for (const auto& x: input_map) { var_names_ord[x.second] = x.first; } - //print first parameters in formula taking input std::vector::const_iterator it; std::vector par_names, par_names_nf; all_var_names = var_names; int par_index = 0; + const ListOfSpecies* species_list = model.getListOfSpecies(); + double avog_num = 6.02214e+23; + bool concentration_used = false; + for (it = var_names_ord.cbegin(); it < var_names_ord.cend(); it++){ rrdit = rate_rules_dep_map.find(*it); evassigit = ev_assign.find(*it); @@ -1139,10 +1257,58 @@ void generate_cxx_code::print_reaction_rates( all_var_names_it = all_var_names.find(*itf); if (all_var_names_it == all_var_names.cend()) { genfile << " " << Real << " " << *itf << " = __input[" << par_index++ << "];\n"; + if (species_list->get(*itf) != NULL) { + if (species_list->get(*itf)->isSetInitialConcentration()) { + double comp_unit; + double sub_unit; + std::string compart_name; + double convert_factor = count2Concentration_converter(model, *species_list->get(*itf), compart_name, comp_unit, sub_unit); + if (model_comp_unit_defined == 0) { + model_comp_unit = model_comp_unit * comp_unit; + model_sub_unit = model_sub_unit * sub_unit; + model_comp_unit_defined = 1; + } else { + if (model_comp_unit != comp_unit) { + WCS_THROW("Compartments with different units"); + } + } + // Add values for the compartment in map + rci = reaction_comp.find(compart_name) ; + if (rci == reaction_comp.end()) { + reaction_comp.insert(std::make_pair(compart_name, comp_unit)); + } + genfile << " " << *itf << " = " << *itf << " / " << convert_factor << ";\n"; + concentration_used = true; + } + } var_names.erase(*itf); } else { if (var_names_it != var_names.cend()) { genfile << " " << Real << " " << *itf << " = __input[" << par_index++ << "];\n"; + if (species_list->get(*itf) != NULL) { + if (species_list->get(*itf)->isSetInitialConcentration()) { + double comp_unit; + double sub_unit; + std::string compart_name; + double convert_factor = count2Concentration_converter(model, *species_list->get(*itf), compart_name, comp_unit, sub_unit); + if (model_comp_unit_defined == 0) { + model_comp_unit = model_comp_unit * comp_unit; + model_sub_unit = model_sub_unit * sub_unit; + model_comp_unit_defined = 1; + } else { + if (model_comp_unit != comp_unit) { + WCS_THROW("Compartments with different units"); + } + } + // Add values for the compartment in map + rci = reaction_comp.find(compart_name) ; + if (rci == reaction_comp.end()) { + reaction_comp.insert(std::make_pair(compart_name, comp_unit)); + } + genfile << " " << *itf << " = " << *itf << " / " << convert_factor << ";\n"; + concentration_used = true; + } + } var_names.erase(*itf); } } @@ -1157,6 +1323,33 @@ void generate_cxx_code::print_reaction_rates( genfile << " wcs_global_var." << *it << " = wcs__rate_" << *it << "(" << function_input << ");\n"; } else { genfile << " " << Real << " " << *it << " = __input[" << par_index++ << "];\n"; + if (species_list->get(*it) != NULL) { + if (species_list->get(*it)->isSetInitialConcentration()) { + double comp_unit; + double sub_unit; + std::string compart_name; + double convert_factor = count2Concentration_converter(model, *species_list->get(*it), compart_name, comp_unit, sub_unit); + if (model_comp_unit_defined == 0) { + model_comp_unit = model_comp_unit * comp_unit; + model_sub_unit = model_sub_unit * sub_unit; + model_comp_unit_defined = 1; + } else { + if (model_comp_unit != comp_unit) { + WCS_THROW("Compartments with different units"); + } + } + // Add values for the compartment in map + rci = reaction_comp.find(compart_name) ; + if (rci == reaction_comp.end()) { + reaction_comp.insert(std::make_pair(compart_name, comp_unit)); + } + genfile << " " << *it << " = " << *it << " / " << convert_factor << ";\n"; + concentration_used = true; + } + } + if (reaction.getReactant(*it) != NULL) { + sum_stoich = sum_stoich + reaction.getReactant(*it)->getStoichiometry(); + } par_names.push_back(*it); } var_names.erase(*it); @@ -1199,7 +1392,6 @@ void generate_cxx_code::print_reaction_rates( par_names_nf.push_back(x); } } - // put input parameters into maps dep_params_f.insert(std::make_pair(reaction.getIdAttribute(), par_names)); @@ -1232,8 +1424,9 @@ void generate_cxx_code::print_reaction_rates( } math = *reaction.getKineticLaw()->getMath(); update_scope_ast_node(math, wcs_all_var, wcs_all_const, localpset); + std::string new_math = FunctionInterfaceForJIT::SBML_formulaToString_wcs(&math); genfile << " " << Real << " " << reaction.getIdAttribute() << " = " - << SBML_formulaToString( &math) + << new_math //SBML_formulaToString( &math) << ";\n"; genfile << " if (!isfinite(" << reaction.getIdAttribute() << ")) {\n"; // find denominators of the dependent expressions of the reaction rate formula @@ -1289,7 +1482,16 @@ void generate_cxx_code::print_reaction_rates( genfile << " }\n"; //genfile << " printf(\" Result of %s : %f \", \"" << reaction.getIdAttribute() << "\"," // << reaction.getIdAttribute() << ");\n"; - genfile << " return " << reaction.getIdAttribute() << ";\n"; + if (!concentration_used) { + genfile << " return " << reaction.getIdAttribute() << ";\n"; + } else { + // Calculate the compartment unit for all the compartments + // double comp_unit = 1.0; + // for ( auto rci = reaction_comp.begin(); rci != reaction_comp.end(); ++rci ){ + // comp_unit = comp_unit * rci->second; + // } + genfile << " return " << reaction.getIdAttribute() << " * " << avog_num << " * " << model_comp_unit << " * " << model_sub_unit << ";\n"; + } genfile << "}\n\n"; } } @@ -1674,6 +1876,7 @@ void generate_cxx_code::write_header( << "#include \n" << "#include \n" << "#include \n" + << "#include \n" << "#include \"utils/exception.hpp\"\n\n" << "//Include the text of the SBML in here.\n" << "//Use \"xxd -i\" to get the text as a C symbol.\n" diff --git a/src/utils/graph_factory.hpp b/src/utils/graph_factory.hpp index 46e8b4f..b6a7504 100644 --- a/src/utils/graph_factory.hpp +++ b/src/utils/graph_factory.hpp @@ -303,7 +303,6 @@ GraphFactory::convert_to( { g.m_vertices[vd].m_in_edges.reserve(num_reactants); } - // Add reactants species for (unsigned int si = 0u; si < num_reactants; si++) { const auto &reactant = *(reaction.getReactant(si)); @@ -315,7 +314,8 @@ GraphFactory::convert_to( reaction_in_species.insert(s_label); if (it == smap.cend()) { - wcs::Vertex vs(*species_list->get(reactant.getSpecies()), g); + wcs::Vertex vs(model, *species_list->get(reactant.getSpecies()), g); + //WCS_THROW("TESTT"); vds = boost::add_vertex(vs, g); smap.insert(std::make_pair(s_label,vds)); } else { @@ -323,27 +323,38 @@ GraphFactory::convert_to( } std::string e_label = g[vds].get_label() + '|' + g[vd].get_label(); - eit = emap.find(e_label); + if (!species_list->get(reactant.getSpecies())->getConstant()) { //Add as a reactant + eit = emap.find(e_label); + + if (eit == emap.cend()) { + const auto ret = boost::add_edge(vds, vd, g); - if (eit == emap.cend()) { + if (!ret.second) { + WCS_THROW("Please check the reactions in your SBML file"); + return; + } + g[ret.first].set_stoichiometry_ratio(reactant.getStoichiometry()); + g[ret.first].set_label(e_label); + emap.insert(std::make_pair(e_label, ret.first)); + } else { + const auto& edge_found = eit->second; + stoic_t new_stoich = g[edge_found].get_stoichiometry_ratio() + + reactant.getStoichiometry(); + g[edge_found].set_stoichiometry_ratio(new_stoich); + } + } else { //Add as a modifier const auto ret = boost::add_edge(vds, vd, g); if (!ret.second) { WCS_THROW("Please check the reactions in your SBML file"); return; } - g[ret.first].set_stoichiometry_ratio(reactant.getStoichiometry()); + g[ret.first].set_stoichiometry_ratio(0); g[ret.first].set_label(e_label); emap.insert(std::make_pair(e_label, ret.first)); - } else { - const auto& edge_found = eit->second; - stoic_t new_stoich = g[edge_found].get_stoichiometry_ratio() - + reactant.getStoichiometry(); - g[edge_found].set_stoichiometry_ratio(new_stoich); - } + } } - // Add modifiers species for (unsigned int si = 0u; si < num_modifiers; si++) { const auto &modifier = *(reaction.getModifier(si)); @@ -355,7 +366,7 @@ GraphFactory::convert_to( reaction_in_species.insert(s_label); if (it == smap.cend()) { - wcs::Vertex vs(*species_list->get(modifier.getSpecies()), g); + wcs::Vertex vs(model, *species_list->get(modifier.getSpecies()), g); vds = boost::add_vertex(vs, g); smap.insert(std::make_pair(s_label,vds)); } else { @@ -396,7 +407,7 @@ GraphFactory::convert_to( reaction_in_species.insert(s_label); if (it == smap.cend()) { - wcs::Vertex vs(*species_list->get(x), g); + wcs::Vertex vs(model, *species_list->get(x), g); vds = boost::add_vertex(vs, g); smap.insert(std::make_pair(s_label,vds)); } else { @@ -434,7 +445,7 @@ GraphFactory::convert_to( if (risit == reaction_in_species.cend()) { reaction_in_species.insert(x); if (it == smap.cend()) { - wcs::Vertex vs(*species_list->get(x), g); + wcs::Vertex vs(model, *species_list->get(x), g); vds = boost::add_vertex(vs, g); smap.insert(std::make_pair(x,vds)); } else { @@ -472,7 +483,7 @@ GraphFactory::convert_to( risit = reaction_in_species.find(x); if (risit == reaction_in_species.cend()) { if (it == smap.cend()) { - wcs::Vertex vs(*species_list->get(x), g); + wcs::Vertex vs(model, *species_list->get(x), g); vds = boost::add_vertex(vs, g); smap.insert(std::make_pair(x,vds)); } else { @@ -506,7 +517,7 @@ GraphFactory::convert_to( v_new_desc_t vds; if (it == smap.cend()) { - wcs::Vertex vs(*species_list->get(product.getSpecies()), g); + wcs::Vertex vs(model, *species_list->get(product.getSpecies()), g); vds = boost::add_vertex(vs, g); smap.insert(std::make_pair(s_label,vds)); } else { @@ -514,22 +525,35 @@ GraphFactory::convert_to( } std::string e_label = g[vd].get_label() + '|' + g[vds].get_label(); - eit = emap.find(e_label); - if (eit == emap.cend()) { - const auto ret = boost::add_edge(vd, vds, g); + + if (!species_list->get(product.getSpecies())->getConstant()) { //Add as a product + eit = emap.find(e_label); + if (eit == emap.cend()) { + const auto ret = boost::add_edge(vd, vds, g); + + if (!ret.second) { + WCS_THROW("Please check the reactions in your SBML file"); + return; + } + g[ret.first].set_stoichiometry_ratio(product.getStoichiometry()); + g[ret.first].set_label( g[vd].get_label() + '|' + g[vds].get_label()); + emap.insert(std::make_pair(e_label, ret.first)); + } else { + const auto& edge_found = eit->second; + stoic_t new_stoich = g[edge_found].get_stoichiometry_ratio() + + product.getStoichiometry(); + g[edge_found].set_stoichiometry_ratio(new_stoich); + } + } else { //Add as a modifier + const auto ret = boost::add_edge(vds, vd, g); if (!ret.second) { WCS_THROW("Please check the reactions in your SBML file"); return; } - g[ret.first].set_stoichiometry_ratio(product.getStoichiometry()); - g[ret.first].set_label( g[vd].get_label() + '|' + g[vds].get_label()); + g[ret.first].set_stoichiometry_ratio(0); + g[ret.first].set_label(e_label); emap.insert(std::make_pair(e_label, ret.first)); - } else { - const auto& edge_found = eit->second; - stoic_t new_stoich = g[edge_found].get_stoichiometry_ratio() - + product.getStoichiometry(); - g[edge_found].set_stoichiometry_ratio(new_stoich); } } } diff --git a/src/utils/sbml_utils.cpp b/src/utils/sbml_utils.cpp index 9efdba7..7599cca 100644 --- a/src/utils/sbml_utils.cpp +++ b/src/utils/sbml_utils.cpp @@ -81,12 +81,22 @@ sbml_utils::find_undeclared_species_in_reaction_formula( } // create an unordered_set for reaction local parameters - const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list - = reaction.getKineticLaw()->getListOfLocalParameters(); - unsigned int num_local_parameters = local_parameter_list->size(); + if (model.getLevel() > 2) { + const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfLocalParameters(); + unsigned int num_local_parameters = local_parameter_list->size(); - for (unsigned int pi = 0u; pi < num_local_parameters; pi++) { - lpset.insert(local_parameter_list->get(pi)->getIdAttribute()); + for (unsigned int pi = 0u; pi < num_local_parameters; pi++) { + lpset.insert(local_parameter_list->get(pi)->getIdAttribute()); + } + } else { + const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfParameters(); + unsigned int num_local_parameters = local_parameter_list->size(); + + for (unsigned int pi = 0u; pi < num_local_parameters; pi++) { + lpset.insert(local_parameter_list->get(pi)->getIdAttribute()); + } } @@ -94,7 +104,11 @@ sbml_utils::find_undeclared_species_in_reaction_formula( unsigned int num_reactants = reaction.getNumReactants(); for (unsigned int ri = 0u; ri < num_reactants; ri++) { const auto &reactant = *(reaction.getReactant(ri)); - rset.insert(species_list->get(reactant.getSpecies())->getIdAttribute()); + if (!species_list->get(reactant.getSpecies())->getConstant()) { //Add as a reactant + rset.insert(species_list->get(reactant.getSpecies())->getIdAttribute()); + } else { + mset.insert(species_list->get(reactant.getSpecies())->getIdAttribute()); + } } // create an unordered set for reaction_modifiers @@ -174,7 +188,9 @@ sbml_utils::get_reaction_parameters( for (unsigned int ri = 0u; ri < num_reactants; ri++) { const auto &reactant = *(reaction.getReactant(ri)); - rset.insert(species_list->get(reactant.getSpecies())->getIdAttribute()); + if (!species_list->get(reactant.getSpecies())->getConstant()) { //Add as a reactant + rset.insert(species_list->get(reactant.getSpecies())->getIdAttribute()); + } } // create an unordered_set for model compartments