COPASI API  4.16.103
CSteadyStateMethod.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2013 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2002 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 /**
16  * CSteadyStateMethod class.
17  * This class describes the interface to all steady state methods.
18  * The various method like Newton have to be derived from
19  * this class.
20  *
21  * Created for COPASI by Stefan Hoops 2002
22  */
23 
24 #include "copasi.h"
25 
27 #include "CSteadyStateMethod.h"
28 #include "CSteadyStateProblem.h"
29 #include "CSteadyStateTask.h"
30 #include "CEigen.h"
31 
32 #include "model/CModel.h"
33 #include "model/CState.h"
34 #include "model/CCompartment.h"
35 
38 {
39  CSteadyStateMethod * pMethod = NULL;
40 
41  switch (subType)
42  {
43  case unset:
44  case Newton:
45  pMethod = new CNewtonMethod();
46  break;
47 
48  default:
49  fatalError();
50  }
51 
52  return pMethod;
53 }
54 
55 /**
56  * Default constructor.
57  */
59  const CCopasiContainer * pParent):
60  CCopasiMethod(CCopasiTask::steadyState, subType, pParent),
61  mpProblem(NULL)
62 {
65 }
66 
67 /**
68  * Copy constructor.
69  * @param "const CSteadyStateMethod &" src
70  */
72  const CCopasiContainer * pParent):
73  CCopasiMethod(src, pParent),
74  mpProblem(src.mpProblem)
75 {
78 }
79 
80 /**
81  * Destructor.
82  */
85 
87 {
88  CCopasiParameter *pParm;
89 
90  assertParameter("Resolution", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-009);
91  mpSSResolution = (C_FLOAT64*)getValue("Resolution").pUDOUBLE;
93 
94  assertParameter("Derivation Factor", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-003);
95  mpDerivationFactor = (C_FLOAT64*)getValue("Derivation Factor").pUDOUBLE;
96 
97  // Check whether we have a method with the old parameter names
98  if ((pParm = getParameter("Newton.DerivationFactor")) != NULL)
99  {
100  setValue("Derivation Factor", *pParm->getValue().pUDOUBLE);
101  removeParameter("Newton.DerivationFactor");
102  }
103 
104  if ((pParm = getParameter("Newton.Resolution")) != NULL)
105  {
106  setValue("Resolution", *pParm->getValue().pUDOUBLE);
107  removeParameter("Newton.Resolution");
108  }
109 }
110 
112 {
114  return true;
115 }
116 
117 /**
118  * This instructs the method to calculate a the steady state
119  * starting with the initialState given.
120  * The steady state is returned in the object pointed to by steadyState.
121  * @param CState * steadyState
122  * @param const CState * initialState
123  * @param C_FLOAT64 * jacobian
124  * @param CEigen * eigenValues
125  * @return CSteadyStateMethod::ReturnCode returnCode
126  */
129  CMatrix< C_FLOAT64 > & jacobianX,
130  CProcessReport * handler)
131 {
132  mpParentTask = dynamic_cast<CSteadyStateTask *>(getObjectParent());
133  assert(mpParentTask);
134 
135  mpSteadyState = pState;
136  mpJacobianX = & jacobianX;
137  mpCallBack = handler;
138 
139  return processInternal();
140 }
141 
142 /**
143  * This function has to be called at the return of any implementation
144  * of the protected function process
145  * @param bool success
146  * @param const C_FLOAT64 & factor
147  * @param const C_FLOAT64 & resolution
148  * @return CSteadyStateMethod::ReturnCode returnCode
149  */
151 CSteadyStateMethod::returnProcess(bool steadyStateFound)
152 {
153  if (!steadyStateFound)
155 
156  if (!allPositive())
158 
161 
163 }
164 
165 /**
166  * This instructs the method to calculate a the steady state
167  * @return CSteadyStateMethod::ReturnCode returnCode
168  */
172 
173 bool CSteadyStateMethod::isEquilibrium(const C_FLOAT64 & resolution) const
174 {
177 
178  for (; it != end; ++it)
179  {
180  const CCompartment * pCompartment = (*it)->getLargestCompartment();
181 
182  if (pCompartment != NULL &&
183  (*it)->getFlux() / pCompartment->getValue() > resolution)
184  return false; //TODO: smallest or largest ?
185  }
186 
187  return true;
188 }
189 
191 {
192  // Assure that all values are updated.
194 
195  CModelEntity *const* ppEntity = mpModel->getStateTemplate().beginIndependent();
196  const C_FLOAT64 * pIt = mpModel->getState().beginIndependent();
197  const C_FLOAT64 * pEnd = mpModel->getState().endDependent();
198 
199  // Skip Model quantities of type ODE
200  for (; pIt != pEnd; ++pIt, ++ppEntity)
201  if (dynamic_cast< const CCompartment *>(*ppEntity) != NULL ||
202  dynamic_cast< const CMetab *>(*ppEntity) != NULL)
203  break;
204 
205  // For all compartments of type ODE we check that the volume is positive
206  for (; pIt != pEnd; ++pIt, ++ppEntity)
207  {
208  if (dynamic_cast< const CCompartment *>(*ppEntity) == NULL)
209  break;
210 
211  if (*pIt < - *mpDerivationResolution)
212  return false;
213  }
214 
215  // We need to check that all metabolites have positive particle numbers
216  // with respect to the given resolution.
217  C_FLOAT64 ParticleResolution =
219 
220  for (; pIt != pEnd; ++pIt, ++ppEntity)
221  {
222  if (dynamic_cast< const CMetab *>(*ppEntity) == NULL)
223  break;
224 
225  if (*pIt < ParticleResolution)
226  return false;
227  }
228 
229  // For all compartments of type ASSIGNMENT we check that the volume is positive
230  for (; pIt != pEnd; ++pIt, ++ppEntity)
231  {
232  if (dynamic_cast< CCompartment *>(*ppEntity) == NULL)
233  break;
234 
235  if (*pIt < - *mpDerivationResolution)
236  return false;
237  }
238 
239  return true;
240 }
241 
242 //virtual
244 {
245  if (!CCopasiMethod::isValidProblem(pProblem)) return false;
246 
247  const CSteadyStateProblem * pP = dynamic_cast<const CSteadyStateProblem *>(pProblem);
248 
249  if (!pP)
250  {
251  //not a TrajectoryProblem
252  CCopasiMessage(CCopasiMessage::EXCEPTION, "Problem is not a steady state problem.");
253  return false;
254  }
255 
256  return true;
257 }
258 
260 {
261  mpProblem = pProblem;
263 
264  return true;
265 }
266 
268  CMatrix< C_FLOAT64 > & jacobianX)
269 {
272 
275 }
276 
278 {
279  //C_FLOAT64* pTmp = (C_FLOAT64*)getValue("Stability Resolution").pUDOUBLE;
280  C_FLOAT64* pTmp = (C_FLOAT64*)getValue("Resolution").pUDOUBLE;
281  assert(pTmp);
282  return *pTmp;
283 }
284 
286 {
290  std::min(*mpDerivationFactor, oldMaxRate),
292 }
293 
295 {
296  return mMethodLog.str();
297 }
C_FLOAT64 * mpDerivationFactor
CMatrix< C_FLOAT64 > * mpJacobianX
virtual bool isValidProblem(const CCopasiProblem *pProblem)
iterator begin()
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
void calculateJacobianX(CMatrix< C_FLOAT64 > &jacobianX, const C_FLOAT64 &derivationFactor, const C_FLOAT64 &resolution)
Definition: CModel.cpp:2082
#define fatalError()
Definition: CState.h:305
std::ostringstream mMethodLog
void calculateJacobian(CMatrix< C_FLOAT64 > &jacobian, const C_FLOAT64 &derivationFactor, const C_FLOAT64 &resolution)
Definition: CModel.cpp:2006
virtual bool isValidProblem(const CCopasiProblem *pProblem)
C_FLOAT64 * endDependent()
Definition: CState.cpp:331
std::string getMethodLog() const
bool removeParameter(const std::string &name)
#define DESTRUCTOR_TRACE
Definition: copasi.h:206
iterator end()
const C_FLOAT64 & getQuantity2NumberFactor() const
Definition: CModel.cpp:2354
void setState(const CState &state)
Definition: CModel.cpp:1785
void calculateJacobianX(const C_FLOAT64 &oldMaxRate)
const Value & getValue() const
bool setValue(const std::string &name, const CType &value)
virtual CSteadyStateMethod::ReturnCode returnProcess(bool steadyStateFound)
CSteadyStateTask * mpParentTask
CModelEntity ** beginIndependent()
Definition: CState.cpp:208
CCopasiParameter * getParameter(const std::string &name)
virtual bool initialize(const CSteadyStateProblem *pProblem)
const CSteadyStateProblem * mpProblem
C_FLOAT64 * mpDerivationResolution
static CSteadyStateMethod * createMethod(CCopasiMethod::SubType subType=CCopasiMethod::Newton)
virtual CSteadyStateMethod::ReturnCode processInternal()
const C_FLOAT64 & getValue() const
#define C_FLOAT64
Definition: copasi.h:92
const CStateTemplate & getStateTemplate() const
Definition: CModel.cpp:1172
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
CSteadyStateMethod::ReturnCode process(CState *pState, CMatrix< C_FLOAT64 > &jacobianX, CProcessReport *handler)
C_FLOAT64 * mpSSResolution
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
C_FLOAT64 getStabilityResolution()
void doJacobian(CMatrix< C_FLOAT64 > &jacobian, CMatrix< C_FLOAT64 > &jacobianX)
const CState & getState() const
Definition: CModel.cpp:1771
bool isEquilibrium(const C_FLOAT64 &resolution) const
CProcessReport * mpCallBack
CModel * getModel() const
#define min(a, b)
Definition: f2c.h:175
CCopasiContainer * getObjectParent() const
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
#define CONSTRUCTOR_TRACE
Definition: copasi.h:202