105 if ((pParm =
getParameter(
"HYBRID.LowerStochLimit")) != NULL)
111 if ((pParm =
getParameter(
"HYBRID.UpperStochLimit")) != NULL)
117 if ((pParm =
getParameter(
"HYBRID.RungeKuttaStepsize")) != NULL)
123 if ((pParm =
getParameter(
"HYBRID.PartitioningInterval")) != NULL)
206 for (i = 0; ((i <
mMaxSteps) && (time < endTime)); i++)
216 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.");
342 while ((dt - integrationTime) >
mStepsize)
354 std::copy(dependents.begin(), dependents.end(),
373 while ((dt - integrationTime) >
mStepsize)
397 std::copy(dependents.begin(), dependents.end(),
498 (*mpReactions)[j->
mIndex]->calculateParticleFlux();
512 deriv[i] += bal * (*mpReactions)[j->
mIndex]->getParticleFlux();
535 target[i] = (*mpMetabolites)[i]->getValue();
555 (*mpMetabolites)[i]->setValue(source[i]);
556 (*mpMetabolites)[i]->refreshConcentration();
612 std::copy(dependents.begin(), dependents.end(),
629 std::set <size_t>::iterator iter, iterEnd;
680 if (
mAmu[rIndex] == 0)
return std::numeric_limits<C_FLOAT64>::infinity();
683 return - 1 * log(rand2) /
mAmu[rIndex];
695 mAmu[rIndex] = (*mpReactions)[rIndex]->calculateParticleFlux();
718 for (
size_t i = 0; i < substrates.size(); i++)
720 num_ident = substrates[i].mMultiplicity;
726 lower_bound = number - num_ident;
727 substrate_factor = substrate_factor * pow((
double) number, (
int) num_ident - 1);
731 while (number > lower_bound)
739 if ((amu == 0) || (substrate_factor == 0))
748 C_FLOAT64 rate_factor = (*mpReactions)[rIndex]->calculateParticleFlux();
752 amu *= rate_factor / substrate_factor;;
757 mAmu[rIndex] = rate_factor;
780 if (
mAmu[rIndex] != 0.0)
808 size_t i, numReactions = mpReactions->
size();
815 for (i = 0; i < numReactions; i++)
818 if ((*mpReactions)[i]->getCompartmentNumber() != 1)
return - 1;
821 if ((*mpReactions)[i]->isReversible() != 0)
return - 2;
827 for (j = 0; j < mStoi.
numRows(); j++)
829 multFloat = mStoi[j][i];
830 multInt =
static_cast<C_INT32>(floor(multFloat + 0.5));
832 if ((multFloat - multInt) >
INT_EPSILON)
return - 3;
858 for (i = 0; i < numReactions; i++)
861 &(*mpReactions)[i]->getChemEq().getBalances();
863 for (j = 0; j < balances->
size(); j++)
868 newElement.
mMultiplicity =
static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity() + 0.5));
878 balances = &(*mpReactions)[i]->getChemEq().getSubstrates();
880 for (j = 0; j < balances->
size(); j++)
885 newElement.
mMultiplicity =
static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity() + 0.5));
903 std::vector< std::set<std::string>* > DependsOn;
904 std::vector< std::set<std::string>* > Affects;
909 for (i = 0; i < numReactions; i++)
923 for (i = 0; i < numReactions; i++)
925 for (j = 0; j < numReactions; j++)
931 std::set<std::string>::iterator iter = Affects[i]->begin();
933 for (; iter != Affects[i]->end(); iter++)
935 if (DependsOn[j]->count(*iter))
950 for (i = 0; i < numReactions; i++)
967 size_t metaboliteIndex;
995 std::vector< std::set<size_t>* > participatesIn;
1003 for (i = 0; i < numReactions; i++)
1009 for (i = 0; i < numReactions; i++)
1012 std::set<size_t>::iterator iter = participatesIn[i]->begin();
1014 for (; iter != participatesIn[i]->end(); iter++)
1018 for (i = 0; i < numReactions; i++)
1020 delete participatesIn[i];
1071 (*mpMetabolites)[i]->refreshConcentration();
1122 if (prevFlag != NULL)
1165 std::set <size_t>::iterator iter, iterEnd;
1192 (*mpMetabolites)[i]->refreshConcentration();
1300 std::set<std::string> *retset =
new std::set<std::string>;
1302 size_t i, imax = (*mpReactions)[rIndex]->getFunctionParameters().size();
1305 for (i = 0; i < imax; ++i)
1311 const std::vector <std::string> & metabKeylist =
1312 (*mpReactions)[rIndex]->getParameterMappings()[i];
1313 jmax = metabKeylist.size();
1315 for (j = 0; j < jmax; ++j)
1317 retset->insert(metabKeylist[j]);
1334 std::set<std::string> *retset =
new std::set<std::string>;
1343 retset->insert(
mLocalBalances[rIndex][i].mpMetabolite->getKey());
1359 std::set<size_t> *retset =
new std::set<size_t>;
1368 static unsigned C_INT32 counter = 0;
1382 os << (*mpMetabolites)[i]->getValue() <<
" ";
1395 os << (*mpMetabolites)[i]->getValue() <<
" ";
1414 std::set <size_t>::iterator iter, iterEnd;
1416 os <<
"outputDebug(" << level <<
") *********************************************** BEGIN" << std::endl;
1424 os <<
"mMaxSteps: " <<
mMaxSteps << std::endl;
1425 os <<
"mMaxBalance: " <<
mMaxBalance << std::endl;
1437 os <<
"mStoi: " << std::endl;
1442 os <<
mStoi[i][j] <<
" ";
1450 os <<
temp[i] <<
" ";
1453 os <<
"curentState: ";
1483 os <<
"mReactionFlags: " << std::endl;
1488 os <<
"mFirstReactionFlag: " << std::endl;
1491 os <<
"mMetabFlags: " << std::endl;
1502 os <<
"mLocalBalances: " << std::endl;
1512 os <<
"mLocalSubstrates: " << std::endl;
1527 os <<
"mStepsize: " <<
mStepsize << std::endl;
1528 os <<
"mMetab2React: " << std::endl;
1540 os <<
"mAmu: " << std::endl;
1543 os <<
mAmu[i] <<
" ";
1546 os <<
"mAmuOld: " << std::endl;
1552 os <<
"mUpdateSet: " << std::endl;
1559 os <<
"mDG: " << std::endl <<
mDG;
1560 os <<
"mPQ: " << std::endl <<
mPQ;
1561 os <<
"Particle numbers: " << std::endl;
1565 os << (*mpMetabolites)[i]->getValue() <<
" ";
1595 os <<
temp[i] <<
" ";
1598 os <<
"currentState: ";
1628 os <<
"mReactionFlags: " << std::endl;
1633 os <<
"mFirstReactionFlag: " << std::endl;
1634 if (mFirstReactionFlag == NULL) os <<
"NULL" << std::endl;
else os << *
mFirstReactionFlag;
1636 os <<
"mMetabFlags: " << std::endl;
1647 os <<
"mAmu: " << std::endl;
1650 os <<
mAmu[i] <<
" ";
1653 os <<
"mAmuOld: " << std::endl;
1659 os <<
"mUpdateSet: " << std::endl;
1665 os <<
"mPQ: " << std::endl <<
mPQ;
1666 os <<
"Particle numbers: " << std::endl;
1670 os << (*mpMetabolites)[i]->getValue() <<
" ";
1683 os <<
"outputDebug(" << level <<
") ************************************************* END" << std::endl;
1689 os <<
"CHybridStochFlag " << std::endl;
1690 os <<
" mIndex: " << d.
mIndex <<
" mValue: " << d.
mValue << std::endl;
1693 os <<
" prevIndex: " << d.
mpPrev->
mIndex <<
" prevPointer: " << d.
mpPrev << std::endl;
1695 os <<
" prevPointer: NULL" << std::endl;
1698 os <<
" nextIndex: " << d.
mpNext->
mIndex <<
" nextPointer: " << d.
mpNext << std::endl;
1700 os <<
" nextPointer: NULL" << std::endl;
1707 os <<
"CHybridBalance" << std::endl;
1737 for (i = 0; i < imax; ++i)
1757 for (i = 0; i < imax; ++i)
1777 for (i = 0; i < imax; ++i)
1817 if (*
getValue(
"Runge Kutta Stepsize").pDOUBLE <= 0.0)
1848 for (i = 0; i < imax; ++i)
1860 for (i = 0; i < imax; ++i)
1872 for (i = 0; i < imax; ++i)
CRandom * mpRandomGenerator
C_FLOAT64 mLowerStochLimit
std::ostream & operator<<(std::ostream &os, const CHybridStochFlag &d)
unsigned C_INT32 mMaxSteps
CCopasiVectorN< CEvent > & getEvents()
void outputData(std::ostream &os, C_INT32 mode)
void getState(std::vector< C_FLOAT64 > &target)
void getStochTimeAndIndex(C_FLOAT64 &ds, size_t &rIndex)
const CCopasiVector< CMetab > & getMetabolites() const
const std::string & getObjectName() const
const CCopasiVectorN< CModelValue > & getModelValues() const
virtual C_FLOAT64 getRandomOO()
virtual size_t size() const
void addDependent(const size_t &node, const size_t &dependent)
void initMethod(C_FLOAT64 time)
void calculateDerivative(std::vector< C_FLOAT64 > &deriv)
static C_INT32 checkModel(CModel *model)
size_t getNumMetabs() const
std::vector< C_FLOAT64 > k3
std::vector< metabStatus > mMetabFlags
void updateSimulatedValues(const bool &updateMoieties)
virtual size_t numRows() const
#define UPPER_STOCH_LIMIT
CHybridStochFlag * mFirstReactionFlag
void setState(std::vector< C_FLOAT64 > &source)
const std::set< size_t > & getDependents(const size_t &node) const
CTrajectoryProblem * mpProblem
#define PARTITIONING_INTERVAL
std::set< size_t > * getParticipatesIn(size_t rIndex)
virtual void setValue(const C_FLOAT64 &value)
size_t getIndex(const CModelEntity *entity) const
void setupMetab2ReactComplete()
CMatrix< C_FLOAT64 > mStoi
CIndexedPriorityQueue mPQ
C_FLOAT64 getKey(const size_t index) const
void setupMetab2ReactPlusModifier()
virtual bool isValidProblem(const CCopasiProblem *pProblem)
virtual Status step(const double &deltaT)
void setTime(const C_FLOAT64 &time)
CHybridStochFlag * mpPrev
void initializeParameter()
size_t getTotSteps() const
void updateTauMu(size_t rIndex, C_FLOAT64 time)
void integrateDeterministicPart(C_FLOAT64 ds)
void removeDeterministicReaction(size_t rIndex)
void rungeKutta(C_FLOAT64 dt)
bool removeParameter(const std::string &name)
const CCopasiVectorNS< CReaction > * mpReactions
#define MCTrajectoryMethod
std::vector< C_FLOAT64 > k4
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
virtual C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)=0
size_t getNumVariableMetabs() const
size_t getNumDependentReactionMetabs() const
void outputDebug(std::ostream &os, size_t level)
void setState(const CState &state)
size_t insertStochReaction(const size_t index, const C_FLOAT64 key)
const C_FLOAT64 & getDuration() const
CHybridMethod(const CHybridMethod &src, const CCopasiContainer *pParent=NULL)
const Value & getValue() const
virtual bool elevateChildren()
C_FLOAT64 generateReactionTime(size_t rIndex)
bool setValue(const std::string &name, const CType &value)
CCopasiParameter * getParameter(const std::string &name)
std::set< std::string > * getAffects(size_t rIndex)
CHybridStochFlag * mpNext
size_t getNumIndependentReactionMetabs() const
size_t removeStochReaction(const size_t index)
CCopasiVector< CMetab > * mpMetabolites
void calculateAmu(size_t rIndex)
std::vector< C_FLOAT64 > mAmu
size_t mStepsAfterPartitionSystem
void updatePriorityQueue(size_t rIndex, C_FLOAT64 time)
CCopasiVectorNS< CCompartment > & getCompartments()
void updateNode(const size_t index, const C_FLOAT64 key)
const C_FLOAT64 & getValue() const
unsigned C_INT32 mPartitioningInterval
const C_FLOAT64 & getTime() const
virtual size_t getIndex(const CCopasiObject *pObject) const
void setupPriorityQueue(C_FLOAT64 startTime=0.0)
virtual bool isValidProblem(const CCopasiProblem *pProblem)
std::vector< C_FLOAT64 > currentState
std::vector< C_FLOAT64 > k2
void insertDeterministicReaction(size_t rIndex)
const CStateTemplate & getStateTemplate() const
std::vector< std::vector< CHybridBalance > > mLocalSubstrates
unsigned C_INT32 mRandomSeed
static CHybridMethod * createHybridMethod()
size_t getNumModelValues() const
const ModelType & getModelType() const
CCopasiVectorNS< CReaction > & getReactions()
std::vector< std::vector< CHybridBalance > > mLocalBalances
std::vector< std::set< size_t > > mMetab2React
void integrateDeterministicPartEuler(C_FLOAT64 ds)
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
virtual void start(const CState *initialState)
C_FLOAT64 mUpperStochLimit
std::vector< C_FLOAT64 > k1
#define RUNGE_KUTTA_STEPSIZE
const CModelEntity::Status & getStatus() const
const CCopasiVector< CMetab > & getMetabolitesX() const
const CMatrix< C_FLOAT64 > & getStoi() const
std::set< size_t > mUpdateSet
std::vector< C_FLOAT64 > mAmuOld
virtual size_t numCols() const
std::vector< C_FLOAT64 > temp
void fireReaction(size_t rIndex)
std::vector< CHybridStochFlag > mReactionFlags
std::string suitableForStochasticSimulation() const
CModel * getModel() const
std::set< std::string > * getDependsOn(size_t rIndex)
size_t mNumVariableMetabs
void initializeIndexPointer(const size_t numberOfReactions)
#define LOWER_STOCH_LIMIT
C_FLOAT64 * beginIndependent()
static bool modelHasAssignments(const CModel *pModel)
virtual void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
void setupDependencyGraph()