COPASI API  4.16.103
Classes | Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | Friends | List of all members
CTrajAdaptiveSA Class Reference

#include <CTrajAdaptiveSA.h>

Inheritance diagram for CTrajAdaptiveSA:
Inheritance graph
[legend]
Collaboration diagram for CTrajAdaptiveSA:
Collaboration graph
[legend]

Classes

class  CReactionDependencies
 

Public Member Functions

 CTrajAdaptiveSA (const CTrajAdaptiveSA &src, const CCopasiContainer *pParent=NULL)
 
virtual bool elevateChildren ()
 
virtual bool isValidProblem (const CCopasiProblem *pProblem)
 
virtual void start (const CState *initialState)
 
virtual Status step (const double &deltaT)
 
 ~CTrajAdaptiveSA ()
 
- Public Member Functions inherited from CTrajectoryMethod
 CTrajectoryMethod (const CTrajectoryMethod &src, const CCopasiContainer *pParent=NULL)
 
const CVector< C_INT > & getRoots () const
 
void setCurrentState (CState *currentState)
 
void setProblem (CTrajectoryProblem *problem)
 
virtual void stateChanged ()
 
 ~CTrajectoryMethod ()
 
- Public Member Functions inherited from CCopasiMethod
 CCopasiMethod (const CCopasiMethod &src, const CCopasiContainer *pParent=NULL)
 
const CCopasiMethod::SubTypegetSubType () const
 
const CCopasiTask::TypegetType () const
 
virtual void load (CReadConfig &configBuffer, CReadConfig::Mode mode=CReadConfig::SEARCH)
 
virtual void print (std::ostream *ostream) const
 
virtual void printResult (std::ostream *ostream) const
 
virtual bool setCallBack (CProcessReport *pCallBack)
 
virtual ~CCopasiMethod ()
 
- Public Member Functions inherited from CCopasiParameterGroup
bool addGroup (const std::string &name)
 
bool addParameter (const CCopasiParameter &parameter)
 
bool addParameter (const std::string &name, const CCopasiParameter::Type type)
 
template<class CType >
bool addParameter (const std::string &name, const CCopasiParameter::Type type, const CType &value)
 
void addParameter (CCopasiParameter *pParameter)
 
CCopasiParameterGroupassertGroup (const std::string &name)
 
