113 if ((pParm =
getParameter(
"HYBRID.LowerStochLimit")) != NULL)
119 if ((pParm =
getParameter(
"HYBRID.UpperStochLimit")) != NULL)
125 if ((pParm =
getParameter(
"HYBRID.PartitioningInterval")) != NULL)
224 for (i = 0; ((i <
mMaxSteps) && (time < endTime)); i++)
234 CCopasiMessage(
CCopasiMessage::WARNING,
"maximum number of reaction events was reached in at least one simulation step.\nThat means time intervals in the output may not be what you requested.");
272 mFirstMetabIndex = mpModel->getStateTemplate().getIndex(mpModel->getMetabolitesX()[0]);
364 if (!pModel)
return 1.0e009;
371 for (i = 0, imax = Compartment.
size(); i < imax; i++)
372 if (Compartment[i]->
getValue() < Volume) Volume = Compartment[i]->
getValue();
456 tdist = fabs(deltaT);
457 d__1 = fabs(
mTime), d__2 = fabs(EndTime);
460 if (tdist < std::numeric_limits< C_FLOAT64 >::epsilon() * 2. * w0)
526 std::copy(dependents.begin(), dependents.end(),
551 (*mpReactions)[j->
mIndex]->calculateParticleFlux();
565 deriv[i] += bal * (*mpReactions)[j->
mIndex]->getParticleFlux();
625 std::copy(dependents.begin(), dependents.end(),
645 std::set <size_t>::iterator iter, iterEnd;
695 if (
mAmu[rIndex] == 0)
return std::numeric_limits<C_FLOAT64>::infinity();
698 return - 1 * log(rand2) /
mAmu[rIndex];
710 mAmu[rIndex] = (*mpReactions)[rIndex]->calculateParticleFlux();
729 const std::vector<CHybridLSODABalance> & substrates =
mLocalSubstrates[rIndex];
733 for (
size_t i = 0; i < substrates.size(); i++)
735 num_ident = substrates[i].mMultiplicity;
740 number =
static_cast<C_INT32>((*mpMetabolites)[substrates[i].mIndex]->getValue());
741 lower_bound = number - num_ident;
742 substrate_factor = substrate_factor * pow((
double) number, (
int) num_ident - 1);
746 while (number > lower_bound)
754 if ((amu == 0) || (substrate_factor == 0))
763 C_FLOAT64 rate_factor = (*mpReactions)[rIndex]->calculateParticleFlux();
767 amu *= rate_factor / substrate_factor;;
772 mAmu[rIndex] = rate_factor;
795 if (
mAmu[rIndex] != 0.0)
824 size_t i, j, numReactions = mpReactions->
size();
828 for (i = 0; i < numReactions; i++)
831 if ((*mpReactions)[i]->getCompartmentNumber() != 1)
return - 1;
834 if ((*mpReactions)[i]->isReversible() != 0)
return - 2;
840 for (j = 0; j < mStoi.
numRows(); j++)
842 multFloat = mStoi[j][i];
843 multInt =
static_cast<C_INT32>(floor(multFloat + 0.5));
845 if ((multFloat - multInt) >
INT_EPSILON)
return - 3;
871 for (i = 0; i < numReactions; i++)
874 &(*mpReactions)[i]->getChemEq().getBalances();
876 for (j = 0; j < balances->
size(); j++)
881 newElement.
mMultiplicity =
static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity() + 0.5));
891 balances = &(*mpReactions)[i]->getChemEq().getSubstrates();
893 for (j = 0; j < balances->
size(); j++)
898 newElement.
mMultiplicity =
static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity() + 0.5));
916 std::vector< std::set<std::string>* > DependsOn;
917 std::vector< std::set<std::string>* > Affects;
922 for (i = 0; i < numReactions; i++)
936 for (i = 0; i < numReactions; i++)
938 for (j = 0; j < numReactions; j++)
944 std::set<std::string>::iterator iter = Affects[i]->begin();
946 for (; iter != Affects[i]->end(); iter++)
948 if (DependsOn[j]->count(*iter))
963 for (i = 0; i < numReactions; i++)
980 size_t metaboliteIndex;
1008 std::vector< std::set<size_t>* > participatesIn;
1016 for (i = 0; i < numReactions; i++)
1022 for (i = 0; i < numReactions; i++)
1025 std::set<size_t>::iterator iter = participatesIn[i]->begin();
1027 for (; iter != participatesIn[i]->end(); iter++)
1031 for (i = 0; i < numReactions; i++)
1033 delete participatesIn[i];
1084 (*mpMetabolites)[i]->refreshConcentration();
1135 if (prevFlag != NULL)
1178 std::set <size_t>::iterator iter, iterEnd;
1207 (*mpMetabolites)[i]->refreshConcentration();
1318 std::set<std::string> *retset =
new std::set<std::string>;
1320 size_t i, imax = (*mpReactions)[rIndex]->getFunctionParameters().size();
1323 for (i = 0; i < imax; ++i)
1329 const std::vector <std::string> & metabKeylist =
1330 (*mpReactions)[rIndex]->getParameterMappings()[i];
1331 jmax = metabKeylist.size();
1333 for (j = 0; j < jmax; ++j)
1335 retset->insert(metabKeylist[j]);
1352 std::set<std::string> *retset =
new std::set<std::string>;
1361 retset->insert(
mLocalBalances[rIndex][i].mpMetabolite->getKey());
1377 std::set<size_t> *retset =
new std::set<size_t>;
1386 static unsigned C_INT32 counter = 0;
1400 os << (*mpMetabolites)[i]->getValue() <<
" ";
1413 os << (*mpMetabolites)[i]->getValue() <<
" ";
1432 std::set< size_t >::iterator iter, iterEnd;
1434 os <<
"outputDebug(" << level <<
") *********************************************** BEGIN" << std::endl;
1442 os <<
"mMaxSteps: " <<
mMaxSteps << std::endl;
1443 os <<
"mMaxBalance: " <<
mMaxBalance << std::endl;
1455 os <<
"mStoi: " << std::endl;
1460 os <<
mStoi[i][j] <<
" ";
1468 os <<
temp[i] <<
" ";
1472 os <<
"mReactionFlags: " << std::endl;
1477 os <<
"mFirstReactionFlag: " << std::endl;
1480 os <<
"mMetabFlags: " << std::endl;
1491 os <<
"mLocalBalances: " << std::endl;
1501 os <<
"mLocalSubstrates: " << std::endl;
1516 os <<
"mStepsize: " <<
mStepsize << std::endl;
1517 os <<
"mMetab2React: " << std::endl;
1529 os <<
"mAmu: " << std::endl;
1532 os <<
mAmu[i] <<
" ";
1535 os <<
"mAmuOld: " << std::endl;
1541 os <<
"mUpdateSet: " << std::endl;
1548 os <<
"mDG: " << std::endl <<
mDG;
1549 os <<
"mPQ: " << std::endl <<
mPQ;
1550 os <<
"Particle numbers: " << std::endl;
1554 os << (*mpMetabolites)[i]->getValue() <<
" ";
1584 os <<
temp[i] <<
" ";
1588 os <<
"mReactionFlags: " << std::endl;
1593 os <<
"mFirstReactionFlag: " << std::endl;
1594 if (mFirstReactionFlag == NULL) os <<
"NULL" << std::endl;
else os << *
mFirstReactionFlag;
1596 os <<
"mMetabFlags: " << std::endl;
1607 os <<
"mAmu: " << std::endl;
1610 os <<
mAmu[i] <<
" ";
1613 os <<
"mAmuOld: " << std::endl;
1619 os <<
"mUpdateSet: " << std::endl;
1625 os <<
"mPQ: " << std::endl <<
mPQ;
1626 os <<
"Particle numbers: " << std::endl;
1630 os << (*mpMetabolites)[i]->getValue() <<
" ";
1643 os <<
"outputDebug(" << level <<
") ************************************************* END" << std::endl;
1649 os <<
"CHybridLSODAStochFlag " << std::endl;
1650 os <<
" mIndex: " << d.
mIndex <<
" mValue: " << d.
mValue << std::endl;
1653 os <<
" prevIndex: " << d.
mpPrev->
mIndex <<
" prevPointer: " << d.
mpPrev << std::endl;
1655 os <<
" prevPointer: NULL" << std::endl;
1658 os <<
" nextIndex: " << d.
mpNext->
mIndex <<
" nextPointer: " << d.
mpNext << std::endl;
1660 os <<
" nextPointer: NULL" << std::endl;
1667 os <<
"CHybridLSODABalance" << std::endl;
1690 for (i = 0; i < imax; ++i)
1702 for (i = 0; i < imax; ++i)
1714 for (i = 0; i < imax; ++i)
1754 for (i = 0; i < imax; ++i)
1766 for (i = 0; i < imax; ++i)
1778 for (i = 0; i < imax; ++i)
C_FLOAT64 getDefaultAtol(const CModel *pModel) const
CVector< C_FLOAT64 > mDWork
void outputData(std::ostream &os, C_INT32 mode)
void getStochTimeAndIndex(C_FLOAT64 &ds, size_t &rIndex)
CCopasiVectorN< CEvent > & getEvents()
unsigned C_INT32 mPartitioningInterval
std::vector< std::set< size_t > > mMetab2React
std::vector< C_FLOAT64 > mAmu
const CCopasiVector< CMetab > & getMetabolites() const
const std::string & getObjectName() const
unsigned C_INT32 mMaxSteps
CIndexedPriorityQueue mPQ
const CCopasiVectorN< CModelValue > & getModelValues() const
virtual C_FLOAT64 getRandomOO()
virtual size_t size() const
void setupDependencyGraph()
void addDependent(const size_t &node, const size_t &dependent)
static C_INT32 checkModel(CModel *model)
void updateTauMu(size_t rIndex, C_FLOAT64 time)
CRandom * mpRandomGenerator
size_t getNumMetabs() const
unsigned C_INT32 mStepsAfterPartitionSystem
void setOstream(std::ostream &os)
void setupMetab2ReactComplete()
void updateSimulatedValues(const bool &updateMoieties)
virtual size_t numRows() const
virtual bool elevateChildren()
virtual bool isValidProblem(const CCopasiProblem *pProblem)
#define UPPER_STOCH_LIMIT
static bool modelHasAssignments(const CModel *pModel)
C_FLOAT64 mUpperStochLimit
const std::set< size_t > & getDependents(const size_t &node) const
CTrajectoryProblem * mpProblem
#define PARTITIONING_INTERVAL
void resize(size_t size, const bool ©=false)
virtual void setValue(const C_FLOAT64 &value)
C_FLOAT64 generateReactionTime(size_t rIndex)
CHybridLSODAStochFlag * mpNext
C_FLOAT64 getKey(const size_t index) const
std::ostringstream mErrorMsg
virtual bool isValidProblem(const CCopasiProblem *pProblem)
void updatePriorityQueue(size_t rIndex, C_FLOAT64 time)
unsigned C_INT32 mRandomSeed
void setTime(const C_FLOAT64 &time)
CCopasiVector< CMetab > * mpMetabolites
CHybridMethodLSODA * pMethod
CMatrix< C_FLOAT64 > mStoi
static CHybridMethodLSODA * createHybridMethodLSODA()
void setupPriorityQueue(C_FLOAT64 startTime=0.0)
size_t mNumVariableMetabs
bool removeParameter(const std::string &name)
#define MCTrajectoryMethod
virtual Status step(const double &deltaT)
void integrateDeterministicPart(C_FLOAT64 ds)
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
const C_FLOAT64 & getQuantity2NumberFactor() const
size_t getNumVariableMetabs() const
size_t getNumDependentReactionMetabs() const
void setState(const CState &state)
void setupMetab2ReactPlusModifier()
size_t insertStochReaction(const size_t index, const C_FLOAT64 key)
const C_FLOAT64 & getDuration() const
CVector< C_FLOAT64 > mYdot
std::vector< std::vector< CHybridLSODABalance > > mLocalBalances
const Value & getValue() const
#define PARTITIONING_STEPSIZE
virtual C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)=0
C_INT32 mMaxIntBeforeStep
bool setValue(const std::string &name, const CType &value)
std::vector< C_FLOAT64 > temp
CCopasiParameter * getParameter(const std::string &name)
static void EvalF(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
CHybridMethodLSODA(const CHybridMethodLSODA &src, const CCopasiContainer *pParent=NULL)
std::set< size_t > * getParticipatesIn(size_t rIndex)
size_t getNumIndependentReactionMetabs() const
std::vector< CHybridLSODAStochFlag > mReactionFlags
size_t removeStochReaction(const size_t index)
const C_FLOAT64 & getTime() const
void removeDeterministicReaction(size_t rIndex)
void evalF(const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
void setTime(const C_FLOAT64 &time)
std::set< std::string > * getDependsOn(size_t rIndex)
CCopasiVectorNS< CCompartment > & getCompartments()
void updateNode(const size_t index, const C_FLOAT64 key)
const C_FLOAT64 & getValue() const
virtual void start(const CState *initialState)
const C_FLOAT64 & getTime() const
virtual size_t getIndex(const CCopasiObject *pObject) const
std::ostream & operator<<(std::ostream &os, const CHybridLSODAStochFlag &d)
void initMethod(C_FLOAT64 time)
void outputDebug(std::ostream &os, size_t level)
CHybridLSODAStochFlag * mpPrev
size_t getNumModelValues() const
bool addParameter(const CCopasiParameter ¶meter)
CCopasiVectorNS< CReaction > & getReactions()
const CCopasiVectorNS< CReaction > * mpReactions
std::set< std::string > * getAffects(size_t rIndex)
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
std::vector< metabStatus > mMetabFlags
const CState & getState() const
const CModelEntity::Status & getStatus() const
const CCopasiVector< CMetab > & getMetabolitesX() const
const CMatrix< C_FLOAT64 > & getStoi() const
virtual size_t numCols() const
C_FLOAT64 mLowerStochLimit
std::string suitableForStochasticSimulation() const
CModel * getModel() const
void initializeParameter()
void insertDeterministicReaction(size_t rIndex)
std::vector< std::vector< CHybridLSODABalance > > mLocalSubstrates
std::vector< C_FLOAT64 > mAmuOld
void calculateDerivative(std::vector< C_FLOAT64 > &deriv)
CHybridLSODAStochFlag * mFirstReactionFlag
void initializeIndexPointer(const size_t numberOfReactions)
#define LOWER_STOCH_LIMIT
void fireReaction(size_t rIndex)
C_FLOAT64 * beginIndependent()
virtual void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
std::set< size_t > mUpdateSet
void calculateAmu(size_t rIndex)