COPASI API  4.16.103
CScanMethod.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2015 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) 2003 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 /**
16  * CScanMethod class.
17  * This class describes the Scan method
18  *
19  * Created for COPASI by Rohan Luktuke 2002
20  */
21 
22 #include <string>
23 #include <cmath>
24 
25 #include "copasi.h"
26 #include "model/CModel.h"
27 #include "model/CState.h"
28 #include "utilities/CReadConfig.h"
30 //#include "utilities/CWriteConfig.h"
31 #include "CScanProblem.h"
32 #include "CScanMethod.h"
33 #include "CScanTask.h"
34 
38 // this will have to be defined somewhere else with the
39 // values of other distribution types
40 //#define SD_UNIFORM 0
41 //#define SD_GAUSS 1
42 //#define SD_BOLTZ 2
43 //#define SD_REGULAR 3
44 
45 //**************** CScanItem classes ***************************
46 
47 //static
49  CRandom* rg,
50  const bool & continueFromCurrentState)
51 {
52  if (!si) return NULL;
53 
54  CScanProblem::Type type = *(CScanProblem::Type*)(si->getValue("Type").pUINT);
55 
56  CScanItem* tmp = NULL;
57 
58  if (type == CScanProblem::SCAN_REPEAT)
59  tmp = new CScanItemRepeat(si, continueFromCurrentState);
60 
61  if (type == CScanProblem::SCAN_LINEAR)
62  tmp = new CScanItemLinear(si, continueFromCurrentState);
63 
64  if (type == CScanProblem::SCAN_RANDOM)
65  tmp = new CScanItemRandom(si, rg, continueFromCurrentState);
66 
67  /* if (type == CScanProblem::SCAN_BREAK)
68  tmp = new CScanItemBreak(si, st);*/
69 
70  return tmp;
71 }
72 
74  mNumSteps(0),
75  mpObject(NULL),
76  mpInitialObject(NULL),
77  mStoreValue(0.0),
78  mIndex(0),
79  mFlagFinished(false),
80  mIsStateVariable(false)
81 {}
82 
84  const bool & continueFromCurrentState):
85  mNumSteps(0),
86  mpObject(NULL),
87  mpInitialObject(NULL),
88  mStoreValue(0.0),
89  mIndex(0),
90  mFlagFinished(false),
91  mIsStateVariable(false)
92 {
93  assert(si != NULL);
94 
95  mNumSteps = * si->getValue("Number of steps").pUINT;
96 
97  std::string tmpString = * si->getValue("Object").pCN;
98  CCopasiDataModel* pDataModel = si->getObjectDataModel();
99  assert(pDataModel != NULL);
100  CCopasiObject * tmpObject = pDataModel->getDataObject(tmpString);
101 
102  if (tmpObject == NULL || !tmpObject->isValueDbl())
103  {
104  return;
105  }
106 
107  mpInitialObject = tmpObject;
108 
109  if (continueFromCurrentState)
110  {
111  // Determine whether the object is the initial value of a state variable;
112  const CModel * pModel = static_cast< CModel * >(mpInitialObject->getObjectAncestor("Model"));
113 
114  if (pModel != NULL)
115  {
116  mIsStateVariable = pModel->isStateVariable(tmpObject);
117 
118  if (!mIsStateVariable)
119  {
120  // We need to find the object for the transient value if it exists so that we can update while continuing.
121  mpObject = pModel->getCorrespondingTransientObject(tmpObject);
122  }
123  }
124  }
125  else
126  {
128  }
129 }
130 
131 size_t CScanItem::getNumSteps() const {return mNumSteps;};
132 
134 {
135  if (mpObject)
137 };
138 
140 {
141  if (mpObject)
142  {
143  assert(mpObject->isValueDbl());
145  }
146 };
147 
149 {
150  mIndex = 0;
151  mFlagFinished = false;
152  this->step(); //purely virtual
153 }
154 
155 bool CScanItem::isFinished() const {return mFlagFinished;};
156 
157 bool CScanItem::isValidScanItem(const bool & continueFromCurrentState)
158 {
159  if (continueFromCurrentState &&
161  mpObject == NULL)
162  {
164  return true;
165  }
166 
167  if (mpInitialObject == NULL)
168  {
169  CCopasiMessage(CCopasiMessage::ERROR, "Invalid or missing scan parameter.");
170  return false;
171  }
172 
173  return true;
174 }
175 
177 {
178  return mpObject;
179 }
180 
181 //*******
182 
184  const bool & continueFromCurrentState):
185  CScanItem(si, continueFromCurrentState)
186 {
187  if (mNumSteps >= 1)
188  --mNumSteps; // for the repeat item mNumSteps is the number of iterations, not of intervals
189 }
190 
192 {
193  //do something ...
194 
195  //the index
196  if (mIndex > mNumSteps)
197  mFlagFinished = true;
198 
199  ++mIndex;
200 }
201 
202 bool CScanItemRepeat::isValidScanItem(const bool & /* continueFromCurrentState */)
203 {
204  return true;
205 }
206 
207 //*******
208 
210  const bool & continueFromCurrentState):
211  CScanItem(si, continueFromCurrentState),
212  mLog(false)
213 {
214  mLog = * si->getValue("log").pBOOL;
215  mMin = * si->getValue("Minimum").pDOUBLE;
216  mMax = * si->getValue("Maximum").pDOUBLE;
217 
218  if (mLog)
219  {
220  mMin = log(mMin);
221  mMax = log(mMax);
222  }
223 
224  mFaktor = (mMax - mMin) / mNumSteps;
225 
226  //TODO: log scanning of negative values?
227 }
228 
230 {
231  //do something ...
232  C_FLOAT64 Value = mMin + mIndex * mFaktor;
233 
234  if (mLog)
235  Value = exp(Value);
236 
237  //the index
238  if (mIndex > mNumSteps)
239  mFlagFinished = true;
240 
241  if (mpObject) mpObject->setObjectValue(Value);
242 
244 
245  ++mIndex;
246 }
247 
248 bool CScanItemLinear::isValidScanItem(const bool & continueFromCurrentState)
249 {
250  if (!CScanItem::isValidScanItem(continueFromCurrentState)) return false;
251 
252  if (mLog)
253  {
255  {
256  //not a valid range for log
257  CCopasiMessage(CCopasiMessage::ERROR, "Only positive values for min and max are possible for a logarithmic scan.");
258  return false;
259  }
260  }
261 
262  return true;
263 }
264 
265 //*******
266 
268  const bool & continueFromCurrentState):
269  CScanItem(si, continueFromCurrentState),
270  mRg(rg),
271  mRandomType(0),
272  mLog(false)
273 {
274  mRandomType = * si->getValue("Distribution type").pUINT;
275 
276  mLog = * si->getValue("log").pBOOL;
277  mMin = * si->getValue("Minimum").pDOUBLE;
278  mMax = * si->getValue("Maximum").pDOUBLE;
279 
280  if (mLog && mRandomType == 0)
281  {
282  mMin = log(mMin);
283  mMax = log(mMax);
284  }
285 
286  mNumSteps = 0;
287  mFaktor = (mMax - mMin);
288 }
289 
291 {
292  C_FLOAT64 Value;
293 
294  //the index
295  if (mIndex > mNumSteps)
296  mFlagFinished = true;
297  else
298  {
299  C_FLOAT64 tmpF;
300 
301  switch (mRandomType)
302  {
303  case 0: // uniform
304  Value = mMin + mRg->getRandomCC() * mFaktor;
305 
306  if (mLog)
307  Value = exp(Value);
308 
309  break;
310 
311  case 1: // normal
312  tmpF = mRg->getRandomNormal01();
313  Value = mMin + tmpF * mMax;
314 
315  if (mLog)
316  Value = exp(Value);
317 
318  break;
319 
320  case 2: // poisson
321 
322  if (mMin < 0)
323  CCopasiMessage(CCopasiMessage::WARNING, "Invalid ScanItem: Requested Poisson random variable for negative argument: %lf", mMin);
324 
325  Value = mRg->getRandomPoisson(mMin);
326 
327  break;
328 
329  case 3: // gamma
330  Value = mRg->getRandomGamma(mMin, mMax);
331 
332  if (mLog)
333  Value = exp(Value);
334 
335  break;
336  }
337  }
338 
339  if (mpObject) mpObject->setObjectValue(Value);
340 
342 
343  ++mIndex;
344 }
345 
346 //**************** CScanMethod class ***************************
347 
349 {
350  return new CScanMethod();
351 }
352 
354  CCopasiMethod(CCopasiTask::scan, CCopasiMethod::scanMethod),
355  mpProblem(NULL),
356  mpTask(NULL),
357  mpRandomGenerator(NULL),
358  mTotalSteps(1),
359  mLastNestingItem(C_INVALID_INDEX),
360  mContinueFromCurrentState(false)
361 {
363 }
364 
366 {
368  delete mpRandomGenerator;
369  mpRandomGenerator = NULL;
370 }
371 
373 {
374  if (!mpProblem) return false;
375 
376  size_t i, imax = mScanItems.size();
377 
378  for (i = 0; i < imax; ++i) if (mScanItems[i]) delete mScanItems[i];
379 
380  mScanItems.clear();
381  return true;
382 }
383 
385 {
386  if (!mpProblem) return false;
387 
388  mpTask = dynamic_cast< CScanTask * >(getObjectParent());
389 
390  if (mpTask == NULL) return false;
391 
393  mInitialRefreshes.clear();
394  //mTransientRefreshes.clear();
395 
396  mTotalSteps = 1;
397  std::set< const CCopasiObject * > InitialObjectSet;
398  std::set< const CCopasiObject * > TransientObjectSet;
399 
400  size_t i, imax = mpProblem->getNumberOfScanItems();
402 
403  for (i = 0; i < imax; ++i)
404  {
408 
409  if (pItem == NULL)
410  {
411  continue;
412  }
413 
414  mScanItems.push_back(pItem);
415  mTotalSteps *= pItem->getNumSteps() + 1;
416 
417  if (pItem->getObject() != NULL)
418  {
419  InitialObjectSet.insert(pItem->getObject());
420  }
421  }
422 
424  {
426  }
427  else
428  {
430  }
431 
432  //set mLastNestingItem
434 
435  if (imax != 0)
436  {
437  //search from the end
438  size_t j;
439 
440  for (j = mScanItems.size() - 1; j != C_INVALID_INDEX; --j)
441  {
442  if (mScanItems[j]->isNesting())
443  {
444  mLastNestingItem = j;
445  break;
446  }
447  }
448  }
449 
450  return true;
451 }
452 
454 {
455  if (!mpProblem) return false;
456 
457  //a hack to ensure that the first subtask is run with initial conditions
458  //pDataModel->getModel()->setState(&mpProblem->getInitialState());
459 
460  bool success = true;
461 
462  size_t i, imax = mScanItems.size();
463 
464  //store old parameter values
465  for (i = 0; i < imax; ++i)
466  mScanItems[i]->storeValue();
467 
468  //Do the scan...
469  if (imax) //there are scan items
470  success = loop(0);
471  else
472  success = calculate(); //nothing to scan, only one call to the subtask
473 
474  //restore old parameter values
475  for (i = 0; i < imax; ++i)
476  mScanItems[i]->restoreValue();
477 
478  return success;
479 }
480 
481 bool CScanMethod::loop(size_t level)
482 {
483  bool isLastMasterItem = (level == (mScanItems.size() - 1)); //TODO
484 
485  CScanItem* currentSI = mScanItems[level];
486 
487  for (currentSI->reset(); !currentSI->isFinished(); currentSI->step())
488  {
489  //TODO: handle slave SIs
490 
491  if (isLastMasterItem)
492  {
493  if (!calculate()) return false;
494  }
495  else
496  {
497  if (!loop(level + 1)) return false;
498  } //TODO
499 
500  //separator needs to be handled slightly differently if we are at the last item
501  if (currentSI->isNesting())
502  ((CScanTask*)(getObjectParent()))->outputSeparatorCallback(level == mLastNestingItem);
503  }
504 
505  return true;
506 }
507 
509 {
510  std::vector< Refresh * >::iterator it = mInitialRefreshes.begin();
511  std::vector< Refresh * >::iterator end = mInitialRefreshes.end();
512 
513  while (it != end)
514  (**it++)();
515 
516  return mpTask->processCallback();
517 }
518 
520 {mpProblem = problem;}
521 
522 //virtual
524 {
525  if (!CCopasiMethod::isValidProblem(pProblem)) return false;
526 
527  const CScanProblem * pP = dynamic_cast<const CScanProblem *>(pProblem);
528 
529  if (!pP)
530  {
531  //not a TrajectoryProblem
532  CCopasiMessage(CCopasiMessage::EXCEPTION, "Problem is not a Scan problem.");
533  return false;
534  }
535 
537 
538  size_t i, imax = pP->getNumberOfScanItems();
539 
540  if (imax <= 0)
541  {
542  //no scan items
543  CCopasiMessage(CCopasiMessage::WARNING, "There is nothing to scan.");
544  return false;
545  }
546 
547  for (i = 0; i < imax; ++i)
548  {
552 
553  if (!si)
554  {
555  //parameter group could not be interpreted
556  CCopasiMessage(CCopasiMessage::ERROR, "Internal problem with scan definition.");
557  return false;
558  }
559 
561  {
562  //the self check of the scan item failed.
563  //the message should be generated by the isValidScanItem() method.
564  delete si;
565  return false;
566  }
567 
568  delete si;
569  //mTotalSteps *= mScanItems[i]->getNumSteps();
570  }
571 
572  return true;
573 }
CCopasiDataModel * getObjectDataModel()
CCopasiContainer * getObjectAncestor(const std::string &type) const
CCopasiObject * getDataObject(const CCopasiObjectName &CN) const
virtual std::string getObjectDisplayName(bool regular=true, bool richtext=false) const
CCopasiObject * getCorrespondingTransientObject(const CCopasiObject *pObject) const
Definition: CModel.cpp:4032
C_FLOAT64 mMax
Definition: CScanMethod.h:131
void setObjectValue(const C_FLOAT64 &value)
virtual bool isValidScanItem(const bool &continueFromCurrentState)
void setProblem(CScanProblem *problem)
C_FLOAT64 mMax
Definition: CScanMethod.h:115
static CScanMethod * createMethod(CCopasiMethod::SubType subType=CCopasiMethod::scanMethod)
bool isFinished() const
CCopasiObject * mpInitialObject
Definition: CScanMethod.h:42
virtual bool isValidProblem(const CCopasiProblem *pProblem)
size_t mLastNestingItem
Definition: CScanMethod.h:189
bool mIsStateVariable
Definition: CScanMethod.h:50
bool calculate()
C_FLOAT64 mMin
Definition: CScanMethod.h:115
CScanItemRandom(CCopasiParameterGroup *si, CRandom *rg, const bool &continueFromCurrentState)
virtual void step()
const std::set< const CCopasiObject * > & getUptoDateObjects() const
Definition: CModel.cpp:1178
std::vector< CScanItem * > mScanItems
Definition: CScanMethod.h:177
#define C_INVALID_INDEX
Definition: copasi.h:222
virtual bool isValidScanItem(const bool &continueFromCurrentState)
size_t getNumberOfScanItems() const
bool cleanupScanItems()
unsigned C_INT32 mRandomType
Definition: CScanMethod.h:133
CScanProblem * mpProblem
Definition: CScanMethod.h:165
CRegisteredObjectName * pCN
virtual bool isNesting() const
Definition: CScanMethod.h:70
std::vector< Refresh * > mInitialRefreshes
Definition: CScanMethod.h:179
virtual bool isValidProblem(const CCopasiProblem *pProblem)
virtual bool isValidScanItem(const bool &continueFromCurrentState)
C_FLOAT64 mStoreValue
Definition: CScanMethod.h:44
virtual C_FLOAT64 getRandomNormal01()
Definition: CRandom.cpp:261
virtual void step()
size_t mNumSteps
Definition: CScanMethod.h:39
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
virtual C_FLOAT64 getRandomPoisson(const C_FLOAT64 &mean)
Definition: CRandom.cpp:314
void storeValue()
virtual C_FLOAT64 getRandomCC()
Definition: CRandom.cpp:235
bool processCallback()
Definition: CScanTask.cpp:201
CScanItemRepeat(CCopasiParameterGroup *si, const bool &continueFromCurrentState)
bool mFlagFinished
Definition: CScanMethod.h:48
C_FLOAT64 mMin
Definition: CScanMethod.h:131
#define MCScan
unsigned C_INT32 * pUINT
size_t mIndex
Definition: CScanMethod.h:46
bool getContinueFromCurrentState() const
CScanTask * mpTask
Definition: CScanMethod.h:170
CRandom * mRg
Definition: CScanMethod.h:132
virtual void step()=0
#define C_FLOAT64
Definition: copasi.h:92
CRandom * mpRandomGenerator
Definition: CScanMethod.h:175
CScanItemLinear(CCopasiParameterGroup *si, const bool &continueFromCurrentState)
void reset()
const CCopasiObject * getObject() const
static std::vector< Refresh * > buildUpdateSequence(const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())
bool isValueDbl() const
size_t getNumSteps() const
const CCopasiParameter::Value & getValue(const std::string &name) const
virtual void * getValuePointer() const
virtual void step()
const CCopasiParameterGroup * getScanItem(size_t index) const
virtual C_FLOAT64 getRandomGamma(C_FLOAT64 shape, C_FLOAT64 scale)
Definition: CRandom.cpp:572
Definition: CModel.h:50
std::vector< Refresh * > buildInitialRefreshSequence(std::set< const CCopasiObject * > &changedObjects)
Definition: CModel.cpp:4164
void restoreValue() const
static CScanItem * createScanItemFromParameterGroup(CCopasiParameterGroup *si, CRandom *rg, const bool &continueFromCurrentState)
Definition: CScanMethod.cpp:48
bool loop(size_t level)
C_FLOAT64 mFaktor
Definition: CScanMethod.h:115
CModel * getModel() const
CCopasiObject * mpObject
Definition: CScanMethod.h:41
CCopasiContainer * getObjectParent() const
bool isStateVariable(const CCopasiObject *pObject) const
Definition: CModel.cpp:3983
bool mContinueFromCurrentState
Definition: CScanMethod.h:195
size_t mTotalSteps
Definition: CScanMethod.h:183
C_FLOAT64 mFaktor
Definition: CScanMethod.h:131
#define max(a, b)
Definition: f2c.h:176