template<class CType >
CCopasiParameterassertParameter (const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
 
index_iterator beginIndex () const
 
name_iterator beginName () const
 
 CCopasiParameterGroup (const CCopasiParameterGroup &src, const CCopasiContainer *pParent=NULL)
 
 CCopasiParameterGroup (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &objectType="ParameterGroup")
 
void clear ()
 
index_iterator endIndex () const
 
name_iterator endName () const
 
CCopasiParameterGroupgetGroup (const std::string &name)
 
const CCopasiParameterGroupgetGroup (const std::string &name) const
 
CCopasiParameterGroupgetGroup (const size_t &index)
 
const CCopasiParameterGroupgetGroup (const size_t &index) const
 
size_t getIndex (const std::string &name) const
 
std::string getKey (const std::string &name) const
 
std::string getKey (const size_t &index) const
 
virtual const std::string & getName (const size_t &index) const
 
virtual const CObjectInterfacegetObject (const CCopasiObjectName &cn) const
 
CCopasiParametergetParameter (const std::string &name)
 
const CCopasiParametergetParameter (const std::string &name) const
 
CCopasiParametergetParameter (const size_t &index)
 
const CCopasiParametergetParameter (const size_t &index) const
 
CCopasiParameter::Type getType (const std::string &name) const
 
CCopasiParameter::Type getType (const size_t &index) const
 
std::string getUniqueParameterName (const CCopasiParameter *pParameter) const
 
const CCopasiParameter::ValuegetValue (const std::string &name) const
 
const CCopasiParameter::ValuegetValue (const size_t &index) const
 
CCopasiParameter::ValuegetValue (const std::string &name)
 
CCopasiParameter::ValuegetValue (const size_t &index)
 
CCopasiParameterGroupoperator= (const CCopasiParameterGroup &rhs)
 
bool removeParameter (const std::string &name)
 
bool removeParameter (const size_t &index)
 
template<class CType >
bool setValue (const std::string &name, const CType &value)
 
template<class CType >
bool setValue (const size_t &index, const CType &value)
 
size_t size () const
 
bool swap (const size_t &iFrom, const size_t &iTo)
 
bool swap (index_iterator &from, index_iterator &to)
 
virtual ~CCopasiParameterGroup ()
 
- Public Member Functions inherited from CCopasiParameter
 CCopasiParameter (const CCopasiParameter &src, const CCopasiContainer *pParent=NULL)
 
 CCopasiParameter (const std::string &name, const Type &type, const void *pValue=NULL, const CCopasiContainer *pParent=NULL, const std::string &objectType="Parameter")
 
virtual CCopasiObjectName getCN () const
 
virtual const std::string & getKey () const
 
virtual std::string getObjectDisplayName (bool regular=true, bool richtext=false) const
 
const CCopasiParameter::TypegetType () const
 
const ValuegetValue () const
 
ValuegetValue ()
 
virtual voidgetValuePointer () const
 
CCopasiObjectgetValueReference () const
 
bool isValidValue (const C_FLOAT64 &value) const
 
bool isValidValue (const C_INT32 &value) const
 
bool isValidValue (const unsigned C_INT32 &value) const
 
bool isValidValue (const bool &value) const
 
bool isValidValue (const std::string &value) const
 
bool isValidValue (const CCopasiObjectName &value) const
 
bool isValidValue (const std::vector< CCopasiParameter * > &value) const
 
CCopasiParameteroperator= (const CCopasiParameter &rhs)
 
template<class CType >
bool setValue (const CType &value)
 
bool setValue (const std::vector< CCopasiParameter * > &value)
 
virtual ~CCopasiParameter ()
 
- Public Member Functions inherited from CCopasiContainer
virtual bool add (CCopasiObject *pObject, const bool &adopt=true)
 
 CCopasiContainer (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &type="CN", const unsigned C_INT32 &flag=CCopasiObject::Container)
 
 CCopasiContainer (const CCopasiContainer &src, const CCopasiContainer *pParent=NULL)
 
virtual std::string getChildObjectUnits (const CCopasiObject *pObject) const
 
virtual const objectMapgetObjects () const
 
virtual std::string getUnits () const
 
virtual const CCopasiObjectgetValueObject () const
 
virtual bool remove (CCopasiObject *pObject)
 
virtual ~CCopasiContainer ()
 
- Public Member Functions inherited from CCopasiObject
void addDirectDependency (const CCopasiObject *pObject)
 
 CCopasiObject (const CCopasiObject &src, const CCopasiContainer *pParent=NULL)
 
void clearDirectDependencies ()
 
void clearRefresh ()
 
bool dependsOn (DataObjectSet candidates, const DataObjectSet &context=DataObjectSet()) const
 
void getAllDependencies (DataObjectSet &dependencies, const DataObjectSet &context) const
 
virtual const DataObjectSetgetDirectDependencies (const DataObjectSet &context=DataObjectSet()) const
 
CCopasiContainergetObjectAncestor (const std::string &type) const
 
CCopasiDataModelgetObjectDataModel ()
 
const CCopasiDataModelgetObjectDataModel () const
 
const std::string & getObjectName () const
 
CCopasiContainergetObjectParent () const
 
const std::string & getObjectType () const
 
virtual const
CObjectInterface::ObjectSet
getPrerequisites () const
 
virtual RefreshgetRefresh () const
 
UpdateMethodgetUpdateMethod () const
 
bool hasCircularDependencies (DataObjectSet &candidates, DataObjectSet &verified, const DataObjectSet &context) const
 
bool hasUpdateMethod () const
 
bool isArray () const
 
bool isContainer () const
 
bool isDataModel () const
 
bool isMatrix () const
 
bool isNameVector () const
 
bool isNonUniqueName () const
 
virtual bool isPrerequisiteForContext (const CObjectInterface *pObject, const CMath::SimulationContextFlag &context, const CObjectInterface::ObjectSet &changedObjects) const
 
bool isReference () const
 
bool isRoot () const
 
bool isSeparator () const
 
bool isStaticString () const
 
bool isValueBool () const
 
bool isValueDbl () const
 
bool isValueInt () const
 
bool isValueInt64 () const
 
bool isValueString () const
 
bool isVector () const
 
virtual bool mustBeDeleted (const DataObjectSet &deletedObjects) const
 
void removeDirectDependency (const CCopasiObject *pObject)
 
void setDirectDependencies (const DataObjectSet &directDependencies)
 
bool setObjectName (const std::string &name)
 
virtual bool setObjectParent (const CCopasiContainer *pParent)
 
void setObjectValue (const C_FLOAT64 &value)
 
void setObjectValue (const C_INT32 &value)
 
void setObjectValue (const bool &value)
 
template<class CType >
void setRefresh (CType *pType, void(CType::*method)(void))
 
template<class CType >
void setUpdateMethod (CType *pType, void(CType::*method)(const C_FLOAT64 &))
 
template<class CType >
void setUpdateMethod (CType *pType, void(CType::*method)(const C_INT32 &))
 
template<class CType >
void setUpdateMethod (CType *pType, void(CType::*method)(const bool &))
 
virtual ~CCopasiObject ()
 
- Public Member Functions inherited from CObjectInterface
 CObjectInterface ()
 
virtual ~CObjectInterface ()
 

Static Public Member Functions

static CTrajAdaptiveSAcreateTauLeapMethod ()
 
- Static Public Member Functions inherited from CTrajectoryMethod
static CTrajectoryMethodcreateMethod (CCopasiMethod::SubType subType=CCopasiMethod::deterministic)
 
- Static Public Member Functions inherited from CCopasiObject
static std::vector< Refresh * > buildUpdateSequence (const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())
 
static void setRenameHandler (CRenameHandler *rh)
 

Protected Member Functions

const C_FLOAT64calculateAmu (const size_t &index)
 
 CTrajAdaptiveSA (const CCopasiContainer *pParent=NULL)
 
C_FLOAT64 doSingleSSAStep (const C_FLOAT64 &curTime, const C_FLOAT64 &endTime)
 
C_FLOAT64 doSingleTauLeapStep (const C_FLOAT64 &curTime, const C_FLOAT64 &endTime)
 
- Protected Member Functions inherited from CTrajectoryMethod
 CTrajectoryMethod (const CCopasiMethod::SubType &subType, const CCopasiContainer *pParent=NULL)
 
- Protected Member Functions inherited from CCopasiMethod
 CCopasiMethod (const CCopasiTask::Type &taskType, const SubType &subType, const CCopasiContainer *pParent=NULL)
 
- Protected Member Functions inherited from CCopasiParameterGroup
 CCopasiParameterGroup ()
 
- Protected Member Functions inherited from CCopasiContainer
template<class CType >
CCopasiObjectaddMatrixReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
 
template<class CType >
CCopasiObjectaddObjectReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
 
template<class CType >
CCopasiObjectaddVectorReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
 
void initObjects ()
 
- Protected Member Functions inherited from CCopasiObject
 CCopasiObject ()
 
 CCopasiObject (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &type="CN", const unsigned C_INT32 &flag=0)
 

Protected Attributes

C_FLOAT64 mA0
 
CVector< C_FLOAT64mAmu
 
bool mDoCorrection
 
unsigned C_INT32 mMaxSteps
 
bool mMaxStepsReached
 
CState mMethodState
 
size_t mNextReactionIndex
 
C_FLOAT64 mNextReactionTime
 
size_t mNumReactions
 
size_t mNumReactionSpecies
 
size_t mNumSpecies
 
CVector< const C_FLOAT64 * > mPartitionedAmu
 
CVector< const
CReactionDependencies * > 
mPartitionedDependencies
 
CModelmpModel
 
CRandommpRandomGenerator
 
std::vector
< CReactionDependencies
mReactionDependencies
 
- Protected Attributes inherited from CTrajectoryMethod
CStatempCurrentState
 
CTrajectoryProblemmpProblem
 
CVector< C_INTmRoots
 
- Protected Attributes inherited from CCopasiMethod
CProcessReportmpCallBack
 
- Protected Attributes inherited from CCopasiParameter
std::string mKey
 
CCopasiObjectmpValueReference
 
size_t mSize
 
Value mValue
 
- Protected Attributes inherited from CCopasiContainer
objectMap mObjects
 

Private Member Functions

void initializeParameter ()
 

Private Attributes

CVector< C_FLOAT64mAvgDX
 
C_FLOAT64 mEpsilon
 
size_t mFirstReactionSpeciesIndex
 
CVector< size_t > mMaxReactionFiring
 
CVector< C_FLOAT64 * > mPartitionedReactionFiring
 
C_FLOAT64mpMethodSpecies
 
CVector< C_FLOAT64mReactionFiring
 
CVector< C_FLOAT64mSigDX
 
CVector< C_FLOAT64mSpeciesAfterTau
 
C_FLOAT64 mSSAStepCounter
 
std::vector< Refresh * > mTauCalculations
 

Friends

CTrajectoryMethodCTrajectoryMethod::createMethod (CCopasiMethod::SubType subType)
 

Additional Inherited Members

- Public Types inherited from CTrajectoryMethod
enum  Status { FAILURE = -1, NORMAL = 0, ROOT = 1 }
 
- Public Types inherited from CCopasiMethod
enum  SubType {
  unset = 0, RandomSearch, RandomSearchMaster, SimulatedAnnealing,
  CoranaWalk, DifferentialEvolution, ScatterSearch, GeneticAlgorithm,
  EvolutionaryProgram, SteepestDescent, HybridGASA, GeneticAlgorithmSR,
  HookeJeeves, LevenbergMarquardt, NelderMead, SRES,
  Statistics, ParticleSwarm, Praxis, TruncatedNewton,
  Newton, deterministic, LSODAR, directMethod,
  stochastic, tauLeap, adaptiveSA, hybrid,
  hybridLSODA, hybridODE45, DsaLsodar, tssILDM,
  tssILDMModified, tssCSP, mcaMethodReder, scanMethod,
  lyapWolf, sensMethod, EFMAlgorithm, EFMBitPatternTreeAlgorithm,
  EFMBitPatternAlgorithm, Householder, crossSectionMethod, linearNoiseApproximation
}
 
- Public Types inherited from CCopasiParameterGroup
typedef parameterGroup::iterator index_iterator
 
typedef
CCopasiContainer::objectMap::iterator 
name_iterator
 
typedef std::vector
< CCopasiParameter * > 
parameterGroup
 
- Public Types inherited from CCopasiParameter
enum  Type {
  DOUBLE = 0, UDOUBLE, INT, UINT,
  BOOL, GROUP, STRING, CN,
  KEY, FILE, EXPRESSION, INVALID
}
 
- Public Types inherited from CCopasiContainer
typedef std::multimap
< std::string, CCopasiObject * > 
objectMap
 
- Public Types inherited from CCopasiObject
typedef std::set< const
CCopasiObject * > 
DataObjectSet
 
typedef std::vector< Refresh * > DataUpdateSequence
 
- Public Types inherited from CObjectInterface
typedef std::set< const
CObjectInterface * > 
ObjectSet
 
typedef std::vector
< CObjectInterface * > 
UpdateSequence
 
- Static Public Attributes inherited from CCopasiMethod
static const std::string SubTypeName []
 
static const char * XMLSubType []
 
- Static Public Attributes inherited from CCopasiParameter
static const std::string TypeName []
 
static const char * XMLType []
 
- Static Public Attributes inherited from CCopasiContainer
static const std::vector
< CCopasiContainer * > 
EmptyList
 
- Protected Types inherited from CCopasiObject
enum  Flag {
  Container = 0x1, Vector = 0x2, Matrix = 0x4, NameVector = 0x8,
  Reference = 0x10, ValueBool = 0x20, ValueInt = 0x40, ValueInt64 = 0x80,
  ValueDbl = 0x100, NonUniqueName = 0x200, StaticString = 0x400, ValueString = 0x800,
  Separator = 0x1000, ModelEntity = 0x2000, Array = 0x4000, DataModel = 0x8000,
  Root = 0x10000, Gui = 0x20000
}
 
- Static Protected Attributes inherited from CCopasiObject
static CRenameHandlersmpRenameHandler = NULL
 

Detailed Description

Definition at line 36 of file CTrajAdaptiveSA.h.

Constructor & Destructor Documentation

CTrajAdaptiveSA::CTrajAdaptiveSA ( const CCopasiContainer pParent = NULL)
protected

Default constructor.

Parameters
constCCopasiContainer * pParent (default: NULL)

Definition at line 86 of file CTrajAdaptiveSA.cpp.

References initializeParameter().

Referenced by createTauLeapMethod().

86  :
89  mReactionFiring(0),
91  mAvgDX(0),
92  mSigDX(0),
93  mpMethodSpecies(NULL),
96  mpModel(NULL),
97  mNumReactions(0),
98  mNumSpecies(0),
99  mMaxSteps(1000000),
100  mNextReactionTime(0.0),
102  mDoCorrection(true),
103  mAmu(0),
104  mPartitionedAmu(0),
105  mMethodState(),
108  mA0(0.0),
109  mMaxStepsReached(false)
110 {
112 }
CVector< size_t > mMaxReactionFiring
CRandom * mpRandomGenerator
#define C_INVALID_INDEX
Definition: copasi.h:222
CVector< C_FLOAT64 > mAvgDX
CVector< const C_FLOAT64 * > mPartitionedAmu
CVector< C_FLOAT64 > mReactionFiring
CVector< C_FLOAT64 > mSpeciesAfterTau
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
std::vector< CReactionDependencies > mReactionDependencies
unsigned C_INT32 mMaxSteps
CVector< C_FLOAT64 > mAmu
CVector< const CReactionDependencies * > mPartitionedDependencies
C_FLOAT64 mNextReactionTime
CVector< C_FLOAT64 > mSigDX
CVector< C_FLOAT64 * > mPartitionedReactionFiring
C_FLOAT64 * mpMethodSpecies
CTrajAdaptiveSA::CTrajAdaptiveSA ( const CTrajAdaptiveSA src,
const CCopasiContainer pParent = NULL 
)

Copy constructor.

Parameters
constCTrajAdaptiveSA & src,
constCCopasiContainer * pParent (Default: NULL)

Definition at line 114 of file CTrajAdaptiveSA.cpp.

References initializeParameter().

115  :
116  CTrajectoryMethod(src, pParent),
120  mAvgDX(src.mAvgDX),
121  mSigDX(src.mSigDX),
125  mpModel(src.mpModel),
128  mMaxSteps(src.mMaxSteps),
132  mAmu(src.mAmu),
137  mA0(src.mA0),
139 {
141 }
CVector< size_t > mMaxReactionFiring
CRandom * mpRandomGenerator
CVector< C_FLOAT64 > mAvgDX
CVector< const C_FLOAT64 * > mPartitionedAmu
CVector< C_FLOAT64 > mReactionFiring
CVector< C_FLOAT64 > mSpeciesAfterTau
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
std::vector< CReactionDependencies > mReactionDependencies
unsigned C_INT32 mMaxSteps
CVector< C_FLOAT64 > mAmu
CVector< const CReactionDependencies * > mPartitionedDependencies
C_FLOAT64 mNextReactionTime
CVector< C_FLOAT64 > mSigDX
CVector< C_FLOAT64 * > mPartitionedReactionFiring
C_FLOAT64 * mpMethodSpecies
CTrajAdaptiveSA::~CTrajAdaptiveSA ( )

Destructor.

Definition at line 143 of file CTrajAdaptiveSA.cpp.

References mpRandomGenerator, and pdelete.

144 {
146 }
#define pdelete(p)
Definition: copasi.h:215
CRandom * mpRandomGenerator

Member Function Documentation

const C_FLOAT64 & CTrajAdaptiveSA::calculateAmu ( const size_t &  index)
protected

Calculate the propensity of the indexed reaction

Parameters
constsize_t & index

Definition at line 844 of file CTrajAdaptiveSA.cpp.

References CVectorCore< CType >::array(), C_FLOAT64, mAmu, mDoCorrection, CTrajAdaptiveSA::CReactionDependencies::mMethodSubstrates, CTrajAdaptiveSA::CReactionDependencies::mModelSubstrates, CTrajAdaptiveSA::CReactionDependencies::mpParticleFlux, mReactionDependencies, CTrajAdaptiveSA::CReactionDependencies::mSubstrateMultiplier, and CVectorCore< CType >::size().

Referenced by doSingleSSAStep(), doSingleTauLeapStep(), and start().

845 {
846  const CReactionDependencies & Dependencies = mReactionDependencies[index];
847  C_FLOAT64 & Amu = mAmu[index];
848 
849  Amu = *Dependencies.mpParticleFlux;
850 
851  if (Amu < 0.0)
852  {
853  // TODO CRITICAL Create a warning message
854  Amu = 0.0;
855  }
856 
857  if (!mDoCorrection)
858  {
859  return Amu;
860  }
861 
862  C_FLOAT64 SubstrateMultiplier = 1.0;
863  C_FLOAT64 SubstrateDevisor = 1.0;
864  C_FLOAT64 Multiplicity;
865  C_FLOAT64 LowerBound;
866  C_FLOAT64 Number;
867 
868  bool ApplyCorrection = false;
869 
870  const C_FLOAT64 * pMultiplicity = Dependencies.mSubstrateMultiplier.array();
871  const C_FLOAT64 * pEndMultiplicity = pMultiplicity + Dependencies.mSubstrateMultiplier.size();
872  C_FLOAT64 *const* ppLocalSubstrate = Dependencies.mMethodSubstrates.array();
873  C_FLOAT64 *const* ppModelSubstrate = Dependencies.mModelSubstrates.array();
874 
875  for (; pMultiplicity != pEndMultiplicity; ++pMultiplicity, ++ppLocalSubstrate, ++ppModelSubstrate)
876  {
877  Multiplicity = *pMultiplicity;
878 
879  // TODO We should check the error introduced through rounding.
880  **ppLocalSubstrate = floor(**ppModelSubstrate + 0.5);
881 
882  if (Multiplicity > 1.01)
883  {
884  ApplyCorrection = true;
885 
886  Number = **ppLocalSubstrate;
887 
888  LowerBound = Number - Multiplicity;
889  SubstrateDevisor *= pow(Number, Multiplicity - 1.0); //optimization
890  Number -= 1.0;
891 
892  while (Number > LowerBound)
893  {
894  SubstrateMultiplier *= Number;
895  Number -= 1.0;
896  }
897  }
898  }
899 
900  // at least one substrate particle number is zero
901  if (SubstrateMultiplier < 0.5 || SubstrateDevisor < 0.5)
902  {
903  Amu = 0.0;
904  }
905  else if (ApplyCorrection)
906  {
907  Amu *= SubstrateMultiplier / SubstrateDevisor;
908  }
909 
910  return Amu;
911 }
std::vector< CReactionDependencies > mReactionDependencies
CVector< C_FLOAT64 > mAmu
#define C_FLOAT64
Definition: copasi.h:92
CTrajAdaptiveSA * CTrajAdaptiveSA::createTauLeapMethod ( )
static

Chooses a stochastic method adequate for the problem

Definition at line 77 of file CTrajAdaptiveSA.cpp.

References CTrajAdaptiveSA().

78 {
79  CTrajAdaptiveSA * pMethod = NULL;
80 
81  pMethod = new CTrajAdaptiveSA();
82 
83  return pMethod;
84 }
CTrajAdaptiveSA(const CCopasiContainer *pParent=NULL)
C_FLOAT64 CTrajAdaptiveSA::doSingleSSAStep ( const C_FLOAT64 curTime,
const C_FLOAT64 endTime 
)
protected

Fire the next reaction if it fire before the endTime

Parameters
constC_FLOAT64 & curTime
constC_FLOAT64 & endTime
Returns
C_FLOAT64 timeAfterStep

Definition at line 753 of file CTrajAdaptiveSA.cpp.

References CVectorCore< CType >::array(), C_FLOAT64, C_INVALID_INDEX, calculateAmu(), CCopasiMessage::EXCEPTION, CRandom::getRandomOO(), mA0, mAmu, CTrajAdaptiveSA::CReactionDependencies::mCalculations, MCTrajectoryMethod, CTrajAdaptiveSA::CReactionDependencies::mDependentReactions, CTrajAdaptiveSA::CReactionDependencies::mMethodSpecies, CTrajAdaptiveSA::CReactionDependencies::mModelSpecies, mNextReactionIndex, mNextReactionTime, mNumReactions, mpModel, mpRandomGenerator, mReactionDependencies, CTrajAdaptiveSA::CReactionDependencies::mSpeciesMultiplier, CModel::setTime(), and CVectorCore< CType >::size().

Referenced by step().

754 {
756  {
757  if (mA0 == 0)
758  {
759  return endTime - curTime;
760  }
761 
762  // We need to throw an exception if mA0 is NaN
763  if (isnan(mA0))
764  {
766  }
767 
768  mNextReactionTime = curTime - log(mpRandomGenerator->getRandomOO()) / mA0;
769  }
770 
771  if (mNextReactionTime >= endTime)
772  {
773  return endTime - curTime;
774  }
775 
776  // Update the model time (for explicitly time dependent models)
778 
779  // We are sure that we have at least 1 reaction
780  mNextReactionIndex = 0;
781 
782  C_FLOAT64 sum = 0.0;
784 
785  C_FLOAT64 * pAmu = mAmu.array();
786  C_FLOAT64 * endAmu = pAmu + mNumReactions;
787 
788  for (; (sum < rand) && (pAmu != endAmu); ++pAmu, ++mNextReactionIndex)
789  {
790  sum += *pAmu;
791  }
792 
794 
795  CReactionDependencies & Dependencies = mReactionDependencies[mNextReactionIndex];
796 
797  // Update the method internal and model species numbers
798  C_FLOAT64 ** ppModelSpecies = Dependencies.mModelSpecies.array();
799  C_FLOAT64 ** ppLocalSpecies = Dependencies.mMethodSpecies.array();
800  C_FLOAT64 * pMultiplicity = Dependencies.mSpeciesMultiplier.array();
801  C_FLOAT64 * pMultiplicityEnd = pMultiplicity + Dependencies.mSpeciesMultiplier.size();
802  //if(mAmu.array()[mNextReactionIndex]==0) //Important check
803  // printf("Error 0\n");
804 
805  for (; pMultiplicity != pMultiplicityEnd; ++pMultiplicity, ++ppLocalSpecies, ++ppModelSpecies)
806  {
807  **ppLocalSpecies += *pMultiplicity;
808  **ppModelSpecies = **ppLocalSpecies;
809  }
810 
811  // Calculate all values which depend on the firing reaction
812  std::vector< Refresh * >::const_iterator itCalcualtion = Dependencies.mCalculations.begin();
813  std::vector< Refresh * >::const_iterator endCalcualtion = Dependencies.mCalculations.end();
814 
815  while (itCalcualtion != endCalcualtion)
816  {
817  (**itCalcualtion++)();
818  }
819 
820  // calculate the propensities which depend on the firing reaction
821  size_t * pDependentReaction = Dependencies.mDependentReactions.array();
822  size_t * endDependentReactions = pDependentReaction + Dependencies.mDependentReactions.size();
823 
824  for (; pDependentReaction != endDependentReactions; ++pDependentReaction)
825  {
826  calculateAmu(*pDependentReaction);
827  }
828 
829  // calculate the total propensity
830  pAmu = mAmu.array();
831 
832  mA0 = 0.0;
833 
834  for (; pAmu != endAmu; ++pAmu)
835  {
836  mA0 += *pAmu;
837  }
838 
840 
841  return mNextReactionTime - curTime;
842 }
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
CRandom * mpRandomGenerator
#define C_INVALID_INDEX
Definition: copasi.h:222
const C_FLOAT64 & calculateAmu(const size_t &index)
#define MCTrajectoryMethod
std::vector< CReactionDependencies > mReactionDependencies
void setTime(const C_FLOAT64 &time)
Definition: CModel.cpp:1187
CVector< C_FLOAT64 > mAmu
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
C_FLOAT64 mNextReactionTime
C_FLOAT64 CTrajAdaptiveSA::doSingleTauLeapStep ( const C_FLOAT64 curTime,
const C_FLOAT64 endTime 
)
protected

Fire the next reaction if it fire before the endTime

Parameters
constC_FLOAT64 & curTime
constC_FLOAT64 & endTime
Returns
C_FLOAT64 timeAfterStep

Definition at line 482 of file CTrajAdaptiveSA.cpp.

References CVectorCore< CType >::array(), C_FLOAT64, C_INVALID_INDEX, calculateAmu(), CCopasiMessage::EXCEPTION, CRandom::getRandomOO(), CRandom::getRandomPoisson(), mA0, mAmu, mAvgDX, max, MCTrajectoryMethod, mEpsilon, mMaxReactionFiring, CTrajAdaptiveSA::CReactionDependencies::mMethodSpeciesIndex, mMethodState, mNumReactions, mNumReactionSpecies, mPartitionedAmu, mPartitionedDependencies, mPartitionedReactionFiring, mpMethodSpecies, mpModel, mpRandomGenerator, mReactionDependencies, mReactionFiring, mSigDX, mSpeciesAfterTau, CTrajAdaptiveSA::CReactionDependencies::mSpeciesMultiplier, mSSAStepCounter, CModel::setState(), CState::setTime(), CVectorCore< CType >::size(), SSA_MULTIPLE, SSA_UPPER_NUM, CModel::updateSimulatedValues(), and UPPER_LIMIT.

Referenced by step().

483 {
484  std::vector< CReactionDependencies >::const_iterator itDependencies = mReactionDependencies.begin();
485  std::vector< CReactionDependencies >::const_iterator endDependencies = mReactionDependencies.end();
486  size_t * pMaxReactionFiring = mMaxReactionFiring.array();
487 
488  for (; itDependencies != endDependencies; ++itDependencies, ++pMaxReactionFiring)
489  {
490  C_FLOAT64 *const* ppModelSpecies = itDependencies->mModelSpecies.array();
491  const C_FLOAT64 * pMultiplicity = itDependencies->mSpeciesMultiplier.array();
492  const C_FLOAT64 * pMultiplicityEnd = pMultiplicity + itDependencies->mSpeciesMultiplier.size();
493 
494  *pMaxReactionFiring = std::numeric_limits< size_t >::max(); // Assigned a maximum value
495 
496  for (; pMultiplicity != pMultiplicityEnd; pMultiplicity++, ppModelSpecies++)
497  {
498  if (*pMultiplicity < 0)
499  {
500  size_t TmpMax;
501 
502  if ((TmpMax = (size_t) fabs(**ppModelSpecies / *pMultiplicity)) < *pMaxReactionFiring)
503  {
504  *pMaxReactionFiring = TmpMax;
505  }
506  }
507  }
508  }
509 
510  size_t NonCriticalReactions = 0;
511  size_t InsertCritical = mNumReactions - 1;
512 
513  pMaxReactionFiring = mMaxReactionFiring.array();
514  const C_FLOAT64 * pAmu = mAmu.array();
515 
516  itDependencies = mReactionDependencies.begin();
517  C_FLOAT64 * pReactionFiring = mReactionFiring.array();
518 
519  for (; itDependencies != endDependencies; ++itDependencies, ++pAmu, ++pMaxReactionFiring, ++pReactionFiring)
520  {
521  *pReactionFiring = 0;
522 
523  if (*pAmu == 0 ||
524  *pMaxReactionFiring > UPPER_LIMIT)
525  {
526  mPartitionedDependencies[NonCriticalReactions] = &(*itDependencies);
527  mPartitionedAmu[NonCriticalReactions] = pAmu;
528  mPartitionedReactionFiring[NonCriticalReactions] = pReactionFiring;
529  NonCriticalReactions++;
530  }
531  else
532  {
533  mPartitionedDependencies[InsertCritical] = &(*itDependencies);
534  mPartitionedAmu[InsertCritical] = pAmu;
535  mPartitionedReactionFiring[InsertCritical] = pReactionFiring;
536  InsertCritical--;
537  }
538  }
539 
540  mAvgDX = 0.0;
541  mSigDX = 0.0;
542 
543  const CReactionDependencies **ppOrderedReactions = mPartitionedDependencies.array();
544  const C_FLOAT64 ** ppOrderedAmu = mPartitionedAmu.array();
545  const C_FLOAT64 ** ppOrderedAmuEnd = ppOrderedAmu + NonCriticalReactions;
546 
547  for (; ppOrderedAmu != ppOrderedAmuEnd; ++ppOrderedReactions, ++ppOrderedAmu)
548  {
549  const C_FLOAT64 * pMultiplicity = (*ppOrderedReactions)->mSpeciesMultiplier.array();
550  const C_FLOAT64 * pMultiplicityEnd = pMultiplicity + (*ppOrderedReactions)->mSpeciesMultiplier.size();
551  const size_t * pIndex = (*ppOrderedReactions)->mMethodSpeciesIndex.array();
552 
553  for (; pMultiplicity != pMultiplicityEnd; pMultiplicity++, pIndex++)
554  {
555  mAvgDX[*pIndex] += **ppOrderedAmu * *pMultiplicity;
556  mSigDX[*pIndex] += **ppOrderedAmu * *pMultiplicity * *pMultiplicity;
557  }
558  }
559 
560  C_FLOAT64 MaxTau;
561 
562  if (NonCriticalReactions == 0)
563  {
564  MaxTau = 0;
565  }
566  else
567  {
568  C_FLOAT64 ex, t1, t2;
569  t1 = t2 = std::numeric_limits< C_FLOAT64 >::infinity();
570 
571  for (size_t i = 0; i < mNumReactionSpecies; i++)
572  {
573  C_FLOAT64 t3, t4, t5, t6;
574 
575  ex = (*(mpMethodSpecies + i) * mEpsilon) + 1.0;
576  t3 = fabs(mAvgDX[i]);
577  t4 = mSigDX[i];
578  t5 = ex / t3;
579  t6 = ex * ex / t4;
580 
581  if (t3 != 0 && t5 < t1) t1 = t5;
582 
583  if (t4 != 0 && t6 < t2) t2 = t6;
584  }
585 
586  if (t1 < t2)
587  MaxTau = t1;
588  else
589  MaxTau = t2;
590  }
591 
592  if (MaxTau < (SSA_MULTIPLE / mA0) ||
593  MaxTau == std::numeric_limits< C_FLOAT64 >::infinity())
594  {
598 
599  return 0;
600  }
601 
602  C_FLOAT64 AmuCritical = 0;
603  ppOrderedAmu = mPartitionedAmu.array() + NonCriticalReactions;
604  ppOrderedAmuEnd = mPartitionedAmu.array() + mNumReactions;
605 
606  for (; ppOrderedAmu != ppOrderedAmuEnd; ++ppOrderedAmu)
607  {
608  AmuCritical += **ppOrderedAmu;
609  }
610 
611  C_FLOAT64 CriticalReactionTau;
612 
613  if (NonCriticalReactions == mNumReactions || AmuCritical == 0)
614  CriticalReactionTau = std::numeric_limits< C_FLOAT64 >::infinity();
615  else
616  CriticalReactionTau = -log(mpRandomGenerator->getRandomOO()) / AmuCritical;
617 
618  bool isUpdated = false;
619 
620  C_FLOAT64 Tau;
621 
622  while (!isUpdated)
623  {
624  Tau = MaxTau;
625 
626  if (CriticalReactionTau < Tau || Tau == 0) Tau = CriticalReactionTau;
627 
628  if ((endTime - curTime) < Tau) Tau = (endTime - curTime);
629 
630  // Calculate the firing of non critical reactions
631  ppOrderedAmu = mPartitionedAmu.array();
632  ppOrderedAmuEnd = mPartitionedAmu.array() + NonCriticalReactions;
633 
634  C_FLOAT64 **ppOrderedReactionFiring = mPartitionedReactionFiring.array();
635 
636  for (; ppOrderedAmu != ppOrderedAmuEnd; ++ppOrderedAmu, ++ppOrderedReactionFiring)
637  {
638  C_FLOAT64 Lambda = **ppOrderedAmu * Tau;
639 
640  if (Lambda < 0.0)
642  else if (Lambda > 2.0e9)
644 
645  **ppOrderedReactionFiring = mpRandomGenerator->getRandomPoisson(Lambda);
646  }
647 
648  size_t CriticalReactionIndex = C_INVALID_INDEX;
649 
650  if (CriticalReactionTau <= Tau)
651  {
652  // Determine the critical reaction which fires.
653  C_FLOAT64 sum = 0;
654  C_FLOAT64 rand = mpRandomGenerator->getRandomOO() * AmuCritical;
655 
656  ppOrderedAmu = mPartitionedAmu.array() + NonCriticalReactions;
657  ppOrderedAmuEnd = mPartitionedAmu.array() + mNumReactions;
658 
659  CriticalReactionIndex = NonCriticalReactions;
660 
661  for (; (sum < rand) && (ppOrderedAmu != ppOrderedAmuEnd); ++ppOrderedAmu, ++CriticalReactionIndex)
662  {
663  sum += **ppOrderedAmu;
664  }
665 
666  CriticalReactionIndex--;
667  }
668 
669  C_FLOAT64 * pSpeciesAfterTau = mSpeciesAfterTau.array();
670  C_FLOAT64 * pSpeciesAfterTauEnd = pSpeciesAfterTau + mNumReactionSpecies;
671  C_FLOAT64 * pMethodSpecies = mpMethodSpecies;
672 
673  for (; pSpeciesAfterTau != pSpeciesAfterTauEnd; ++pSpeciesAfterTau, ++pMethodSpecies)
674  {
675  *pSpeciesAfterTau = *pMethodSpecies;
676  }
677 
678  ppOrderedReactions = mPartitionedDependencies.array();
679  const CReactionDependencies **ppOrderedReactionsEnd = ppOrderedReactions + NonCriticalReactions;
680  ppOrderedReactionFiring = mPartitionedReactionFiring.array();
681 
682  for (; ppOrderedReactions != ppOrderedReactionsEnd; ++ppOrderedReactions, ++ppOrderedReactionFiring)
683  {
684  const C_FLOAT64 * pMultiplicity = (*ppOrderedReactions)->mSpeciesMultiplier.array();
685  const C_FLOAT64 * pMultiplicityEnd = pMultiplicity + (*ppOrderedReactions)->mSpeciesMultiplier.size();
686  const size_t *pIndex = (*ppOrderedReactions)->mMethodSpeciesIndex.array();
687 
688  for (; pMultiplicity != pMultiplicityEnd; pMultiplicity++, pIndex++)
689  {
690  mSpeciesAfterTau[*pIndex] += (*pMultiplicity * **ppOrderedReactionFiring);
691  }
692 
693  **ppOrderedReactionFiring = 0.0;
694  }
695 
696  if (CriticalReactionIndex != C_INVALID_INDEX)
697  {
698  const CReactionDependencies * pCriticalReaction = mPartitionedDependencies[CriticalReactionIndex];
699 
700  const C_FLOAT64 * pMultiplicity = pCriticalReaction->mSpeciesMultiplier.array();
701  const C_FLOAT64 * pMultiplicityEnd = pMultiplicity + pCriticalReaction->mSpeciesMultiplier.size();
702  const size_t *pIndex = pCriticalReaction->mMethodSpeciesIndex.array();
703 
704  for (; pMultiplicity != pMultiplicityEnd; pMultiplicity++, pIndex++)
705  {
706  mSpeciesAfterTau[*pIndex] += *pMultiplicity;
707  }
708  }
709 
710  pSpeciesAfterTau = mSpeciesAfterTau.array();
711 
712  for (; pSpeciesAfterTau != pSpeciesAfterTauEnd; ++pSpeciesAfterTau)
713  {
714  if (*pSpeciesAfterTau < 0)
715  break;
716  }
717 
718  if (pSpeciesAfterTau != pSpeciesAfterTauEnd)
719  {
720  isUpdated = false;
721  MaxTau = MaxTau / 2.0;
722  }
723  else
724  {
725  isUpdated = true;
726 
727  pSpeciesAfterTau = mSpeciesAfterTau.array();
728  pMethodSpecies = mpMethodSpecies;
729 
730  for (; pSpeciesAfterTau != pSpeciesAfterTauEnd; ++pSpeciesAfterTau, ++pMethodSpecies)
731  {
732  *pMethodSpecies = *pSpeciesAfterTau;
733  }
734  }
735  }
736 
737  // Update the model time (for explicitly time dependent models)
738  mMethodState.setTime(curTime + Tau);
741 
742  mA0 = 0;
743  size_t i;
744 
745  for (i = 0; i < mNumReactions; i++)
746  {
747  mA0 += calculateAmu(i);
748  }
749 
750  return Tau;
751 }
#define UPPER_LIMIT
#define SSA_MULTIPLE
CVector< size_t > mMaxReactionFiring
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
size_t mNumReactionSpecies
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
CRandom * mpRandomGenerator
#define C_INVALID_INDEX
Definition: copasi.h:222
CVector< C_FLOAT64 > mAvgDX
CVector< const C_FLOAT64 * > mPartitionedAmu
const C_FLOAT64 & calculateAmu(const size_t &index)
CVector< C_FLOAT64 > mReactionFiring
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
C_FLOAT64 mSSAStepCounter
#define MCTrajectoryMethod
CVector< C_FLOAT64 > mSpeciesAfterTau
virtual C_FLOAT64 getRandomPoisson(const C_FLOAT64 &mean)
Definition: CRandom.cpp:314
std::vector< CReactionDependencies > mReactionDependencies
void setState(const CState &state)
Definition: CModel.cpp:1785
CVector< C_FLOAT64 > mAmu
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
#define SSA_UPPER_NUM
CVector< const CReactionDependencies * > mPartitionedDependencies
CVector< C_FLOAT64 > mSigDX
CVector< C_FLOAT64 * > mPartitionedReactionFiring
C_FLOAT64 * mpMethodSpecies
#define max(a, b)
Definition: f2c.h:176
bool CTrajAdaptiveSA::elevateChildren ( )
virtual

This methods must be called to elevate subgroups to derived objects. The default implementation does nothing.

Returns
bool success

Reimplemented from CCopasiParameterGroup.

Definition at line 156 of file CTrajAdaptiveSA.cpp.

References initializeParameter().

157 {
159  return true;
160 }
void CTrajAdaptiveSA::initializeParameter ( )
private

Initialize the method parameter

Definition at line 148 of file CTrajAdaptiveSA.cpp.

References CCopasiParameterGroup::assertParameter(), CCopasiParameter::BOOL, C_FLOAT64, C_INT32, CCopasiParameter::DOUBLE, EPS, CCopasiParameter::INT, and CCopasiParameter::UINT.

Referenced by CTrajAdaptiveSA(), and elevateChildren().

149 {
151  assertParameter("Max Internal Steps", CCopasiParameter::INT, (C_INT32) 1000000);
152  assertParameter("Use Random Seed", CCopasiParameter::BOOL, false);
153  assertParameter("Random Seed", CCopasiParameter::UINT, (unsigned C_INT32) 1);
154 }
#define C_INT32
Definition: copasi.h:90
#define C_FLOAT64
Definition: copasi.h:92
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
#define EPS
bool CTrajAdaptiveSA::isValidProblem ( const CCopasiProblem pProblem)
virtual

Check if the method is suitable for this problem

Returns
bool suitability of the method

Reimplemented from CTrajectoryMethod.

Definition at line 411 of file CTrajAdaptiveSA.cpp.

References CStateTemplate::beginIndependent(), CStateTemplate::endIndependent(), CCopasiMessage::EXCEPTION, CTrajectoryProblem::getDuration(), CCopasiProblem::getModel(), CModel::getStateTemplate(), CModel::getTotSteps(), CCopasiParameter::getValue(), CTrajectoryMethod::isValidProblem(), MCTrajectoryMethod, CModelEntity::ODE, CCopasiParameter::Value::pINT, and CModel::suitableForStochasticSimulation().

412 {
413  if (!CTrajectoryMethod::isValidProblem(pProblem)) return false;
414 
415  const CTrajectoryProblem * pTP = dynamic_cast<const CTrajectoryProblem *>(pProblem);
416 
417  if (pTP->getDuration() < 0.0)
418  {
419  //back integration not possible
421  return false;
422  }
423 
424  if (pTP->getModel()->getTotSteps() < 1)
425  {
426  //at least one reaction necessary
428  return false;
429  }
430 
431  // check for ODEs
432  const CStateTemplate & StateTemplate = pTP->getModel()->getStateTemplate();
433  CModelEntity *const* ppEntity = StateTemplate.beginIndependent();
434  CModelEntity *const* ppEntityEnd = StateTemplate.endIndependent();
435 
436  for (; ppEntity != ppEntityEnd; ++ppEntity)
437  {
438  if ((*ppEntity)->getStatus() == CModelEntity::ODE)
439  {
440  if (dynamic_cast<const CModelValue *>(*ppEntity) != NULL)
441  {
442  // global quantity ode rule found
444  return false;
445  }
446  else if (dynamic_cast<const CCompartment *>(*ppEntity) != NULL)
447  {
448  // compartment ode rule found
450  return false;
451  }
452  else
453  {
454  // species ode rule found
456  return false;
457  }
458  }
459  }
460 
461  //TODO: rewrite CModel::suitableForStochasticSimulation() to use
462  // CCopasiMessage
463  std::string message = pTP->getModel()->suitableForStochasticSimulation();
464 
465  if (message != "")
466  {
467  //model not suitable, message describes the problem
468  CCopasiMessage(CCopasiMessage::EXCEPTION, message.c_str());
469  return false;
470  }
471 
472  if (* getValue("Max Internal Steps").pINT <= 0)
473  {
474  //max steps should be at least 1
476  return false;
477  }
478 
479  return true;
480 }
virtual bool isValidProblem(const CCopasiProblem *pProblem)
size_t getTotSteps() const
Definition: CModel.cpp:1136
#define MCTrajectoryMethod
const C_FLOAT64 & getDuration() const
const Value & getValue() const
CModelEntity ** beginIndependent()
Definition: CState.cpp:208
const CStateTemplate & getStateTemplate() const
Definition: CModel.cpp:1172
std::string suitableForStochasticSimulation() const
Definition: CModel.cpp:3901
CModel * getModel() const
CModelEntity ** endIndependent()
Definition: CState.cpp:209
void CTrajAdaptiveSA::start ( const CState initialState)
virtual

This instructs the method to prepare for integration starting with the initialState given.

Parameters
const CState *initialState

Reimplemented from CTrajectoryMethod.

Definition at line 194 of file CTrajAdaptiveSA.cpp.

References CVectorCore< CType >::array(), CCopasiVector< T >::begin(), CStateTemplate::beginIndependent(), CState::beginIndependent(), CCopasiObject::buildUpdateSequence(), C_FLOAT64, C_INT32, C_INVALID_INDEX, calculateAmu(), CModel::deterministic, CCopasiVector< T >::end(), CStateTemplate::endFixed(), CStateTemplate::getIndex(), CModel::getMetabolitesX(), CCopasiProblem::getModel(), CModel::getModelType(), CModel::getNumDependentReactionMetabs(), CModel::getNumIndependentReactionMetabs(), CModel::getReactions(), CModel::getState(), CModel::getStateTemplate(), CModelEntity::getStatus(), CState::getTime(), CCopasiParameter::getValue(), CCopasiObject::getValuePointer(), CModelEntity::getValueReference(), CRandom::initialize(), mA0, mAmu, mAvgDX, mDoCorrection, mEpsilon, mFirstReactionSpeciesIndex, mMaxReactionFiring, mMaxSteps, mMaxStepsReached, mMethodState, mNextReactionIndex, mNextReactionTime, mNumReactions, mNumReactionSpecies, mNumSpecies, mPartitionedAmu, mPartitionedDependencies, mPartitionedReactionFiring, CTrajectoryMethod::mpCurrentState, mpMethodSpecies, mpModel, CTrajectoryMethod::mpProblem, mpRandomGenerator, mReactionDependencies, mReactionFiring, mSigDX, mSpeciesAfterTau, mSSAStepCounter, CCopasiParameter::Value::pBOOL, CCopasiParameter::Value::pDOUBLE, CCopasiParameter::Value::pINT, CCopasiParameter::Value::pUINT, CModelEntity::REACTIONS, CVector< CType >::resize(), CModel::setState(), CCopasiVector< T >::size(), and CModel::updateSimulatedValues().

195 {
196  /* get configuration data */
197  mMaxSteps = * getValue("Max Internal Steps").pINT;
198  mEpsilon = * getValue("Epsilon").pDOUBLE;
199 
200  bool useRandomSeed = * getValue("Use Random Seed").pBOOL;
201  unsigned C_INT32 randomSeed = * getValue("Random Seed").pUINT;
202 
203  if (useRandomSeed) mpRandomGenerator->initialize(randomSeed);
204 
205  //mpCurrentState is initialized. This state is not used internally in the
206  //stochastic solver, but it is used for returning the result after each step.
207  *mpCurrentState = *initialState;
208 
210  assert(mpModel);
211 
213  mDoCorrection = true;
214  else
215  mDoCorrection = false;
216 
217  //initialize the vector of ints that contains the particle numbers
218  //for the discrete simulation. This also floors all particle numbers in the model.
219 
222 
228 
231  mAmu = 0.0;
232 
234 
236 
239 
240  // Initialize the
241  // Create a local copy of the state where the particle number are rounded to integers
243 
244  const CStateTemplate & StateTemplate = mpModel->getStateTemplate();
245 
246  CModelEntity *const* ppEntity = StateTemplate.beginIndependent();
247  CModelEntity *const* endEntity = StateTemplate.endFixed();
249 
251  size_t Index = 1;
252 
253  C_FLOAT64 * pSpeciesAfterTau = mSpeciesAfterTau.array();
254 
255  for (; ppEntity != endEntity; ++ppEntity, ++pValue, ++Index)
256  {
257  if (dynamic_cast< const CMetab * >(*ppEntity) != NULL)
258  {
259  *pValue = floor(*pValue + 0.5);
260 
261  if ((*ppEntity)->getStatus() == CModelEntity::REACTIONS &&
262  (*ppEntity)->isUsed())
263  {
264  *pSpeciesAfterTau = *pValue;
265  pSpeciesAfterTau++;
266 
268  {
270  mpMethodSpecies = pValue;
271  }
272  }
273  }
274  }
275 
276  // Update the model state so that the species are all represented by integers.
278  mpModel->updateSimulatedValues(false); //for assignments
279 
280  // TODO handle species of type ASSIGNMENT.
281  // These need to be checked whether they are sufficiently close to an integer
282 
283  C_FLOAT64 * pMethodStateValue = mMethodState.beginIndependent() - 1;
284 
285  // Build the reaction dependencies
286  size_t NumReactions = 0;
287 
290  std::vector< CReactionDependencies >::iterator itDependencies = mReactionDependencies.begin();
291 
292  for (; it != end; ++it)
293  {
294  const CCopasiVector<CChemEqElement> & Balances = (*it)->getChemEq().getBalances();
295  const CCopasiVector<CChemEqElement> & Substrates = (*it)->getChemEq().getSubstrates();
296 
297  // This reactions does not change anything we ignore it
298  if (Balances.size() == 0 && Substrates.size() == 0)
299  {
300  continue;
301  }
302 
303  itDependencies->mpParticleFlux = (C_FLOAT64 *)(*it)->getParticleFluxReference()->getValuePointer();
304 
305  itDependencies->mMethodSpeciesIndex.resize(Balances.size());
306  itDependencies->mSpeciesMultiplier.resize(Balances.size());
307  itDependencies->mMethodSpecies.resize(Balances.size());
308  itDependencies->mModelSpecies.resize(Balances.size());
309 
310  std::set< const CCopasiObject * > changed;
311 
312  // The time is always updated
313  changed.insert(mpModel->getValueReference());
314 
317 
318  size_t Index = 0;
319 
320  for (; itBalance != endBalance; ++itBalance)
321  {
322  const CMetab * pMetab = (*itBalance)->getMetabolite();
323 
324  if (pMetab->getStatus() == CModelEntity::REACTIONS)
325  {
326  itDependencies->mMethodSpeciesIndex[Index] = StateTemplate.getIndex(pMetab) - mFirstReactionSpeciesIndex;
327  itDependencies->mSpeciesMultiplier[Index] = floor((*itBalance)->getMultiplicity() + 0.5);
328  itDependencies->mMethodSpecies[Index] = pMethodStateValue + StateTemplate.getIndex(pMetab);
329  itDependencies->mModelSpecies[Index] = (C_FLOAT64 *) pMetab->getValueReference()->getValuePointer();
330 
331  changed.insert(pMetab->getValueReference());
332 
333  Index++;
334  }
335  }
336 
337  // Correct allocation for metabolites which are not determined by reactions
338  itDependencies->mMethodSpeciesIndex.resize(Index, true);
339  itDependencies->mSpeciesMultiplier.resize(Index, true);
340  itDependencies->mMethodSpecies.resize(Index, true);
341  itDependencies->mModelSpecies.resize(Index, true);
342 
343  itDependencies->mSubstrateMultiplier.resize(Substrates.size());
344  itDependencies->mMethodSubstrates.resize(Substrates.size());
345  itDependencies->mModelSubstrates.resize(Substrates.size());
346 
347  CCopasiVector< CChemEqElement >::const_iterator itSubstrate = Substrates.begin();
348  CCopasiVector< CChemEqElement >::const_iterator endSubstrate = Substrates.end();
349 
350  Index = 0;
351 
352  for (; itSubstrate != endSubstrate; ++itSubstrate, ++Index)
353  {
354  const CMetab * pMetab = (*itSubstrate)->getMetabolite();
355 
356  itDependencies->mSubstrateMultiplier[Index] = floor((*itSubstrate)->getMultiplicity() + 0.5);
357  itDependencies->mMethodSubstrates[Index] = pMethodStateValue + StateTemplate.getIndex(pMetab);
358  itDependencies->mModelSubstrates[Index] = (C_FLOAT64 *) pMetab->getValueReference()->getValuePointer();
359  }
360 
361  std::set< const CCopasiObject * > dependend;
362 
364  itDependencies->mDependentReactions.resize(mNumReactions);
365 
366  Index = 0;
367  size_t Count = 0;
368 
369  for (; itReaction != end; ++itReaction, ++Index)
370  {
371  if ((*itReaction)->getParticleFluxReference()->dependsOn(changed))
372  {
373  dependend.insert((*itReaction)->getParticleFluxReference());
374  itDependencies->mDependentReactions[Count] = Index;
375 
376  Count++;
377  }
378  }
379 
380  itDependencies->mDependentReactions.resize(Count, true);
381 
382  itDependencies->mCalculations = CCopasiObject::buildUpdateSequence(dependend, changed);
383 
384  ++itDependencies;
385  ++NumReactions;
386  }
387 
388  mNumReactions = NumReactions;
389 
391  mAmu.resize(mNumReactions, true);
392 
393  mA0 = 0;
394  size_t i;
395 
396  for (i = 0; i < mNumReactions; i++)
397  {
398  mA0 += calculateAmu(i);
399  }
400 
401  mMaxStepsReached = false;
404 
405  mSSAStepCounter = 0;
406 
407  return;
408 }
CVector< size_t > mMaxReactionFiring
virtual size_t size() const
size_t mNumReactionSpecies
iterator begin()
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
CRandom * mpRandomGenerator
CTrajectoryProblem * mpProblem
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
size_t mFirstReactionSpeciesIndex
#define C_INVALID_INDEX
Definition: copasi.h:222
CVector< C_FLOAT64 > mAvgDX
size_t getIndex(const CModelEntity *entity) const
Definition: CState.cpp:231
#define C_INT32
Definition: copasi.h:90
CVector< const C_FLOAT64 * > mPartitionedAmu
Definition: CMetab.h:178
const C_FLOAT64 & calculateAmu(const size_t &index)
CVector< C_FLOAT64 > mReactionFiring
C_FLOAT64 mSSAStepCounter
CVector< C_FLOAT64 > mSpeciesAfterTau
iterator end()
std::vector< CType * >::const_iterator const_iterator
Definition: CCopasiVector.h:57
size_t getNumDependentReactionMetabs() const
Definition: CModel.cpp:1133
std::vector< CReactionDependencies > mReactionDependencies
CModelEntity ** endFixed()
Definition: CState.cpp:213
void setState(const CState &state)
Definition: CModel.cpp:1785
const Value & getValue() const
unsigned C_INT32 * pUINT
CModelEntity ** beginIndependent()
Definition: CState.cpp:208
size_t getNumIndependentReactionMetabs() const
Definition: CModel.cpp:1130
unsigned C_INT32 mMaxSteps
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
CVector< C_FLOAT64 > mAmu
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
const CStateTemplate & getStateTemplate() const
Definition: CModel.cpp:1172
const ModelType & getModelType() const
Definition: CModel.cpp:2339
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
static std::vector< Refresh * > buildUpdateSequence(const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())
virtual void * getValuePointer() const
const CState & getState() const
Definition: CModel.cpp:1771
const CModelEntity::Status & getStatus() const
CVector< const CReactionDependencies * > mPartitionedDependencies
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
C_FLOAT64 mNextReactionTime
CVector< C_FLOAT64 > mSigDX
CModel * getModel() const
CVector< C_FLOAT64 * > mPartitionedReactionFiring
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
virtual void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
Definition: CRandom.cpp:162
CCopasiObject * getValueReference() const
C_FLOAT64 * mpMethodSpecies
CTrajectoryMethod::Status CTrajAdaptiveSA::step ( const double &  deltaT)
virtual

This instructs the method to calculate a time step of deltaT starting with the current state, i.e., the result of the previous step. The new state (after deltaT) is expected in the current state. The return value is the actual timestep taken.

Parameters
const double &deltaT
Returns
Status status

Reimplemented from CTrajectoryMethod.

Definition at line 162 of file CTrajAdaptiveSA.cpp.

References C_FLOAT64, doSingleSSAStep(), doSingleTauLeapStep(), CCopasiMessage::EXCEPTION, CCopasiProblem::getModel(), CModel::getState(), CState::getTime(), MCTrajectoryMethod, mMaxSteps, CTrajectoryMethod::mpCurrentState, CTrajectoryMethod::mpProblem, mSSAStepCounter, CTrajectoryMethod::NORMAL, and CState::setTime().

163 {
164  // do several steps
165  C_FLOAT64 Time = mpCurrentState->getTime();
166  C_FLOAT64 EndTime = Time + deltaT;
167 
168  size_t Steps = 0;
169 
170  while (Time < EndTime)
171  {
172  if (mSSAStepCounter > 0)
173  {
174  Time += doSingleSSAStep(Time, EndTime);
175  mSSAStepCounter--;
176  }
177  else
178  {
179  Time += doSingleTauLeapStep(Time, EndTime);
180  }
181 
182  if (++Steps > mMaxSteps)
183  {
185  }
186  }
187 
189  mpCurrentState->setTime(Time);
190 
191  return NORMAL;
192 }
CTrajectoryProblem * mpProblem
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
C_FLOAT64 mSSAStepCounter
#define MCTrajectoryMethod
C_FLOAT64 doSingleTauLeapStep(const C_FLOAT64 &curTime, const C_FLOAT64 &endTime)
unsigned C_INT32 mMaxSteps
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
#define C_FLOAT64
Definition: copasi.h:92
const CState & getState() const
Definition: CModel.cpp:1771
C_FLOAT64 doSingleSSAStep(const C_FLOAT64 &curTime, const C_FLOAT64 &endTime)
CModel * getModel() const

Friends And Related Function Documentation

Member Data Documentation

C_FLOAT64 CTrajAdaptiveSA::mA0
protected

Total propensity (sum over mAmu[i])

Definition at line 337 of file CTrajAdaptiveSA.h.

Referenced by doSingleSSAStep(), doSingleTauLeapStep(), and start().

CVector< C_FLOAT64 > CTrajAdaptiveSA::mAmu
protected

A vector of reaction propensities

Definition at line 312 of file CTrajAdaptiveSA.h.

Referenced by calculateAmu(), doSingleSSAStep(), doSingleTauLeapStep(), and start().

CVector< C_FLOAT64 > CTrajAdaptiveSA::mAvgDX
private

The mean and variance of species

Definition at line 233 of file CTrajAdaptiveSA.h.

Referenced by doSingleTauLeapStep(), and start().

bool CTrajAdaptiveSA::mDoCorrection
protected

A boolean flag indicating whether correction for higher order reactions need to be applied

Definition at line 302 of file CTrajAdaptiveSA.h.

Referenced by calculateAmu(), and start().

C_FLOAT64 CTrajAdaptiveSA::mEpsilon
private

The tolerance ratio x(t+t')< eps*x(t)

Definition at line 212 of file CTrajAdaptiveSA.h.

Referenced by doSingleTauLeapStep(), and start().

size_t CTrajAdaptiveSA::mFirstReactionSpeciesIndex
private

The Ordered reaction

Definition at line 257 of file CTrajAdaptiveSA.h.

Referenced by start().

CVector< size_t > CTrajAdaptiveSA::mMaxReactionFiring
private

The upper fires of the j-th reactions

Definition at line 222 of file CTrajAdaptiveSA.h.

Referenced by doSingleTauLeapStep(), and start().

unsigned C_INT32 CTrajAdaptiveSA::mMaxSteps
protected

max number of single stochastic steps to do in one step()

Definition at line 287 of file CTrajAdaptiveSA.h.

Referenced by start(), and step().

bool CTrajAdaptiveSA::mMaxStepsReached
protected

A boolean flag indicating whether the maximum steps have been reached. This is used to avoid multiple messages.

Definition at line 343 of file CTrajAdaptiveSA.h.

Referenced by start().

CState CTrajAdaptiveSA::mMethodState
protected

The method internal state which contains particle rounded particle numbers.

Definition at line 322 of file CTrajAdaptiveSA.h.

Referenced by doSingleTauLeapStep(), and start().

size_t CTrajAdaptiveSA::mNextReactionIndex
protected

The index of the next reaction which fires

Definition at line 297 of file CTrajAdaptiveSA.h.

Referenced by doSingleSSAStep(), and start().

C_FLOAT64 CTrajAdaptiveSA::mNextReactionTime
protected

The time the next reaction fires

Definition at line 292 of file CTrajAdaptiveSA.h.

Referenced by doSingleSSAStep(), and start().

size_t CTrajAdaptiveSA::mNumReactions
protected

The particle and reaction numbers

Definition at line 282 of file CTrajAdaptiveSA.h.

Referenced by doSingleSSAStep(), doSingleTauLeapStep(), and start().

size_t CTrajAdaptiveSA::mNumReactionSpecies
protected

Number of variable metabolites.

Definition at line 307 of file CTrajAdaptiveSA.h.

Referenced by doSingleTauLeapStep(), and start().

size_t CTrajAdaptiveSA::mNumSpecies
protected

Definition at line 282 of file CTrajAdaptiveSA.h.

Referenced by start().

CVector< const C_FLOAT64 * > CTrajAdaptiveSA::mPartitionedAmu
protected

The ordered propensity function

Definition at line 317 of file CTrajAdaptiveSA.h.

Referenced by doSingleTauLeapStep(), and start().

CVector< const CReactionDependencies * > CTrajAdaptiveSA::mPartitionedDependencies
protected

The Ordered reaction

Definition at line 332 of file CTrajAdaptiveSA.h.

Referenced by doSingleTauLeapStep(), and start().

CVector< C_FLOAT64 * > CTrajAdaptiveSA::mPartitionedReactionFiring
private

Definition at line 228 of file CTrajAdaptiveSA.h.

Referenced by doSingleTauLeapStep(), and start().

C_FLOAT64* CTrajAdaptiveSA::mpMethodSpecies
private

The species pointer for average, variance, and population (ordered)

Definition at line 242 of file CTrajAdaptiveSA.h.

Referenced by doSingleTauLeapStep(), and start().

CModel* CTrajAdaptiveSA::mpModel
protected

A pointer to the instance of CModel being used.

Definition at line 277 of file CTrajAdaptiveSA.h.

Referenced by doSingleSSAStep(), doSingleTauLeapStep(), and start().

CRandom* CTrajAdaptiveSA::mpRandomGenerator
protected

The random number generator

Definition at line 272 of file CTrajAdaptiveSA.h.

Referenced by doSingleSSAStep(), doSingleTauLeapStep(), start(), and ~CTrajAdaptiveSA().

std::vector< CReactionDependencies > CTrajAdaptiveSA::mReactionDependencies
protected

A vector containing dependency information to minimize the required updates.

Definition at line 327 of file CTrajAdaptiveSA.h.

Referenced by calculateAmu(), doSingleSSAStep(), doSingleTauLeapStep(), and start().

CVector< C_FLOAT64 > CTrajAdaptiveSA::mReactionFiring
private

The number of fires

Definition at line 227 of file CTrajAdaptiveSA.h.

Referenced by doSingleTauLeapStep(), and start().

CVector< C_FLOAT64 > CTrajAdaptiveSA::mSigDX
private

Definition at line 234 of file CTrajAdaptiveSA.h.

Referenced by doSingleTauLeapStep(), and start().

CVector< C_FLOAT64 > CTrajAdaptiveSA::mSpeciesAfterTau
private

The temporary species

Definition at line 247 of file CTrajAdaptiveSA.h.

Referenced by doSingleTauLeapStep(), and start().

C_FLOAT64 CTrajAdaptiveSA::mSSAStepCounter
private

The counter to count the time of doing single SSA

Definition at line 217 of file CTrajAdaptiveSA.h.

Referenced by doSingleTauLeapStep(), start(), and step().

std::vector< Refresh * > CTrajAdaptiveSA::mTauCalculations
private

Vector of refresh methods which need to be executed to update all values required for simulation

Definition at line 252 of file CTrajAdaptiveSA.h.


The documentation for this class was generated from the following files: