COPASI API  4.16.103
CDimension.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) 2006 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include <sstream>
19 #include "copasi/model/CModel.h"
20 #include "copasi/model/CChemEq.h"
21 
23  : mD1(0), mD2(0), mD3(0), mD4(0), mD5(0),
24  mUnknown(true),
25  mContradiction(false)
26 {}
27 
29 {
30  mUnknown = true;
31  mContradiction = false;
32 }
34 {
35  return mUnknown;
36 }
37 
39 {
40  mUnknown = false;
41  mContradiction = true;
42 }
44 {
45  return mContradiction;
46 }
47 
48 void CDimension::setDimension(const C_FLOAT64 & d1, const C_FLOAT64 & d2, const C_FLOAT64 & d3,
49  const C_FLOAT64 & d4, const C_FLOAT64 & d5)
50 {
51  mUnknown = false;
52  mContradiction = false;
53  mD1 = d1;
54  mD2 = d2;
55  mD3 = d3;
56  mD4 = d4;
57  mD5 = d5;
58 }
59 
60 //static
61 std::string CDimension::constructDisplayElement(const std::string & base, C_FLOAT64 exponent)
62 {
63  if (exponent <= 0) return "";
64 
65  if (exponent == 1.0) return base;
66 
67  std::ostringstream ss;
68  ss << base << "^" << exponent;
69  return ss.str();
70 }
71 
72 std::string CDimension::getDisplayString(const CModel * pModel) const
73 {
74  if (isUnknown())
75  return "?";
76 
77  if (isContradiction())
78  return "EEE";
79 
80  assert(pModel != NULL);
81 
82  std::string vol = pModel->getVolumeUnitName();
83  std::string time = pModel->getTimeUnitName();
84  std::string quan = pModel->getQuantityUnitName();
85  std::string area = pModel->getAreaUnitName();
86  std::string len = pModel->getLengthUnitName();
87 
88  std::string tmp;
89 
90  //positive exponents
91  std::string s1;
92  s1 = constructDisplayElement(quan, mD1);
93 
94  tmp = constructDisplayElement(vol, mD2);
95 
96  if ((s1 != "") && (tmp != "")) s1 += "*";
97 
98  s1 += tmp;
99 
100  tmp = constructDisplayElement(time, mD3);
101 
102  if ((s1 != "") && (tmp != "")) s1 += "*";
103 
104  s1 += tmp;
105 
106  tmp = constructDisplayElement(area, mD4);
107 
108  if ((s1 != "") && (tmp != "")) s1 += "*";
109 
110  s1 += tmp;
111 
112  tmp = constructDisplayElement(len, mD5);
113 
114  if ((s1 != "") && (tmp != "")) s1 += "*";
115 
116  s1 += tmp;
117 
118  //negative exponents
119  std::string s2;
120  bool parflag = false;
121  s2 = constructDisplayElement(quan, -mD1);
122 
123  tmp = constructDisplayElement(vol, -mD2);
124 
125  if ((s2 != "") && (tmp != ""))
126  {s2 += "*"; parflag = true;}
127 
128  s2 += tmp;
129 
130  tmp = constructDisplayElement(time, -mD3);
131 
132  if ((s2 != "") && (tmp != ""))
133  {s2 += "*"; parflag = true;}
134 
135  s2 += tmp;
136 
137  tmp = constructDisplayElement(area, -mD4);
138 
139  if ((s2 != "") && (tmp != ""))
140  {s2 += "*"; parflag = true;}
141 
142  s2 += tmp;
143 
144  tmp = constructDisplayElement(len, -mD5);
145 
146  if ((s2 != "") && (tmp != ""))
147  {s2 += "*"; parflag = true;}
148 
149  s2 += tmp;
150 
151  if (parflag) s2 = "(" + s2 + ")";
152 
153  //put it together..
154  if ((s1 == "") && (s2 == ""))
155  return "1";
156 
157  if (s1 == "")
158  return "1/" + s2;
159 
160  if (s2 == "")
161  return s1;
162 
163  return s1 + "/" + s2;
164 }
165 
166 bool CDimension::operator==(const CDimension & rhs) const
167 {
168  return (mUnknown == rhs.mUnknown)
169  && (mContradiction == rhs.mContradiction)
170  && (mD1 == rhs.mD1)
171  && (mD2 == rhs.mD2)
172  && (mD3 == rhs.mD3)
173  && (mD4 == rhs.mD4)
174  && (mD5 == rhs.mD5);
175 }
176 
178 {
179  CDimension result;
180 
181  if (isContradiction() || rhs.isContradiction())
182  result.setContradiction();
183  else if (isUnknown() || rhs.isUnknown())
184  result.setUnknown();
185  else
186  result.setDimension(mD1 + rhs.mD1, mD2 + rhs.mD2, mD3 + rhs.mD3, mD4 + rhs.mD4, mD5 + rhs.mD5);
187 
188  return result;
189 }
190 
192 {
193  CDimension result;
194 
195  if (isContradiction() || rhs.isContradiction())
196  result.setContradiction();
197  else if (isUnknown() || rhs.isUnknown())
198  result.setUnknown();
199  else
200  result.setDimension(mD1 - rhs.mD1, mD2 - rhs.mD2, mD3 - rhs.mD3, mD4 - rhs.mD4, mD5 - rhs.mD5);
201 
202  return result;
203 }
204 
206 {
207  CDimension result;
208 
209  if (isContradiction())
210  result.setContradiction();
211  else if (isUnknown())
212  result.setUnknown();
213  else
214  result.setDimension(mD1 * rhs, mD2 * rhs, mD3 * rhs, mD4 * rhs, mD5 * rhs);
215 
216  return result;
217 }
218 
220 {
221  CDimension result;
222 
223  if (this->isContradiction() || rhs.isContradiction())
224  result.setContradiction();
225  else if (*this == rhs)
226  result = *this;
227  else if (this->isUnknown())
228  result = rhs;
229  else if (rhs.isUnknown())
230  result = *this;
231  else
232  result.setContradiction();
233 
234  return result;
235 }
236 
237 /**
238  * Disabled becuase the in order to generate output, the dimensions instance
239  * needs the datamodel.
240 std::ostream & operator<<(std::ostream &os, const CDimension & d)
241 {
242  if (d.mUnknown) os << "Dim: unknown";
243  else if (d.mContradiction) os << "Dim: conctradiction";
244  else os << "Dim: (" << d.mD1 << ", " << d.mD2 << ", " << d.mD3 << ") " << d.getDisplayString();
245 
246  return os;
247 }
248  */
249 
250 std::string CDimension::print(const CModel* pModel) const
251 {
252  std::ostringstream os;
253 
254  if (this->mUnknown) os << "Dim: unknown";
255  else if (this->mContradiction) os << "Dim: contradiction";
256  else os << "Dim: (" << this->mD1 << ", " << this->mD2 << ", " << this->mD3
257  << ", " << this->mD4 << ", " << this->mD5 << ") "
258  << this->getDisplayString(pModel);
259 
260  return os.str();
261 }
262 
263 void CDimension::fixDimensionless(bool d1, bool d2, bool d3, bool d4, bool d5)
264 {
265  if (d1)
266  mD1 = 0;
267 
268  if (d2)
269  mD2 = 0;
270 
271  if (d3)
272  mD3 = 0;
273 
274  if (d4)
275  mD4 = 0;
276 
277  if (d5)
278  mD5 = 0;
279 }
280 
281 //*************************************************************************
282 
283 #include "function/CFunction.h"
284 
285 CFindDimensions::CFindDimensions(const CFunction* function, bool d1, bool d2, bool d3, bool d4, bool d5)
286  : mpFunction(function),
287  mRootDimension(),
288  mUseHeuristics(false),
289  mM1(-1.0), mM2(-1.0),
290  mD1(d1), mD2(d2), mD3(d3), mD4(d4), mD5(d5)
291 {
292  setupDimensions();
293 }
294 
296 {
297  if (!mpFunction) return;
298 
300 
301  size_t i, imax = mpFunction->getVariables().size();
302 
303  for (i = 0; i < imax; ++i)
304  {
305  switch (mpFunction->getVariables()[i]->getUsage())
306  {
310  mDimensions[i].setDimension(1, -1, 0, 0, 0); //concentration
311  break;
312 
314  mDimensions[i].setUnknown(); // TODO Dimension(0, 1, 0); //volume
315  break;
316 
318  mDimensions[i].setDimension(0, 0, 1, 0, 0); //time
319  break;
320 
321  default:
322  mDimensions[i].setUnknown();
323  break;
324  }
325 
326  mDimensions[i].fixDimensionless(mD1, mD2, mD3, mD4, mD5);
327  }
328 }
329 
331 {
333 }
334 
335 const std::vector<CDimension> & CFindDimensions::getDimensions() const
336 {
337  return mDimensions;
338 }
339 
341 {
342  mRootDimension = rootDim;
343  findDimensions();
344 }
345 
346 void CFindDimensions::findDimensions(bool isMulticompartment)
347 {
348  if (isMulticompartment)
349  mRootDimension.setDimension(1, 0, -1, 0, 0); //amount of subs/time
350  else
351  mRootDimension.setDimension(1, -1, -1, 0, 0); //TODO !!! conc/time
352 
354 
355  findDimensions();
356 }
357 
358 std::vector<std::string> CFindDimensions::findDimensionsBoth(const CModel* pModel)
359 {
360  //first for single compartment
361  findDimensions(false);
362  std::vector<CDimension> store = mDimensions;
363 
364  //next for multiple compartments
365  setupDimensions();
366  findDimensions(true);
367 
368  //compare...
369  std::vector<std::string> ret;
370  std::vector<CDimension>::const_iterator it1, it2, it1end = store.end();
371 
372  for (it1 = store.begin(), it2 = mDimensions.begin(); it1 != it1end; ++it1, ++it2)
373  {
374  if (*it1 == *it2)
375  ret.push_back(it1->getDisplayString(pModel));
376  else
377  ret.push_back(it1->getDisplayString(pModel) + " or " + it2->getDisplayString(pModel));
378  }
379 
380  return ret;
381 }
382 
384 {
385  if (!mpFunction) return;
386 
387  if (dynamic_cast<const CMassAction*>(mpFunction))
388  {
390  return;
391  }
392 
393  size_t i, imax = mpFunction->getVariables().size();
394 
395  for (i = 0; i < imax; ++i)
396  if (mDimensions[i].isUnknown())
397  findDimension(i);
398 
399  for (i = 0; i < imax; ++i)
400  if (mDimensions[i].isUnknown())
401  findDimension(i);
402 
403  for (i = 0; i < imax; ++i)
404  if (mDimensions[i].isUnknown())
405  findDimension(i);
406 
407  //TODO: conistency check for known dimensions?
408 }
409 
411 {
412  //mpChemEq = eq;
413  if (!eq)
414  {
415  mM1 = 0.0; mM2 = 0.0;
416  return;
417  }
418 
421 }
422 
424  const size_t & m2)
425 {
426  mM1 = (m1 != C_INVALID_INDEX) ? m1 : -1.0;
427  mM2 = (m2 != C_INVALID_INDEX) ? m2 : -1.0;
428 }
429 
431 {
432  if (mM1 < 0) return;
433 
434  CDimension conc; conc.setDimension(1.0, -1.0, 0.0, 0, 0); //TODO!!!
435 
437  conc.fixDimensionless(mD1, mD2, mD3, mD4, mD5);
438 
439  if (mDimensions[0].isUnknown())
440  {
441  mDimensions[0] = mRootDimension - conc * mM1;
442  }
443 
444  if (mDimensions.size() == 2) return; //irreversible
445 
446  if (mDimensions[2].isUnknown())
447  {
448  mDimensions[2] = mRootDimension - conc * mM2;
449  }
450 }
451 
452 //find dim for one parameter
454 {
455  if (!mpFunction) return;
456 
457  if (index >= mDimensions.size()) return;
458 
459  CDimension result;
460 
461  //find all variable nodes with given index
462  std::vector<const CEvaluationNode*> nodes;
463  const std::vector< CEvaluationNode * > & allnodes = mpFunction->getNodeList();
464  std::vector< CEvaluationNode * >::const_iterator it, itEnd = allnodes.end();
465 
466  for (it = allnodes.begin(); it != itEnd; ++it)
467  {
468  if ((*it)->getType() == CEvaluationNode::VARIABLE)
469  if (dynamic_cast<const CEvaluationNodeVariable*>(*it)->getIndex() == index)
470  {
471  nodes.push_back(*it);
472  }
473  }
474 
475  //find dimension for all nodes and compare results
476  std::vector<const CEvaluationNode * >::const_iterator it2, it2End = nodes.end();
477 
478  for (it2 = nodes.begin(); it2 != it2End; ++it2)
479  {
480  result = result.compare(findDimension(*it2));
481  }
482 
483  mDimensions[index] = result;
484 }
485 
487  const CEvaluationNode * requestingNode)
488 {
489  CDimension result;
490 
491  //variable node
492  const CEvaluationNodeVariable* varnode = dynamic_cast<const CEvaluationNodeVariable*>(node);
493 
494  if (varnode)
495  {
496  //is dim known?
497  if (!mDimensions[varnode->getIndex()].isUnknown())
498  result = mDimensions[varnode->getIndex()];
499 
500  else if (requestingNode) //do not try to evaluate dim recursively if asked by parent node
501  {
502  if (node->getParent())
503  result.setUnknown();
504  else //no parent
505  fatalError();
506  }
507  else // !requestingNode; method is called from outside, evaluate recursively
508  {
509  if (node->getParent())
510  {
511  const CEvaluationNode* parent = dynamic_cast<const CEvaluationNode*>(node->getParent());
512 
513  if (parent) result = findDimension(parent, node);
514  }
515  else //no parent, use root dimension
516  {
517  result = mRootDimension;
518  }
519  }
520  }
521 
522  //operator node
523  const CEvaluationNodeOperator* opnode = dynamic_cast<const CEvaluationNodeOperator*>(node);
524 
525  if (opnode)
526  {
527  assert(requestingNode); //should not be called from outside
528 
529  switch (opnode->getType() & 0x00FFFFFF)
530  {
533  {
534  CDimension r1, r2;
535 
536  if (requestingNode == node->getParent()) //called by parent
537  {
538  r1 = findDimension(opnode->getLeft(), node);
539  r2 = findDimension(opnode->getRight(), node);
540  }
541  else //called by one of the children
542  {
543  const CEvaluationNode* parent = dynamic_cast<const CEvaluationNode*>(node->getParent());
544 
545  if (parent) r1 = findDimension(parent, node);
546  else r1 = mRootDimension;
547 
548  if (requestingNode == opnode->getLeft())
549  r2 = findDimension(opnode->getRight(), node);
550 
551  if (requestingNode == opnode->getRight())
552  r2 = findDimension(opnode->getLeft(), node);
553  }
554 
555  //r1 and r2 are now the dimensions of the two nodes that are not the
556  //requesting node.
557  /*if (r1.isContradiction() || r2.isContradiction())
558  result.setContradiction();
559  else if (r1 == r2)
560  result = r1;
561  else if (r1.isUnknown())
562  result = r2;
563  else if (r2.isUnknown())
564  result =r1;
565  else
566  result.setContradiction();*/
567  result = r1.compare(r2);
568  }
569  break;
570 
572  {
573  CDimension r1, r2;
574 
575  if (requestingNode == node->getParent()) //called by parent
576  {
577  r1 = findDimension(opnode->getLeft(), node);
578  r2 = findDimension(opnode->getRight(), node);
579 
580  result = r1 + r2;
581  }
582  else //called by one of the children
583  {
584  //parent node
585  const CEvaluationNode* parent = dynamic_cast<const CEvaluationNode*>(node->getParent());
586 
587  if (parent) r1 = findDimension(parent, node);
588  else r1 = mRootDimension;
589 
590  //other child
591  if (requestingNode == opnode->getLeft())
592  r2 = findDimension(opnode->getRight(), node);
593 
594  if (requestingNode == opnode->getRight())
595  r2 = findDimension(opnode->getLeft(), node);
596 
597  result = r1 - r2;
598  }
599  }
600  break;
601 
603  {
604  CDimension r1, r2;
605 
606  if (requestingNode == node->getParent()) //called by parent
607  {
608  r1 = findDimension(opnode->getLeft(), node);
609  r2 = findDimension(opnode->getRight(), node);
610 
611  result = r1 - r2;
612  }
613  else //called by one of the children
614  {
615  //parent node
616  const CEvaluationNode* parent = dynamic_cast<const CEvaluationNode*>(node->getParent());
617 
618  if (parent) r1 = findDimension(parent, node);
619  else r1 = mRootDimension;
620 
621  //other child
622  if (requestingNode == opnode->getLeft())
623  {
624  r2 = findDimension(opnode->getRight(), node);
625  result = r1 + r2;
626  }
627 
628  if (requestingNode == opnode->getRight())
629  {
630  r2 = findDimension(opnode->getLeft(), node);
631  result = r2 - r1;
632  }
633  }
634  }
635  break;
636 
637  default:
638  break;
639  }
640  }
641 
642  //number node
643  const CEvaluationNodeNumber* numnode = dynamic_cast<const CEvaluationNodeNumber*>(node);
644 
645  if (numnode)
646  {
647  //heuristics!
648  if (mUseHeuristics && (numnode->getValue() == 1.0))
649  result.setDimension(0, 0, 0, 0, 0);
650  }
651 
652  return result;
653 }
654 
655 #ifdef COPASI_DEBUG
656 void CFindDimensions::printDebugOutput(const CModel* pModel) const
657 {
658  std::cout << "mDimensions " << mDimensions.size() << std::endl;
659  size_t i, imax = mDimensions.size();
660 
661  for (i = 0; i < imax; ++i)
662  std::cout << i << ": " << mDimensions[i].print(pModel) << std::endl;
663 }
664 #endif // COPASI_DEBUG
C_FLOAT64 mD3
exponent of time base unit
Definition: CDimension.h:78
bool mContradiction
Definition: CDimension.h:85
CDimension compare(const CDimension &rhs) const
Definition: CDimension.cpp:219
const C_FLOAT64 & getValue() const
CDimension mRootDimension
Definition: CDimension.h:165
#define fatalError()
const Type & getType() const
std::string getTimeUnitName() const
Definition: CModel.cpp:2238
void fixDimensionless(bool d1, bool d2, bool d3, bool d4, bool d5)
Definition: CDimension.cpp:263
bool mUnknown
Definition: CDimension.h:84
#define C_INVALID_INDEX
Definition: copasi.h:222
static std::string constructDisplayElement(const std::string &base, C_FLOAT64 exponent)
Definition: CDimension.cpp:61
const std::vector< CDimension > & getDimensions() const
Definition: CDimension.cpp:335
std::string getDisplayString(const CModel *pModel) const
Definition: CDimension.cpp:72
void setUseHeuristics(bool flag)
Definition: CDimension.cpp:330
C_FLOAT64 mD5
exponent of length base unit
Definition: CDimension.h:82
C_FLOAT64 mM2
Definition: CDimension.h:169
CDimension operator-(const CDimension &rhs) const
Definition: CDimension.cpp:191
void setContradiction()
Definition: CDimension.cpp:38
std::string getLengthUnitName() const
Definition: CModel.cpp:2215
void setChemicalEquation(const CChemEq *eq)
Definition: CDimension.cpp:410
void findDimension(size_t index)
find dim for one parameter
Definition: CDimension.cpp:453
std::string print(const CModel *pModel) const
Definition: CDimension.cpp:250
std::vector< CDimension > mDimensions
Definition: CDimension.h:164
C_FLOAT64 mD4
exponent of area base unit
Definition: CDimension.h:80
const CFunction * mpFunction
Definition: CDimension.h:163
long int flag
Definition: f2c.h:52
std::string getAreaUnitName() const
Definition: CModel.cpp:2193
bool operator==(const CDimension &rhs) const
Definition: CDimension.cpp:166
void setUnknown()
Definition: CDimension.cpp:28
std::string getQuantityUnitName() const
Definition: CModel.cpp:2321
C_FLOAT64 mD2
exponent of volume base unit
Definition: CDimension.h:76
#define C_FLOAT64
Definition: copasi.h:92
void setupDimensions()
Definition: CDimension.cpp:295
void findDimensions()
find dim for all parameters
Definition: CDimension.cpp:383
CCopasiNode< Data > * getParent()
Definition: CCopasiNode.h:139
void setDimension(const C_FLOAT64 &d1, const C_FLOAT64 &d2, const C_FLOAT64 &d3, const C_FLOAT64 &d4, const C_FLOAT64 &d5)
Definition: CDimension.cpp:48
The class for handling a chemical kinetic function.
Definition: CFunction.h:29
C_FLOAT64 mD1
exponent of quantity base unit
Definition: CDimension.h:74
void findDimensionsMassAction()
Definition: CDimension.cpp:430
CDimension operator+(const CDimension &rhs) const
Definition: CDimension.cpp:177
Definition: CModel.h:50
size_t getMolecularity(const MetaboliteRole role) const
Definition: CChemEq.cpp:210
std::string getVolumeUnitName() const
Definition: CModel.cpp:2170
bool isUnknown() const
Definition: CDimension.cpp:33
bool isContradiction() const
Definition: CDimension.cpp:43
std::vector< std::string > findDimensionsBoth(const CModel *pModel)
Definition: CDimension.cpp:358
const std::vector< CEvaluationNode * > & getNodeList() const
void setMolecularitiesForMassAction(const size_t &m1, const size_t &m2)
Definition: CDimension.cpp:423
CFunctionParameters & getVariables()
Definition: CFunction.cpp:148
CDimension operator*(const C_FLOAT64 &rhs) const
Definition: CDimension.cpp:205
C_FLOAT64 mM1
Definition: CDimension.h:168