Skip to content

Commit

Permalink
Bug 3240: Automatic step size implemented in Hybrid LSODA & Runge-Kutta.
Browse files Browse the repository at this point in the history
  • Loading branch information
@shoops
shoops committed May 23, 2024
1 parent b73a7c4 commit 8ea6faa
Show file tree
Hide file tree
Showing 4 changed files with 113 additions and 27 deletions.
101 changes: 89 additions & 12 deletions copasi/trajectory/CHybridMethod.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// Copyright (C) 2019 - 2023 by Pedro Mendes, Rector and Visitors of the
// Copyright (C) 2019 - 2024 by Pedro Mendes, Rector and Visitors of the
// University of Virginia, University of Heidelberg, and University
// of Connecticut School of Medicine.
// All rights reserved.
Expand Down Expand Up @@ -68,15 +68,71 @@
*/
CHybridMethod::CHybridMethod(const CDataContainer * pParent,
const CTaskEnum::Method & methodType,
const CTaskEnum::Task & taskType):
CTrajectoryMethod(pParent, methodType, taskType)
const CTaskEnum::Task & taskType)
: CTrajectoryMethod(pParent, methodType, taskType)
, mNumVariableMetabs(C_INVALID_INDEX)
, mFirstMetabIndex(C_INVALID_INDEX)
, mpFirstMetabValue(nullptr)
, mMaxSteps(MAX_STEPS)
, mMaxStepsReached(false)
, mUseRandomSeed(USE_RANDOM_SEED)
, mRandomSeed(RANDOM_SEED)
, mMaxBalance()
, mMaxIntBeforeStep()
, mReactions()
, mCurrentState()
, mSpeciesRates()
, mRateOffset(C_INVALID_INDEX)
, mUpdateSequences()
, mReactionFlags()
, mFirstReactionFlag(nullptr)
, mMetabFlags()
, mLowerStochLimit(LOWER_STOCH_LIMIT)
, mUpperStochLimit(UPPER_STOCH_LIMIT)
, mPartitioningInterval(PARTITIONING_INTERVAL)
, mStepsAfterPartitionSystem(0)
, mMetab2React()
, mAmu()
, mAmuOld()
, mpRandomGenerator(nullptr)
, mDG()
, mPQ()
, mAutomaticStepSize(false)
{
initializeParameter();
}

CHybridMethod::CHybridMethod(const CHybridMethod & src,
const CDataContainer * pParent):
CTrajectoryMethod(src, pParent)
const CDataContainer * pParent)
: CTrajectoryMethod(src, pParent)
, mNumVariableMetabs(src.mNumVariableMetabs)
, mFirstMetabIndex(src.mFirstMetabIndex)
, mpFirstMetabValue(nullptr)
, mMaxSteps(src.mMaxSteps)
, mMaxStepsReached(false)
, mUseRandomSeed(src.mUseRandomSeed)
, mRandomSeed(src.mRandomSeed)
, mMaxBalance(src.mMaxBalance)
, mMaxIntBeforeStep()
, mReactions()
, mCurrentState()
, mSpeciesRates()
, mRateOffset()
, mUpdateSequences()
, mReactionFlags()
, mFirstReactionFlag(nullptr)
, mMetabFlags()
, mLowerStochLimit(src.mLowerStochLimit)
, mUpperStochLimit(src.mUpperStochLimit)
, mPartitioningInterval(src.mPartitioningInterval)
, mStepsAfterPartitionSystem(0)
, mMetab2React()
, mAmu()
, mAmuOld()
, mpRandomGenerator(nullptr)
, mDG()
, mPQ()
, mAutomaticStepSize(src.mAutomaticStepSize)
{
initializeParameter();
}
Expand Down Expand Up @@ -161,6 +217,9 @@ CTrajectoryMethod::Status CHybridMethod::step(const double & deltaT,
for (i = 0; ((i < mMaxSteps) && (time < endTime)); i++)
{
time = doSingleStep(time, endTime);

if (mAutomaticStepSize)
break;
}

*mpContainerStateTime = time;
Expand Down Expand Up @@ -222,6 +281,8 @@ void CHybridMethod::start()

mMaxStepsReached = false;

mAutomaticStepSize = mpProblem->getAutomaticStepSize();

return;
}

Expand Down Expand Up @@ -255,15 +316,25 @@ C_FLOAT64 CHybridMethod::doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)

if (ds <= endTime) // ds is an absolute time value!
{
bool doNotFire = false;

// if there are deterministic reactions
if (mFirstReactionFlag != NULL) // there is at least one deterministic reaction
{
integrateDeterministicPart(ds - currentTime);
doNotFire = *mpContainerStateTime != ds;
}

mReactions[rIndex].fire();
*mpContainerStateTime = ds;
stateChange(CMath::eStateChange::State);
if (doNotFire)
{
ds = *mpContainerStateTime;
}
else
{
mReactions[rIndex].fire();
*mpContainerStateTime = ds;
stateChange(CMath::eStateChange::State);
}

if (++mStepsAfterPartitionSystem >= mPartitioningInterval)
{
Expand All @@ -281,9 +352,12 @@ C_FLOAT64 CHybridMethod::doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)
if (mFirstReactionFlag != NULL) // there is at least one deterministic reaction
{
integrateDeterministicPart(ds - currentTime);
ds = *mpContainerStateTime;
}
else
{
*mpContainerStateTime = ds;
}

*mpContainerStateTime = ds;

if (++mStepsAfterPartitionSystem >= mPartitioningInterval)
{
Expand All @@ -300,9 +374,12 @@ C_FLOAT64 CHybridMethod::doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)
if (mFirstReactionFlag != NULL) // there is at least one deterministic reaction
{
integrateDeterministicPart(ds - currentTime);
ds = *mpContainerStateTime;
}
else
{
*mpContainerStateTime = ds;
}

*mpContainerStateTime = ds;

if (++mStepsAfterPartitionSystem >= mPartitioningInterval)
{
Expand Down
4 changes: 3 additions & 1 deletion copasi/trajectory/CHybridMethod.h
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// Copyright (C) 2019 by Pedro Mendes, Rector and Visitors of the
// Copyright (C) 2019 - 2024 by Pedro Mendes, Rector and Visitors of the
// University of Virginia, University of Heidelberg, and University
// of Connecticut School of Medicine.
// All rights reserved.
Expand Down Expand Up @@ -462,6 +462,8 @@ class CHybridMethod : public CTrajectoryMethod
*/
CIndexedPriorityQueue mPQ;

bool mAutomaticStepSize;

private:
};

Expand Down
3 changes: 2 additions & 1 deletion copasi/trajectory/CHybridNextReactionLSODAMethod.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// Copyright (C) 2019 - 2020 by Pedro Mendes, Rector and Visitors of the
// Copyright (C) 2019 - 2024 by Pedro Mendes, Rector and Visitors of the
// University of Virginia, University of Heidelberg, and University
// of Connecticut School of Medicine.
// All rights reserved.
Expand Down Expand Up @@ -102,6 +102,7 @@ void CHybridNextReactionLSODAMethod::start()
}
}

mLSODA.setProblem(mpProblem);
mLSODA.setMathContainer(mpContainer);
mLSODA.start();
}
Expand Down
32 changes: 19 additions & 13 deletions copasi/trajectory/CHybridNextReactionRKMethod.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// Copyright (C) 2019 by Pedro Mendes, Rector and Visitors of the
// Copyright (C) 2019 - 2024 by Pedro Mendes, Rector and Visitors of the
// University of Virginia, University of Heidelberg, and University
// of Connecticut School of Medicine.
// All rights reserved.
Expand Down Expand Up @@ -43,6 +43,7 @@
#include "CHybridNextReactionRKMethod.h"
#include "copasi/math/CMathContainer.h"
#include "copasi/math/CMathReaction.h"
#include "copasi/sbml/SBMLImporter.h"

CHybridNextReactionRKMethod::CHybridNextReactionRKMethod(const CDataContainer * pParent,
const CTaskEnum::Method & methodType,
Expand Down Expand Up @@ -89,19 +90,24 @@ void CHybridNextReactionRKMethod::initializeParameter()
*/
void CHybridNextReactionRKMethod::integrateDeterministicPart(C_FLOAT64 dt)
{
C_FLOAT64 integrationTime = 0.0;
size_t stepCounter = 0;
C_FLOAT64 totalStepsTaken = stepCounter * mStepsize;

CHybridStochFlag * react = NULL;

// This method uses a 4th order RungeKutta-method to integrate the deterministic part of the system. Maybe a better numerical method (adaptive stepsize, lsoda, ...) should be introduced here later on

while ((dt - integrationTime) > mStepsize)
while (!SBMLImporter::areApproximatelyEqual(dt, totalStepsTaken, 100.0 * std::numeric_limits< C_FLOAT64 >::epsilon()))
{
rungeKutta(mStepsize); // for the deterministic part of the system
integrationTime += mStepsize;
}
rungeKutta(std::min(dt - totalStepsTaken, mStepsize)); // for the deterministic part of the system
totalStepsTaken = stepCounter * mStepsize + std::min(dt - totalStepsTaken, mStepsize);

rungeKutta(dt - integrationTime);
if (mAutomaticStepSize)
break;

stepCounter += 1;
}

*mpContainerStateTime += totalStepsTaken;
mpContainer->updateSimulatedValues(false);

return;
Expand All @@ -119,7 +125,7 @@ void CHybridNextReactionRKMethod::rungeKutta(C_FLOAT64 dt)
{
size_t i;

CVector< C_FLOAT64 > CurrentState = mCurrentState;
CVector< C_FLOAT64 > PreviousState = mCurrentState;

/* k1 step: k1 = dt*f(x(n)) */
calculateDerivative(temp); // systemState == x(n)
Expand All @@ -132,7 +138,7 @@ void CHybridNextReactionRKMethod::rungeKutta(C_FLOAT64 dt)
/* k2 step: k2 = dt*f(x(n) + k1/2) */
for (i = 0; i < mNumVariableMetabs; i++)
{
temp[i] = k1[i] / 2.0 + CurrentState[i];
temp[i] = k1[i] / 2.0 + PreviousState[i];
}

mCurrentState = temp;
Expand All @@ -146,7 +152,7 @@ void CHybridNextReactionRKMethod::rungeKutta(C_FLOAT64 dt)
/* k3 step: k3 = dt*f(x(n) + k2/2) */
for (i = 0; i < mNumVariableMetabs; i++)
{
temp[i] = k2[i] / 2.0 + CurrentState[i];
temp[i] = k2[i] / 2.0 + PreviousState[i];
}

mCurrentState = temp;
Expand All @@ -160,7 +166,7 @@ void CHybridNextReactionRKMethod::rungeKutta(C_FLOAT64 dt)
/* k4 step: k4 = dt*f(x(n) + k3); */
for (i = 0; i < mNumVariableMetabs; i++)
{
temp[i] = k3[i] + CurrentState[i];
temp[i] = k3[i] + PreviousState[i];
}

mCurrentState = temp;
Expand All @@ -174,7 +180,7 @@ void CHybridNextReactionRKMethod::rungeKutta(C_FLOAT64 dt)
/* Find next position: x(n+1) = x(n) + 1/6*(k1 + 2*k2 + 2*k3 + k4) */
for (i = 0; i < mNumVariableMetabs; i++)
{
temp[i] = CurrentState[i] + (1.0 / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
temp[i] = PreviousState[i] + (1.0 / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
}

mCurrentState = temp;
Expand Down

0 comments on commit 8ea6faa

Please sign in to comment.