COPASI API  4.16.103
CSBMLunitInterface.cpp
Go to the documentation of this file.
1 // Copyright (C) 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 Sven Sahle and University of Heidelberg
7 // All rights reserved.
8 
9 #include "CSBMLunitInterface.h"
10 #include <sbml/Model.h>
11 
12 #include "copasi/copasi.h"
13 
14 CSBMLunitInterface::CSBMLunitInterface(Model * model, bool unitsFromModel) :
15  mpModel(model),
16  mSBMLLevel(2),
17  mSBMLVersion(4),
18  mAssumeDimensionlessOne(false),
19  mpSBMLTimeUnit(NULL),
20  mpSBMLAmountUnit(NULL),
21  mpSBMLVolumeUnit(NULL),
22  mpSBMLAreaUnit(NULL),
23  mpSBMLLengthUnit(NULL),
24  mpSBMLConflictUnit(NULL)
25 {
26  if (this->mpModel != NULL)
27  {
28  this->mSBMLLevel = this->mpModel->getLevel();
29  this->mSBMLVersion = this->mpModel->getVersion();
30  }
31 
33  initializeFromSBMLModel(unitsFromModel);
34 }
35 
36 /**
37  * Destructor.
38  */
40 {
41  pdelete(this->mpSBMLTimeUnit);
44  pdelete(this->mpSBMLAreaUnit);
47 }
48 
50 {
51  UnitDefinition tmpTime(this->mSBMLLevel, this->mSBMLVersion);
52  Unit tmpUnitTime(this->mSBMLLevel, this->mSBMLVersion);
53  tmpUnitTime.setKind(UNIT_KIND_SECOND);
54  tmpTime.addUnit(&tmpUnitTime);
56 
57  UnitDefinition tmpAmount(this->mSBMLLevel, this->mSBMLVersion);
58  Unit tmpUnitAmount(this->mSBMLLevel, this->mSBMLVersion);
59  tmpUnitAmount.setKind(UNIT_KIND_MOLE);
60  tmpAmount.addUnit(&tmpUnitAmount);
62 
63  UnitDefinition tmpVol(this->mSBMLLevel, this->mSBMLVersion);
64  Unit tmpUnitVol(this->mSBMLLevel, this->mSBMLVersion);
65  tmpUnitVol.setKind(UNIT_KIND_LITRE);
66  tmpVol.addUnit(&tmpUnitVol);
68 
69  UnitDefinition tmpAr(this->mSBMLLevel, this->mSBMLVersion);
70  Unit tmpUnitAr(this->mSBMLLevel, this->mSBMLVersion);
71  tmpUnitAr.setKind(UNIT_KIND_METRE);
72  tmpUnitAr.setExponent(2);
73  tmpAr.addUnit(&tmpUnitAr);
75 
76  UnitDefinition tmpL(this->mSBMLLevel, this->mSBMLVersion);
77  Unit tmpUnitL(this->mSBMLLevel, this->mSBMLVersion);
78  tmpUnitL.setKind(UNIT_KIND_METRE);
79  tmpL.addUnit(&tmpUnitL);
81 
83 }
84 
86 {
87  if (!this->mpModel) return;
88 
89  //initialize global units
90  if (unitsFromModel)
91  {
92  if (this->mpModel->getUnitDefinition("time"))
93  (*mpSBMLTimeUnit) = CSBMLunitInformation(this->mpModel->getUnitDefinition("time"), CSBMLunitInformation::GLOBAL);
94 
95  if (this->mpModel->getUnitDefinition("substance"))
96  (*mpSBMLAmountUnit) = CSBMLunitInformation(this->mpModel->getUnitDefinition("substance"), CSBMLunitInformation::GLOBAL);
97 
98  if (this->mpModel->getUnitDefinition("volume"))
99  (*mpSBMLVolumeUnit) = CSBMLunitInformation(this->mpModel->getUnitDefinition("volume"), CSBMLunitInformation::GLOBAL);
100 
101  if (this->mpModel->getUnitDefinition("area"))
102  (*mpSBMLAreaUnit) = CSBMLunitInformation(this->mpModel->getUnitDefinition("area"), CSBMLunitInformation::GLOBAL);
103 
104  if (this->mpModel->getUnitDefinition("length"))
105  (*mpSBMLLengthUnit) = CSBMLunitInformation(this->mpModel->getUnitDefinition("length"), CSBMLunitInformation::GLOBAL);
106  }
107 
108  mSBMLObjectsMap.clear();
109 
110  //initialize list of species
111  unsigned int i;
112 
113  for (i = 0; i < this->mpModel->getNumSpecies(); i++)
114  {
115  Species * s = this->mpModel->getSpecies(i);
116 
117  if (unitsFromModel && s->isSetUnits())
118  {
119  std::map<std::string, CSBMLunitInformation>::iterator pos = this->mSBMLObjectsMap.find(s->getId());
120 
121  if (pos != this->mSBMLObjectsMap.end())
122  {
123  pos->second = CSBMLunitInformation(s->getDerivedUnitDefinition(), CSBMLunitInformation::PROVIDED);
124  }
125  else
126  {
127  mSBMLObjectsMap.insert(std::pair<std::string, CSBMLunitInformation>(s->getId() , CSBMLunitInformation(s->getDerivedUnitDefinition(), CSBMLunitInformation::PROVIDED)));
128  }
129  }
130  else //take info from global units
131  {
132  if (s->getHasOnlySubstanceUnits())
133  {
134  std::map<std::string, CSBMLunitInformation>::iterator pos = this->mSBMLObjectsMap.find(s->getId());
135 
136  if (pos != this->mSBMLObjectsMap.end())
137  {
138  pos->second = (*mpSBMLAmountUnit);
139  }
140  else
141  {
142  mSBMLObjectsMap.insert(std::pair<std::string, CSBMLunitInformation>(s->getId(), (*mpSBMLAmountUnit)));
143  }
144  }
145  else //concentration
146  {
147  CSBMLunit tmp(this->mpModel->getLevel(), this->mpModel->getVersion());
148  Compartment* comp = this->mpModel->getCompartment(s->getCompartment());
149 
150  if (!comp)
151  {
152  std::map<std::string, CSBMLunitInformation>::iterator pos = this->mSBMLObjectsMap.find(s->getId());
153 
154  if (pos != this->mSBMLObjectsMap.end())
155  {
156  pos->second = (*mpSBMLConflictUnit);
157  }
158  else
159  {
160  mSBMLObjectsMap.insert(std::pair<std::string, CSBMLunitInformation>(s->getId(), (*mpSBMLConflictUnit)));
161  }
162 
163  continue;
164  }
165 
166  switch (comp->getSpatialDimensions())
167  {
168  case 0:
169  break;
170  case 1:
171  tmp = (*mpSBMLLengthUnit);
172  break;
173  case 2:
174  tmp = (*mpSBMLAreaUnit);
175  break;
176  case 3:
177  tmp = (*mpSBMLVolumeUnit);
178  break;
179  default:
180  break;
181  };
182 
183  tmp.invert();
184 
185  tmp.multiply(*mpSBMLAmountUnit);
186 
187  std::map<std::string, CSBMLunitInformation>::iterator pos = this->mSBMLObjectsMap.find(s->getId());
188 
189  if (pos != this->mSBMLObjectsMap.end())
190  {
192  }
193  else
194  {
195  mSBMLObjectsMap.insert(std::pair<std::string, CSBMLunitInformation>(s->getId(), CSBMLunitInformation(tmp, CSBMLunitInformation::GLOBAL)));
196  }
197 
198  //TODO: it should be DEFAULT rather than GLOBAL if both amount unit and size unit are DEFAULT.
199  }
200  }
201  }
202 
203  //initialize list of compartments
204  for (i = 0; i < this->mpModel->getNumCompartments(); i++)
205  {
206  Compartment *c = this->mpModel->getCompartment(i);
207 
208  if (unitsFromModel && c->isSetUnits())
209  {
210 
211  std::map<std::string, CSBMLunitInformation>::iterator pos = this->mSBMLObjectsMap.find(c->getId());
212 
213  if (pos != this->mSBMLObjectsMap.end())
214  {
215  pos->second = CSBMLunitInformation(c->getDerivedUnitDefinition(), CSBMLunitInformation::PROVIDED);
216  }
217  else
218  {
219  mSBMLObjectsMap.insert(std::pair<std::string, CSBMLunitInformation>(c->getId(), CSBMLunitInformation(c->getDerivedUnitDefinition(), CSBMLunitInformation::PROVIDED)));
220  }
221  }
222  else //take info from global units
223  {
224  CSBMLunitInformation tmp(this->mSBMLLevel, this->mSBMLVersion);
225 
226  switch (c->getSpatialDimensions())
227  {
228  case 0:
229  break;
230  case 1:
231  tmp = (*mpSBMLLengthUnit);
232  break;
233  case 2:
234  tmp = (*mpSBMLAreaUnit);
235  break;
236  case 3:
237  tmp = (*mpSBMLVolumeUnit);
238  break;
239  default:
240  break;
241  };
242 
243  std::map<std::string, CSBMLunitInformation>::iterator pos = this->mSBMLObjectsMap.find(c->getId());
244 
245  if (pos != this->mSBMLObjectsMap.end())
246  {
247  pos->second = tmp;
248  }
249  else
250  {
251  mSBMLObjectsMap.insert(std::pair<std::string, CSBMLunitInformation>(c->getId(), tmp));
252  }
253  }
254  }
255 
256  //initialize list of global parameters
257  for (i = 0; i < this->mpModel->getNumParameters(); i++)
258  {
259  Parameter *p = this->mpModel->getParameter(i);
260 
261  if (unitsFromModel && p->isSetUnits())
262  {
263  std::map<std::string, CSBMLunitInformation>::iterator pos = this->mSBMLObjectsMap.find(p->getId());
264 
265  if (pos != this->mSBMLObjectsMap.end())
266  {
267  pos->second = CSBMLunitInformation(p->getDerivedUnitDefinition(), CSBMLunitInformation::PROVIDED);
268  }
269  else
270  {
271  mSBMLObjectsMap.insert(std::pair<std::string, CSBMLunitInformation>(p->getId(), CSBMLunitInformation(p->getDerivedUnitDefinition(), CSBMLunitInformation::PROVIDED)));
272  }
273  }
274  else
275  {
276  std::map<std::string, CSBMLunitInformation>::iterator pos = this->mSBMLObjectsMap.find(p->getId());
277 
278  if (pos != this->mSBMLObjectsMap.end())
279  {
281  }
282  else
283  {
284  mSBMLObjectsMap.insert(std::pair<std::string, CSBMLunitInformation>(p->getId(), CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion, CSBMLunitInformation::UNKNOWN)));
285  }
286  }
287  }
288 
289  //initialize list of local parameters
290  mSBMLLocalParametersMap.clear();
291  unsigned int j;
292 
293  for (j = 0; j < this->mpModel->getNumReactions(); j++)
294  {
295  Reaction * reaction = this->mpModel->getReaction(j);
296  std::map<std::string, CSBMLunitInformation> tmpMap;
297 
298  if (reaction->getKineticLaw() != NULL)
299  {
300  for (i = 0; i < reaction->getKineticLaw()->getNumParameters(); i++)
301  {
302  Parameter *p = reaction->getKineticLaw()->getParameter(i);
303 
304  if (unitsFromModel && p->isSetUnits())
305  {
306  //tmpMap[p->getId()]=CSBMLunitInformation(p->getDerivedUnitDefinition(), CSBMLunitInformation::PROVIDED);
307  UnitDefinition* tmpUD = this->mpModel->getUnitDefinition(p->getUnits());
308 
309  if (tmpUD)
310  {
311  std::map<std::string, CSBMLunitInformation>::iterator pos = tmpMap.find(p->getId());
312 
313  if (pos != tmpMap.end())
314  {
316  }
317  else
318  {
319  tmpMap.insert(std::pair<std::string, CSBMLunitInformation>(p->getId(), CSBMLunitInformation(tmpUD, CSBMLunitInformation::PROVIDED)));
320  }
321  }
322  else
323  {
324  //this is just a workaround
325  if (p->getUnits() == "dimensionless")
326  {
327  std::map<std::string, CSBMLunitInformation>::iterator pos = tmpMap.find(p->getId());
328 
329  if (pos != tmpMap.end())
330  {
332  }
333  else
334  {
335  tmpMap.insert(std::pair<std::string, CSBMLunitInformation>(p->getId(), CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion, CSBMLunitInformation::PROVIDED)));
336  }
337  }
338  else
339  {
340  std::map<std::string, CSBMLunitInformation>::iterator pos = tmpMap.find(p->getId());
341 
342  if (pos != tmpMap.end())
343  {
345  }
346  else
347  {
348  tmpMap.insert(std::pair<std::string, CSBMLunitInformation>(p->getId(), CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion, CSBMLunitInformation::UNKNOWN)));
349  }
350 
351  // "Could not resolve unit ID " << p->getUnits() << " for " << p->getId() << std::endl;
352  //TODO
353  }
354  }
355  }
356  else //take info from global units
357  {
358  std::map<std::string, CSBMLunitInformation>::iterator pos = tmpMap.find(p->getId());
359 
360  if (pos != tmpMap.end())
361  {
363  }
364  else
365  {
366  tmpMap.insert(std::pair<std::string, CSBMLunitInformation>(p->getId(), CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion, CSBMLunitInformation::UNKNOWN)));
367  }
368  }
369  }
370  }
371 
372  mSBMLLocalParametersMap[reaction->getId()] = tmpMap;
373  }
374 
375  //determine units of stoichiometries
376 
377  //**** generate list of all mathematical expressions that need to be considered
378 
379  mSBMLExpressions.clear();
380 
381  //stoichiometry math TODO
382 
383  //kinetic laws
384  //construct the units for amount/time (for kinetic laws)
385  CSBMLunitInformation amountPerTimeUnit = (*mpSBMLTimeUnit);
386  amountPerTimeUnit.invert();
387  amountPerTimeUnit.multiply(*mpSBMLAmountUnit);
388 
390  amountPerTimeUnit.setInfo(CSBMLunitInformation::GLOBAL);
391 
392  for (i = 0; i < this->mpModel->getNumReactions(); i++)
393  {
394  Reaction * reaction = this->mpModel->getReaction(i);
395 
396  if (!reaction->getKineticLaw()) continue;
397 
398  if (!reaction->getKineticLaw()->isSetMath()) continue;
399 
401  tmp.mTypeDescription = "Kinetic Law";
402  tmp.mObjectDisplayString = reaction->isSetName() ? reaction->getName() : reaction->getId();
403  tmp.mpExpression = reaction->getKineticLaw()->getMath();
404  tmp.mRootUnit = amountPerTimeUnit;
405  tmp.mReactionId = reaction->getId();
406  mSBMLExpressions.push_back(tmp);
407  }
408 
409  //rate rules
410  //assignment rules
411  //algebraic rules
412  for (i = 0; i < this->mpModel->getNumRules(); i++)
413  {
414  Rule* rule = this->mpModel->getRule(i);
415 
416  if (!rule->isSetMath()) continue;
417 
419  tmp.mpExpression = rule->getMath();
420 
421  if (rule->isRate())
422  {
423  tmp.mTypeDescription = "Rate rule";
424  tmp.mObjectDisplayString = rule->getVariable(); //TODO
425  tmp.mRootObject = rule->getVariable();
426  tmp.mPerTime = true;
427  }
428  else if (rule->isAssignment())
429  {
430  tmp.mTypeDescription = "Assignment rule";
431  tmp.mObjectDisplayString = rule->getVariable(); //TODO
432  tmp.mRootObject = rule->getVariable();
433  tmp.mPerTime = false;
434  }
435  else if (rule->isAlgebraic())
436  {
437  tmp.mTypeDescription = "Algebraic rule";
438  //nothing to be done. UNKNOWN is the default for a CSBMLunitInformation
439  //tmp.mRootUnit=UNKNOWN;
440  }
441 
442  mSBMLExpressions.push_back(tmp);
443  }
444 
445  //Events
446  for (i = 0; i < this->mpModel->getNumEvents(); i++)
447  {
448  Event* event = this->mpModel->getEvent(i);
449 
450  //trigger
451  if (event->isSetTrigger() && event->getTrigger()->isSetMath())
452  {
454  tmp.mTypeDescription = "Event trigger";
455  tmp.mObjectDisplayString = event->isSetName() ? event->getName() : event->getId();
456  tmp.mpExpression = event->getTrigger()->getMath();
457  mSBMLExpressions.push_back(tmp);
458  }
459 
460  //delay
461  if (event->isSetDelay() && event->getDelay()->isSetMath())
462  {
464  tmp.mTypeDescription = "Event delay";
465  tmp.mObjectDisplayString = event->isSetName() ? event->getName() : event->getId();
466  tmp.mpExpression = event->getDelay()->getMath();
467  tmp.mRootUnit = (*mpSBMLTimeUnit);
468  mSBMLExpressions.push_back(tmp);
469  }
470 
471  //event assignments
472  for (j = 0; j < event->getNumEventAssignments(); ++j)
473  {
474  EventAssignment* ea = event->getEventAssignment(j);
475 
476  if (ea->isSetMath())
477  {
479  tmp.mTypeDescription = "Event assignment";
480  tmp.mObjectDisplayString = ea->getVariable(); //TODO
481  tmp.mpExpression = ea->getMath();
482  tmp.mRootObject = ea->getVariable();
483  tmp.mPerTime = false;
484  mSBMLExpressions.push_back(tmp);
485  }
486  }
487  }
488 
489  //initial assignments
490 
491  //constraints?
492 
493  mSBMLNumbersMap.clear();
495 }
496 
498 {
499  //**** global objects *****
500 
501  //species
502  unsigned int i;
503 
504  //do nothing for level 2 since the species units are never unknown
505  // for (i = 0; i < mpModel->getNumSpecies(); i++)
506  // {
507  // Species * s = mpModel->getSpecies(i);
508  //}
509 
510  //compartments
511  //do nothing for level 2 since the compartment units are never unknown
512 
513  //global parameters
514  for (i = 0; i < mpModel->getNumParameters(); i++)
515  {
516  Parameter *p = mpModel->getParameter(i);
518 
519  //if the unit could be derived and it does not contain a symbolic exponent
520  if (tmp && tmp->getInfo() == CSBMLunitInformation::DERIVED && tmp->getSymbolicExpExp() == 0)
521  {
522  //first try to find an equivalent unit in the model
523  unsigned int j;
524 
525  for (j = 0; j < mpModel->getNumUnitDefinitions(); ++j)
526  if (UnitDefinition::areEquivalent(&tmp->getSBMLUnitDefinition(), mpModel->getUnitDefinition(j)))
527  break;
528 
529  if (j < mpModel->getNumUnitDefinitions())
530  p->setUnits(mpModel->getUnitDefinition(j)->getId());
531  else
532  {
533  //we have to create a new unit definition in the model
534  std::string tmpstring;
535  unsigned int id = 0;
536 
537  do
538  {
539  std::ostringstream tmpid; tmpid << "unit_" << id;
540  tmpstring = tmpid.str();
541  ++id;
542  }
543  while (mpModel->getUnitDefinition(tmpstring));
544 
545  tmp->getSBMLUnitDefinition().setId(tmpstring);
546  tmp->getSBMLUnitDefinition().unsetName();
547  tmp->getSBMLUnitDefinition().unsetMetaId();
548  tmp->getSBMLUnitDefinition().unsetNotes();
549  tmp->getSBMLUnitDefinition().unsetAnnotation();
550  tmp->getSBMLUnitDefinition().unsetSBOTerm();
551  mpModel->addUnitDefinition(&tmp->getSBMLUnitDefinition());
552  p->setUnits(tmpstring);
553  }
554  }
555  }
556 
557  //local parameters
558  unsigned int j;
559 
560  for (j = 0; j < mpModel->getNumReactions(); j++)
561  {
562  Reaction * reaction = mpModel->getReaction(j);
563 
564  if (reaction->getKineticLaw() != NULL)
565  for (i = 0; i < reaction->getKineticLaw()->getNumParameters(); i++)
566  {
567  Parameter *p = reaction->getKineticLaw()->getParameter(i);
568 
569  CSBMLunitInformation * tmp = getMappedUnitFromIdentifier(p->getId(), CEnvironmentInformation(reaction->getId()));
570  //TODO this could be easier directly from the map instead of using getMappedUnitFromIdentifier()
571 
572  //if the unit could be derived and it does not contain a symbolic exponent
573  if (tmp && tmp->getInfo() == CSBMLunitInformation::DERIVED && tmp->getSymbolicExpExp() == 0)
574  {
575  //first try to find an equivalent unit in the model
576  unsigned int j;
577 
578  for (j = 0; j < mpModel->getNumUnitDefinitions(); ++j)
579  if (UnitDefinition::areEquivalent(&tmp->getSBMLUnitDefinition(), mpModel->getUnitDefinition(j)))
580  break;
581 
582  if (j < mpModel->getNumUnitDefinitions())
583  p->setUnits(mpModel->getUnitDefinition(j)->getId());
584  else
585  {
586  //we have to create a new unit definition in the model
587  std::string tmpstring;
588  unsigned int id = 0;
589 
590  do
591  {
592  std::ostringstream tmpid; tmpid << "unit_" << id;
593  tmpstring = tmpid.str();
594  ++id;
595  }
596  while (mpModel->getUnitDefinition(tmpstring));
597 
598  tmp->getSBMLUnitDefinition().setId(tmpstring);
599  tmp->getSBMLUnitDefinition().unsetName();
600  tmp->getSBMLUnitDefinition().unsetMetaId();
601  tmp->getSBMLUnitDefinition().unsetNotes();
602  tmp->getSBMLUnitDefinition().unsetAnnotation();
603  tmp->getSBMLUnitDefinition().unsetSBOTerm();
604  mpModel->addUnitDefinition(&tmp->getSBMLUnitDefinition());
605  p->setUnits(tmpstring);
606  }
607  }
608  }
609  }
610 
611  //numbers TODO (not possible at the moment)
612 }
613 
615 {
616  //std::vector<unsigned int> frequency_global, frequency_local, frequency_numbers, frequency;
618 
619  //CONFLICT is used as index for a unit conflict.
620  const unsigned int CONFLICT = 5;
621  mStatistics.global.resize(CONFLICT + 1);
622  mStatistics.local.resize(CONFLICT + 1);
623  mStatistics.numbers.resize(CONFLICT + 1);
624  mStatistics.all.resize(CONFLICT + 1);
625 
626  //collect statistics from global objects with units
627  std::map<std::string, CSBMLunitInformation>::const_iterator it, itEnd = mSBMLObjectsMap.end();
628 
629  for (it = mSBMLObjectsMap.begin(); it != itEnd; ++it)
630  {
631  if (it->second.isConflict())
632  ++mStatistics.global[CONFLICT];
633  else
634  ++mStatistics.global[it->second.getInfo()];
635  }
636 
637  //collect statistics from localparameter
638  std::map<std::string, std::map<std::string, CSBMLunitInformation> >::const_iterator rit;
639 
640  for (rit = mSBMLLocalParametersMap.begin(); rit != mSBMLLocalParametersMap.end(); ++rit)
641  {
642  for (it = rit->second.begin(); it != rit->second.end(); ++it)
643  {
644  if (it->second.isConflict())
645  ++mStatistics.local[CONFLICT];
646  else
647  ++mStatistics.local[it->second.getInfo()];
648  }
649  }
650 
651  //collect statistics from numbers
652  std::map<const ASTNode*, CSBMLunitInformation>::const_iterator nit, nitEnd = mSBMLNumbersMap.end();
653 
654  for (nit = mSBMLNumbersMap.begin(); nit != nitEnd; ++nit)
655  {
656  if (nit->second.isConflict())
657  ++mStatistics.numbers[CONFLICT];
658  else
659  ++mStatistics.numbers[nit->second.getInfo()];
660  }
661 
662  unsigned int i;
663 
664  for (i = 0; i < mStatistics.all.size(); ++i)
666 }
667 
669 {
670  if (stat.all.size() != 6) return;
671 
672  std::cout << "Global: " << " ? " << stat.global[0]
673  << " default " << stat.global[1]
674  << " glob " << stat.global[2]
675  << " pro " << stat.global[3]
676  << " deriv " << stat.global[4]
677  << (stat.global[5] > 0 ? " conflict " : " ") << stat.global[5] << std::endl;
678  std::cout << "Local: " << " ? " << stat.local[0]
679  << " default " << stat.local[1]
680  << " glob " << stat.local[2]
681  << " pro " << stat.local[3]
682  << " deriv " << stat.local[4]
683  << (stat.local[5] > 0 ? " conflict " : " ") << stat.local[5] << std::endl;
684  std::cout << "Numbers:" << " ? " << stat.numbers[0]
685  << " default " << stat.numbers[1]
686  << " glob " << stat.numbers[2]
687  << " pro " << stat.numbers[3]
688  << " deriv " << stat.numbers[4]
689  << (stat.numbers[5] > 0 ? " conflict " : " ") << stat.numbers[5] << std::endl;
690 
691  if (flag) std::cout << "_";
692 
693  std::cout << "Sum :" << " ? " << stat.all[0]
694  << " default " << stat.all[1]
695  << " glob " << stat.all[2]
696  << " pro " << stat.all[3]
697  << " deriv " << stat.all[4]
698  << (stat.all[5] > 0 ? " CONFLICT_SUM " : " ") << stat.all[5] << std::endl;
699 
700  if (stat.all[0] == 0)
701  std::cout << "***SUCCESS***" << std::endl;
702  else if (stat.global[0] + stat.local[0] == 0)
703  std::cout << "only some numbers are still unknown" << std::endl;
704  else
705  std::cout << "some unknown units left..." << std::endl;
706 
707  //if (stat.all[5]>0)
708  // debugOutput();
709 }
710 
711 std::vector<std::string> CSBMLunitInterface::getListOfObjectsWithGivenUnitStatus(int status) const
712 {
713  std::vector<std::string> ret;
714 
715  std::map<std::string, CSBMLunitInformation>::const_iterator it, itEnd = mSBMLObjectsMap.end();
716 
717  for (it = mSBMLObjectsMap.begin(); it != itEnd; ++it)
718  {
719  if (status == 5 && it->second.isConflict())
720  ret.push_back(it->first);
721 
722  if (status < 5 && (int)it->second.getInfo() == status)
723  ret.push_back(it->first);
724  }
725 
726  return ret;
727 }
728 
729 std::vector<std::pair<std::string, std::string> > CSBMLunitInterface::getListOfLocalParametersWithGivenUnitStatus(int status) const
730 {
731  std::vector<std::pair<std::string, std::string> > ret;
732 
733  std::map<std::string, CSBMLunitInformation>::const_iterator it;//, itEnd = mSBMLObjectsMap.end();
734  std::map<std::string, std::map<std::string, CSBMLunitInformation> >::const_iterator rit;
735 
736  for (rit = mSBMLLocalParametersMap.begin(); rit != mSBMLLocalParametersMap.end(); ++rit)
737  {
738  for (it = rit->second.begin(); it != rit->second.end(); ++it)
739  {
740  if (status == 5 && it->second.isConflict())
741  ret.push_back(std::pair<std::string, std::string>(rit->first, it->first));
742 
743  if (status < 5 && (int)it->second.getInfo() == status)
744  ret.push_back(std::pair<std::string, std::string>(rit->first, it->first));
745  }
746  }
747 
748  return ret;
749 }
750 
752 {
753  std::string ret;
754 
755  ret += "Some objects in the SBML file have unknown units. \n";
756  ret += "Global objects:\n";
757 
758  std::vector<std::string> globalobjects = getListOfObjectsWithGivenUnitStatus(0); //unknown units
759  unsigned int i;
760 
761  for (i = 0; i < globalobjects.size(); ++i)
762  {
763  const Compartment* c = mpModel->getCompartment(globalobjects[i]);
764 
765  if (c)
766  ret += " Compartment " + c->getId() + " " + c->getName() + "\n";
767 
768  const Species* s = mpModel->getSpecies(globalobjects[i]);
769 
770  if (s)
771  ret += " Species " + s->getId() + " " + s->getName() + "\n";
772 
773  const Parameter* p = mpModel->getParameter(globalobjects[i]);
774 
775  if (p)
776  ret += " Global parameter " + p->getId() + " " + p->getName() + "\n";
777  }
778 
779  ret += "\n";
780  ret += "Local parameters:\n";
781 
782  std::vector<std::pair<std::string, std::string> > localparameters =
784 
785  for (i = 0; i < localparameters.size(); ++i)
786  {
787  const Reaction* r = mpModel->getReaction(localparameters[i].first);
788 
789  if (r)
790  {
791  ret += " Reaction " + r->getId() + " " + r->getName()
792  + ", parameter " + localparameters[i].second + "\n";
793  }
794  }
795 
796  //TODO this is not really very nice output. Needs improvement
797  return ret;
798 }
799 
801 {
802  std::cout << "global units:" << std::endl;
803  std::cout << "Time: " << mpSBMLTimeUnit->getDisplayString() << std::endl;
804  std::cout << "Amount: " << mpSBMLAmountUnit->getDisplayString() << std::endl;
805  std::cout << "Volume: " << mpSBMLVolumeUnit->getDisplayString() << std::endl;
806  std::cout << "Area: " << mpSBMLAreaUnit->getDisplayString() << std::endl;
807  std::cout << "Length: " << mpSBMLLengthUnit->getDisplayString() << std::endl;
808 
809  std::cout << std::endl;
810 
811  std::map<std::string, CSBMLunitInformation>::const_iterator it, itEnd = mSBMLObjectsMap.end();
812 
813  for (it = mSBMLObjectsMap.begin(); it != itEnd; ++it)
814  {
815  std::cout << it->first << " " << it->second.getDisplayString() << std::endl;
816  }
817 
818  std::cout << std::endl;
819 
820  std::map<std::string, std::map<std::string, CSBMLunitInformation> >::const_iterator rit;
821 
822  for (rit = mSBMLLocalParametersMap.begin(); rit != mSBMLLocalParametersMap.end(); ++rit)
823  {
824  for (it = rit->second.begin(); it != rit->second.end(); ++it)
825  {
826  std::cout << rit-> first << " "
827  << it->first << " " << it->second.getDisplayString() << std::endl;
828  }
829  }
830 
831  std::cout << std::endl;
832  std::cout << "Expressions: \n\n";
833  unsigned int i;
834 
835  for (i = 0; i < mSBMLExpressions.size(); ++i)
836  {
838  std::cout << std::string(SBML_formulaToString(tmp.mpExpression)) << " :: "
839  << (tmp.mPerTime ? "per time" : "") << " "
840  << tmp.mRootObject << " ";
841 
842  if (tmp.mReactionId != "")
843  std::cout << "env:" << tmp.mReactionId << " ";
844 
845  if (tmp.mRootObject == "")
846  std::cout << tmp.mRootUnit.getDisplayString();
847 
848  if (tmp.mErrorCode) std::cout << "Conflict: " << tmp.mErrorCode;
849 
850  std::cout << std::endl;
851  }
852 }
853 
854 const std::set<const ASTNode *> & CSBMLunitInterface::getListOfConflictingNodes() const
855 {
856  return mConflictingNodes;
857 }
858 
860 {
861  if (!mpModel) return;
862 
863  //clear set of conflicting nodes
864  mConflictingNodes.clear();
865 
866  //first a very crude implementation... TODO
867  std::vector<CExpressionInformation>::iterator it, itEnd = mSBMLExpressions.end();
868  unsigned int count;
869 
870  for (count = 0; count < 3; ++count)
871  {
872  for (it = mSBMLExpressions.begin(); it != itEnd; ++it)
873  {
874  //debug
875  if (false)
876  {
877  const CExpressionInformation& tmp = *it;
878  std::cout << std::string(SBML_formulaToString(tmp.mpExpression)) << " :: "
879  << (tmp.mPerTime ? "per time" : "") << " "
880  << tmp.mRootObject << " ";
881 
882  if (tmp.mReactionId != "")
883  std::cout << "env:" << tmp.mReactionId << " ";
884 
885  if (tmp.mRootObject == "")
886  std::cout << tmp.mRootUnit.getDisplayString();
887 
888  std::cout << std::endl;
889  }
890 
891  handleOneExpression(*it);
892  }
893  }
894 
896 }
897 
899 {
900  mError = 0;
901  CEnvironmentInformation environment;
902 
903  if (ei.mRootObject == "")
904  {
905  //set the environment (for local parameters)
906  environment.mReactionID = ei.mReactionId;
907 
908  //call recursive evaluation
909  recursion(ei.mpExpression, ei.mRootUnit, environment);
910  }
911  else
912  {
913  CSBMLunitInformation * pNodeUnit = getMappedUnitFromIdentifier(ei.mRootObject, environment);
914 
915  if (!pNodeUnit) return;
916 
917  CSBMLunitInformation sourceUnit = *pNodeUnit;
918 
919  if (ei.mPerTime && sourceUnit.getInfo() > CSBMLunitInformation::UNKNOWN)
920  {
921  //devide unit of rule target by time
922  CSBMLunitInformation invTime = (*mpSBMLTimeUnit);
923  invTime.invert();
924  sourceUnit.multiply(invTime);
925  }
926 
927  CSBMLunitInformation tmp = recursion(ei.mpExpression, sourceUnit, environment);
928 
930  {
931  //multiply unit of the rate rule expression by time
932  tmp.multiply(*mpSBMLTimeUnit);
933  }
934 
935  handleTerminalNode(tmp, pNodeUnit, NULL);
936  }
937 
938  //TODO handle mpRootObject case
939 
940  if (mError > ei.mErrorCode) ei.mErrorCode = mError;
941 }
942 
943 //*********************************************************************************
944 
945 //#define UNIT_DEBUG
946 
948  const CSBMLunitInformation & ui,
949  const CEnvironmentInformation & ei)
950 {
951 
952  CSBMLunitInformation ret(this->mSBMLLevel, this->mSBMLVersion);
953 
954  if (!node) return ret;
955 
956  //handle delay first since it is unclear from the documentation in which class of nodes it is included
957  if (node->getType() == AST_FUNCTION_DELAY)
958  {
959  if (node->getNumChildren() != 2)
960  return ret;
961 
962  CSBMLunitInformation tmpTimeUnit = (*mpSBMLTimeUnit);
963  recursion(node->getChild(1), tmpTimeUnit, ei); //second child should have time units
964  return recursion(node->getChild(0), ui, ei); //first child has the same units as the parent
965  }
966 
967  //object node
968  if (isObject(node))
969  {
970  //time nodes
971  if (node->getType() == AST_NAME_TIME)
972  {
973  //TODO: does not work???
974  CSBMLunitInformation tmpTimeUnit = (*mpSBMLTimeUnit);
975  return handleTerminalNode(ui, &tmpTimeUnit, node);
976  }
977 
978  //check if it is a reaction ID
979  //TODO
980 
981  //check if the object is a function variable
982  ASTNode* variable = resolveVariableName(getIdentifier(node), ei);
983 
984  if (variable)
985  {
986  //drop one level from the environment stack
987  CEnvironmentInformation tmpIE = ei;
988  tmpIE.drop();
989  return recursion(variable, ui, tmpIE);
990  }
991 
992  //now it should be certain that the object is a reference to an sbml object with a mapped unit.
994 
995  if (!pNodeUnit) return (*mpSBMLConflictUnit);
996 
997  return handleTerminalNode(ui, pNodeUnit, node);
998  }//is object node
999 
1000  //number node
1001  else if (isNumber(node))
1002  {
1003  CSBMLunitInformation* pNodeUnit = NULL;
1004  std::map<const ASTNode*, CSBMLunitInformation>::iterator pos = mSBMLNumbersMap.find(node);
1005 
1006  if (pos != mSBMLNumbersMap.end())
1007  {
1008  pNodeUnit = &pos->second;
1009  }
1010  else
1011  {
1012  pNodeUnit = &((mSBMLNumbersMap.insert(
1013  std::pair<const ASTNode*, CSBMLunitInformation>
1014  (node, CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion))
1015  )).first->second);
1016  }
1017 
1018 // //check if the node is already in the map
1019 // std::map<const ASTNode*, CSBMLunitInformation>::iterator it = mSBMLNumbersMap.find(node);
1020 // if (it != mSBMLNumbersMap.end())
1021 // pNodeUnit=&it->second;
1022 // else //create new entry in the map
1023 // {
1024 // pNodeUnit=&mSBMLNumbersMap[node];
1025 //}
1026 
1028  {
1029  if (getValueFromNumberNode(node) == 1.0 || getValueFromNumberNode(node) == -1.0)
1031  }
1032 
1033  return handleTerminalNode(ui, pNodeUnit, node);
1034  }
1035 
1036  //operator node
1037  else if (isOperator(node))
1038  {
1039  switch (node->getType())
1040  {
1041  case AST_PLUS:
1042  case AST_MINUS:
1043  case AST_RELATIONAL_EQ:
1044  case AST_RELATIONAL_GEQ:
1045  case AST_RELATIONAL_GT:
1046  case AST_RELATIONAL_LEQ:
1047  case AST_RELATIONAL_LT:
1048  case AST_RELATIONAL_NEQ:
1049  return recursionEqual(node, ui, ei);
1050  break;
1051 
1052  case AST_TIMES:
1053  return recursionTimes(node, ui, ei);
1054  break;
1055 
1056  case AST_DIVIDE:
1057  return recursionDivide(node, ui, ei);
1058  break;
1059 
1060  case AST_POWER:
1061  case AST_FUNCTION_POWER:
1062  return recursionPower(node, ui, ei);
1063  break;
1064 
1065  default:
1066  break;
1067  }
1068  }
1069 
1070  else if (isFunctionCall(node))
1071  {
1072  // function call
1073  //check if the object is a function call
1074  FunctionDefinition *fd = resolveFunctionName(getIdentifier(node));
1075 
1076  if (fd)
1077  {
1078  // could resolve function call
1079  unsigned int i, numArgs = fd->getNumArguments();
1080  assert(numArgs == node->getNumChildren());
1081 
1082  //create mapping from variable ID to the node that should be substituted
1083  std::map<std::string, ASTNode*> tmpMap;
1084 
1085  for (i = 0; i < numArgs; ++i)
1086  tmpMap[fd->getArgument(i)->getName()] = node->getChild(i);
1087 
1088  CEnvironmentInformation tmpEI = ei;
1089  tmpEI.push(tmpMap);
1090  return recursion(fd->getBody(), ui, tmpEI);
1091  }
1092  }
1093 
1094  else if (isBuiltInFunctionCall(node))
1095  {
1096  // built in function
1097  switch (node->getType())
1098  {
1099  case AST_FUNCTION_ARCCOS:
1100  case AST_FUNCTION_ARCCOSH:
1101  case AST_FUNCTION_ARCCOT:
1102  case AST_FUNCTION_ARCCOTH:
1103  case AST_FUNCTION_ARCCSC:
1104  case AST_FUNCTION_ARCCSCH:
1105  case AST_FUNCTION_ARCSEC:
1106  case AST_FUNCTION_ARCSECH:
1107  case AST_FUNCTION_ARCSIN:
1108  case AST_FUNCTION_ARCSINH:
1109  case AST_FUNCTION_ARCTAN:
1110  case AST_FUNCTION_ARCTANH:
1111 
1112  case AST_FUNCTION_COS:
1113  case AST_FUNCTION_COSH:
1114  case AST_FUNCTION_COT:
1115  case AST_FUNCTION_COTH:
1116  case AST_FUNCTION_CSC:
1117  case AST_FUNCTION_CSCH:
1118  case AST_FUNCTION_SEC:
1119  case AST_FUNCTION_SECH:
1120  case AST_FUNCTION_SIN:
1121  case AST_FUNCTION_SINH:
1122  case AST_FUNCTION_TAN:
1123  case AST_FUNCTION_TANH:
1124 
1125  case AST_FUNCTION_EXP:
1126  case AST_FUNCTION_FACTORIAL:
1127  case AST_FUNCTION_LN:
1128  case AST_FUNCTION_LOG:
1131  break;
1132 
1133  case AST_FUNCTION_CEILING:
1134  case AST_FUNCTION_FLOOR:
1135  case AST_FUNCTION_ABS:
1136  return recursion(node->getChild(0), ui, ei);
1137  break;
1138 
1139  default:
1140  break;
1141  }
1142  }
1143 
1144  switch (node->getType())
1145  {
1146  case AST_CONSTANT_E:
1147  case AST_CONSTANT_PI:
1149  break;
1150 
1151  case AST_FUNCTION_PIECEWISE:
1152  return recursionPiecewise(node, ui, ei);
1153 
1154  default:
1155  break;
1156  }
1157 
1158  //TODO numbers, functions, calls, variables, choice, delay, time, ...
1159  //logicals (comparison operators), constants (pi, ...), ...
1160 
1161  //fallback: just to make sure the whole tree is covered
1162  unsigned int i, numChildren = node->getNumChildren();
1163 
1164  for (i = 0; i < numChildren; ++i)
1166 
1167  return ret;
1168 }
1169 
1171 {
1172  //TODO handle case where conflict flag is set before
1173 
1175  {
1176  return *pNodeUnit;
1177  }
1178  else if (pNodeUnit->getInfo() == CSBMLunitInformation::UNKNOWN)
1179  {
1180  *pNodeUnit = ui;
1182  return *pNodeUnit;
1183  }
1184  else //both are known
1185  {
1186  //check for conflict
1187  if (CSBMLunit::isEqual(ui, *pNodeUnit))
1188  return ui;
1189  else //there is a conflict
1190  {
1191  if (ui.getInfo() < pNodeUnit->getInfo()) //the new unit is more reliable
1192  {
1193  *pNodeUnit = ui;
1195  }
1196 
1197  if (!pNodeUnit->isConflict())
1198  {
1199  //only report the conflict for this expression if it was not known before
1200  if (1 > mError) mError = 1;
1201 
1202  mConflictingNodes.insert(node);
1203  }
1204 
1205  pNodeUnit->setConflict(true);
1206  return *pNodeUnit;
1207  }
1208  }
1209 }
1210 
1212  const CSBMLunitInformation & ui,
1213  const CEnvironmentInformation & ei)
1214 {
1215  CSBMLunitInformation ret(this->mSBMLLevel, this->mSBMLVersion);
1216 
1217  if (!node) return ret;
1218 
1219  ret = ui;
1220 
1221  //TODO deal with conflicts
1222 
1223  unsigned int i, numChildren = node->getNumChildren();
1224  std::vector<CSBMLunitInformation> childUnits;
1225  childUnits.resize(numChildren, CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion));
1226 
1227  //first deal with the unit that is passed from above
1229  {
1230  //pass the unit to all child nodes
1231  for (i = 0; i < numChildren; ++i)
1232  childUnits[i] = recursion(node->getChild(i), ui, ei);
1233  }
1234  else
1235  {
1236  //the unit passed from above is unknown. Check if one of the children has a known unit
1237  for (i = 0; i < numChildren; ++i)
1238  {
1239  childUnits[i] = recursion(node->getChild(i), ui, ei);
1240 
1241  if (childUnits[i].getInfo() > CSBMLunitInformation::UNKNOWN)
1242  break;
1243  }
1244 
1245  if (i >= numChildren)
1246  return ret; //do nothing if all units are unknown
1247 
1248  // i is now the index of the first child with known units
1249  // pass this unit information to the return value and all other children
1250  unsigned int tmp = i;
1251  ret = childUnits[tmp];
1252 
1253  for (i = 0; i < numChildren; ++i)
1254  {
1255  if (i == tmp) continue;
1256 
1257  childUnits[i] = recursion(node->getChild(i), ret, ei);
1258  }
1259  }
1260 
1261  //TODO deal with conflicts: if any conflict occured, pass conflict to all children and ret
1262 
1263  return ret;
1264 }
1265 
1267  const CSBMLunitInformation & ui,
1268  const CEnvironmentInformation & ei)
1269 {
1270  CSBMLunitInformation ret(this->mSBMLLevel, this->mSBMLVersion);
1271 
1272  if (!node) return ret;
1273 
1274  ret = ui;
1275 
1276  //TODO deal with conflicts
1277 
1278  unsigned int i, numChildren = node->getNumChildren();
1279  std::vector<CSBMLunitInformation> childUnits;
1280  childUnits.resize(numChildren, CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion));
1281 
1282  // ask all children for their unit
1283  std::vector<int> unknown;
1284  CSBMLunitInformation uu(this->mSBMLLevel, this->mSBMLVersion); //unknown units
1285 
1286  for (i = 0; i < numChildren; ++i) //TODO should stop when we know enough
1287  {
1288  childUnits[i] = recursion(node->getChild(i), uu, ei);
1289 
1290  if (childUnits[i].getInfo() == CSBMLunitInformation::UNKNOWN)
1291  unknown.push_back(i);
1292  }
1293 
1294  //first the case where the parent unit is unknown
1296  {
1297  //if there are children with unknown units we can do nothing
1298  if (unknown.size() > 0)
1299  return ret;
1300 
1301  //determine the units for the parent
1302  ret = childUnits[0];
1303  bool success = true;
1304 
1305  for (i = 1; i < numChildren; ++i)
1306  success &= ret.multiply(childUnits[i]);
1307 
1308  if (success)
1310  else
1312  }
1313  else //parent units are known
1314  {
1315  //if there is more than one child with unknown units we can do nothing
1316  if (unknown.size() > 1)
1317  return ret;
1318 
1319  //if no child has unknown units assume the first (and check for conflicts later)
1320  unsigned int tmp = unknown.size() ? unknown[0] : 0;
1321 
1322  //determine units
1323  CSBMLunitInformation tmpUnit(this->mSBMLLevel, this->mSBMLVersion);
1324  bool success = true;
1325 
1326  for (i = 0; i < numChildren; ++i)
1327  if (i != tmp)
1328  success &= tmpUnit.multiply(childUnits[i]);
1329 
1330  tmpUnit.invert();
1331  success &= tmpUnit.multiply(ui);
1332 
1333  if (success)
1335  else
1337 
1338  //tell child about derived unit
1339  childUnits[tmp] = recursion(node->getChild(tmp), tmpUnit, ei);
1340  //TODO check for conflicts
1341  }
1342 
1343  return ret;
1344 }
1345 
1347  const CSBMLunitInformation & ui,
1348  const CEnvironmentInformation & ei)
1349 {
1350  CSBMLunitInformation ret(this->mSBMLLevel, this->mSBMLVersion);
1351 
1352  if (!node) return ret;
1353 
1354  ret = ui;
1355 
1356  //TODO deal with conflicts
1357 
1358  unsigned int i, numChildren = node->getNumChildren();
1359  assert(numChildren == 2);
1360  std::vector<CSBMLunitInformation> childUnits;
1361  childUnits.resize(numChildren, CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion));
1362 
1363  // ask all children for their unit
1364  std::vector<int> unknown;
1365  CSBMLunitInformation uu(this->mSBMLLevel, this->mSBMLVersion); //unknown units
1366 
1367  for (i = 0; i < numChildren; ++i)
1368  {
1369  childUnits[i] = recursion(node->getChild(i), uu, ei);
1370 
1371  if (childUnits[i].getInfo() == CSBMLunitInformation::UNKNOWN)
1372  unknown.push_back(i);
1373  }
1374 
1375  //first the case where the parent unit is unknown
1377  {
1378  //if there are children with unknown units we can do nothing
1379  if (unknown.size() > 0)
1380  return ret;
1381 
1382  //determine the units for the parent
1383  ret = childUnits[1];
1384  bool success = true;
1385  ret.invert();
1386  success &= ret.multiply(childUnits[0]);
1387 
1388  if (success)
1390  else
1392  }
1393  else
1394  {
1395  //if there is more than one child with unknown units we can do nothing
1396  if (unknown.size() > 1)
1397  return ret;
1398 
1399  //case where nominator is unknown (or assumed unknown)
1400  if (unknown.size() == 0 || unknown[0] == 0)
1401  {
1402  //determine units
1403  CSBMLunitInformation tmpUnit = childUnits[1];
1404  bool success = true;
1405  success &= tmpUnit.multiply(ui);
1406 
1407  if (success)
1409  else
1411 
1412  //tell child about derived unit
1413  childUnits[0] = recursion(node->getChild(0), tmpUnit, ei);
1414  //TODO check for conflicts
1415  }
1416  else //denominator is unknown
1417  {
1418  //determine units
1419  CSBMLunitInformation tmpUnit = ui;
1420  bool success = true;
1421  tmpUnit.invert();
1422  success &= tmpUnit.multiply(childUnits[0]);
1423 
1424  //tmpUnit.multiply(ui);
1425  if (success)
1427  else
1429 
1430  //tell child about derived unit
1431  childUnits[1] = recursion(node->getChild(1), tmpUnit, ei);
1432  }
1433  }
1434 
1435  return ret;
1436 }
1437 
1439  const CSBMLunitInformation & ui,
1440  const CEnvironmentInformation & ei)
1441 {
1442  CSBMLunitInformation ret(this->mSBMLLevel, this->mSBMLVersion);
1443 
1444  if (!node) return ret;
1445 
1446  ret = ui;
1447 
1448  unsigned int numChildren = node->getNumChildren();
1449  assert(numChildren == 2);
1450  std::vector<CSBMLunitInformation> childUnits;
1451  childUnits.resize(numChildren, CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion));
1452 
1453  //the exponent should always be dimensionless
1454  childUnits[1] = recursion(node->getChild(1), CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion, CSBMLunitInformation::DEFAULT), ei);
1455 
1456  //try to evaluate the exponent
1457  EvaluationResult res = evaluate(node->getChild(1));
1458 
1459  if (res.known) //the base of the power could be evaluated
1460  {
1461  //check if the exponent is integer
1462  if (fabs(res.result - floor(res.result + 0.5)) < 1e-100)
1463  {
1464  int intExp = (int) floor(res.result + 0.5);
1465 
1467  {
1468  childUnits[0] = recursion(node->getChild(0), CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion, CSBMLunitInformation::UNKNOWN), ei);
1469  ret = childUnits[0];
1470 
1471  if (ret.getInfo() > CSBMLunitInformation::UNKNOWN)
1472  {//only if the base is known, we can make use of the exponent to calculate the return unit
1473  ret.applyExponent(intExp);
1474  ret.setInfo(CSBMLunitInformation::DERIVED);
1475  }
1476  }
1477  else
1478  {
1479  CSBMLunitInformation tmpUI = ui;
1480  tmpUI.applyExponent(1 / (double)intExp);
1482  childUnits[0] = recursion(node->getChild(0), tmpUI, ei);
1483  }
1484  }
1485  else
1486  {
1487  //TODO perhaps rather conflict???
1488  childUnits[0] = recursion(node->getChild(0), CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion, CSBMLunitInformation::UNKNOWN), ei);
1489  }
1490 
1491  //TODO extend to deal with non integer exponents properly. CSBMLunit needs to support it first.
1492  }
1493  else
1494  {
1495  //the exponent could not be determined as a number
1496 
1497  //special case: exponent is an identifier
1498  if (node->getChild(1)->isName())
1499  {
1501  {
1502  childUnits[0] = recursion(node->getChild(0), CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion, CSBMLunitInformation::UNKNOWN), ei);
1503  ret = childUnits[0];
1504 
1505  if (ret.getInfo() > CSBMLunitInformation::UNKNOWN)
1506  {//only if the base is known, we can make use of the exponent to calculate the return unit
1507  ret.applyExponent(node->getChild(1)->getName(), ei.mFrameStack.size());
1508  ret.setInfo(CSBMLunitInformation::DERIVED);
1509  }
1510  }
1511  else
1512  {
1513  CSBMLunitInformation tmpUI = ui;
1514  tmpUI.applyInverseExponent(node->getChild(1)->getName(), ei.mFrameStack.size());
1516  childUnits[0] = recursion(node->getChild(0), tmpUI, ei);
1517  }
1518  }
1519  else //exponent is neither a number nor a simple identifier. Units of the base are unknown
1520  {
1521  childUnits[0] = recursion(node->getChild(0), CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion, CSBMLunitInformation::UNKNOWN), ei);
1522  }
1523  }
1524 
1525  //heuristics: if the exponent is variable, the base should be dimensionless
1526  //TODO
1527 
1528  //TODO deal with conflicts
1529 
1530  return ret;
1531 }
1532 
1534 {
1535  EvaluationResult ret;
1536  ret.result = 0;
1537  ret.known = false;
1538 
1539  //numbers
1540  if (isNumber(node))
1541  {
1542  ret.result = getValueFromNumberNode(node);
1543  ret.known = true;
1544  }
1545 
1546  //TODO at least plus, minus, uminus
1547 
1548  return ret;
1549 }
1550 
1552  const CSBMLunitInformation & ui,
1553  const CEnvironmentInformation & ei)
1554 {
1555  CSBMLunitInformation ret(this->mSBMLLevel, this->mSBMLVersion);
1556 
1557  if (!node) return ret;
1558 
1559  ret = ui;
1560 
1561  //TODO deal with conflicts
1562 
1563  unsigned int i, numChildren = node->getNumChildren();
1564  std::vector<CSBMLunitInformation> childUnits;
1565  childUnits.resize(numChildren, CSBMLunitInformation(this->mSBMLLevel, this->mSBMLVersion));
1566 
1567  //first do the recursion for the logical parts (the children with uneven index)
1568  for (i = 1; i < numChildren; i += 2)
1570 
1571  //first deal with the unit that is passed from above
1573  {
1574  //pass the unit to all child nodes
1575  for (i = 0; i < numChildren; i += 2)
1576  childUnits[i] = recursion(node->getChild(i), ui, ei);
1577  }
1578  else
1579  {
1580  //the unit passed from above is unknown. Check if one of the children has a known unit
1581  for (i = 0; i < numChildren; i += 2)
1582  {
1583  childUnits[i] = recursion(node->getChild(i), ui, ei);
1584 
1585  if (childUnits[i].getInfo() > CSBMLunitInformation::UNKNOWN)
1586  break;
1587  }
1588 
1589  if (i >= numChildren)
1590  return ret; //do nothing if all units are unknown
1591 
1592  // i is now the index of the first child with known units
1593  // pass this unit information to the return value and all other children
1594  unsigned int tmp = i;
1595  ret = childUnits[tmp];
1596 
1597  for (i = 0; i < numChildren; i += 2)
1598  {
1599  if (i == tmp) continue;
1600 
1601  childUnits[i] = recursion(node->getChild(i), ret, ei);
1602  }
1603  }
1604 
1605  //TODO deal with conflicts: if any conflict occured, pass conflict to all children and ret
1606 
1607  return ret;
1608 }
1609 
1610 //**************************************************
1611 
1613  const CEnvironmentInformation & ei)
1614 {
1615  //try local parameters
1616  if (ei.isReactionScope())
1617  {
1618  std::map<std::string, std::map<std::string, CSBMLunitInformation> >::iterator rit;
1619  rit = mSBMLLocalParametersMap.find(ei.mReactionID);
1620 
1621  if (rit != mSBMLLocalParametersMap.end())
1622  {
1623  std::map<std::string, CSBMLunitInformation>::iterator it;
1624  it = rit->second.find(node);
1625 
1626  if (it != rit->second.end())
1627  return &it->second;
1628  }
1629  else
1630  {
1631  //error
1632  }
1633  }
1634 
1635  //now global objects
1636  std::map<std::string, CSBMLunitInformation>::iterator it;
1637  it = mSBMLObjectsMap.find(node);
1638 
1639  if (it != mSBMLObjectsMap.end())
1640  return &it->second;
1641 
1642  return NULL;
1643 }
1644 
1645 ASTNode* CSBMLunitInterface::resolveVariableName(const std::string & node,
1646  const CEnvironmentInformation & ei)
1647 {
1648  //TODO find ASTNode corresponding to the variable (from top level of frame stack in ei).
1649  if (!ei.mFrameStack.size()) return NULL;
1650 
1651  std::map<std::string, ASTNode*>::const_iterator it;
1652  it = ei.mFrameStack[ei.mFrameStack.size()-1].find(node);
1653 
1654  if (it == ei.mFrameStack[ei.mFrameStack.size()-1].end())
1655  return NULL;
1656  else
1657  return it->second;
1658 }
1659 
1660 FunctionDefinition* CSBMLunitInterface::resolveFunctionName(const std::string & node)
1661 {
1662  if (!mpModel) return NULL;
1663 
1664  return mpModel->getFunctionDefinition(node);
1665 }
1666 
1668 {
1669  //check if the node is in the map
1670  std::map<const ASTNode*, CSBMLunitInformation>::iterator it = mSBMLNumbersMap.find(node);
1671 
1672  if (it != mSBMLNumbersMap.end())
1673  return &it->second;
1674  else
1675  return NULL;
1676 }
1677 
1678 //*************************************************
1679 
1680 bool CSBMLunitInterface::isObject(const ASTNode* node)
1681 {
1682  return (node && node->isName());
1683 }
1684 
1685 bool CSBMLunitInterface::isOperator(const ASTNode* node)
1686 {
1687  //is this a bug in libsbml that power is a function?
1688  return (node &&
1689  (node->isOperator() || node->isRelational() || node->getType() == AST_FUNCTION_POWER));
1690 }
1691 
1692 bool CSBMLunitInterface::isNumber(const ASTNode* node)
1693 {
1694  return (node && (node->isNumber() || node->isRational()));
1695 }
1696 
1697 bool CSBMLunitInterface::isFunctionCall(const ASTNode* node)
1698 {
1699  return (node && (node->getType() == AST_FUNCTION));
1700 }
1701 
1703 {
1704  return (node && (node->getType() > AST_FUNCTION) && (node->getType() <= AST_FUNCTION_TANH));
1705 }
1706 
1708 {
1709  if (!node)
1710  return std::numeric_limits< double >::quiet_NaN();
1711 
1712  switch (node->getType())
1713  {
1714  case AST_INTEGER:
1715  return node->getInteger();
1716  break;
1717 
1718  case AST_REAL:
1719  case AST_REAL_E:
1720  case AST_RATIONAL:
1721  return node->getReal();
1722  break;
1723 
1724  //TODO rational number format
1725  default:
1726  return std::numeric_limits< double >::quiet_NaN();
1727  }
1728 }
1729 
1730 std::string CSBMLunitInterface::getIdentifier(const ASTNode* node)
1731 {
1732  if (node)
1733  return node->getName();
1734  else
1735  return "";
1736 }
static std::string getIdentifier(const ASTNode *node)
EvaluationResult evaluate(const ASTNode *node)
std::vector< unsigned int > local
the units is determined from the sbml defaults
Definition: CSBMLunit.h:108
#define pdelete(p)
Definition: copasi.h:215
std::string mRootObject
if this is not "" then the root of the expression should have the units of the referenced object ...
void handleOneExpression(CExpressionInformation &ei)
try to find out as much as possible about the units from one expression
CSBMLunitInformation * mpSBMLAreaUnit
CSBMLunitInformation * mpSBMLConflictUnit
virtual std::string getDisplayString() const
Definition: CSBMLunit.cpp:210
CSBMLunitInformation recursionPiecewise(const ASTNode *node, const CSBMLunitInformation &ui, const CEnvironmentInformation &ei)
std::vector< CExpressionInformation > mSBMLExpressions
CSBMLunitInformation recursionTimes(const ASTNode *node, const CSBMLunitInformation &ui, const CEnvironmentInformation &ei)
Model * mpModel
the sbml model from which this interface was initialized
void applyExponent(double exp)
apply numeric exponent to the unit.
Definition: CSBMLunit.cpp:104
std::vector< unsigned int > global
std::string mObjectDisplayString
display name of the object that either contains the expression or the the expression applies to...
the units is unknown
Definition: CSBMLunit.h:106
FunctionDefinition * resolveFunctionName(const std::string &node)
static bool isNumber(const ASTNode *node)
determines if the node represents a number (integer, real, or rational)
CSBMLunitInterface(const CSBMLunitInterface &src)
CSBMLunitInformation * mpSBMLLengthUnit
CSBMLunitInformation recursion(const ASTNode *node, const CSBMLunitInformation &ui, const CEnvironmentInformation &ei)
std::vector< std::map< std::string, ASTNode * > > mFrameStack
static bool isEqual(const CSBMLunit &unit1, const CSBMLunit &unit2)
Definition: CSBMLunit.cpp:159
CSBMLunitInformation * getMappedUnitFromIdentifier(const std::string &node, const CEnvironmentInformation &ei)
the units is determined from the model-wide definitions
Definition: CSBMLunit.h:110
void applyInverseExponent(const std::string &id, size_t frame)
apply inverse of symbolic exponent to the unit
Definition: CSBMLunit.cpp:140
bool isReactionScope() const
tells if the current scope involves a reaction (local parameters)
std::map< std::string, CSBMLunitInformation > mSBMLObjectsMap
INFO getInfo() const
get the status information
Definition: CSBMLunit.h:130
std::vector< std::string > getListOfObjectsWithGivenUnitStatus(int status) const
void setConflict(bool c)
set the conflict flag
Definition: CSBMLunit.h:133
bool mPerTime
specifies whether the units of the result of the expression are some unit divided by time ...
CSBMLunitInformation handleTerminalNode(const CSBMLunitInformation &ui, CSBMLunitInformation *pNodeUnit, const ASTNode *node)
int getSymbolicExpExp() const
Definition: CSBMLunit.h:72
std::string mTypeDescription
this indicates the origin of the expression for display purposes, e.g. "kinetic law" ...
the units if provided for a specific object explicitly
Definition: CSBMLunit.h:112
ASTNode * resolveVariableName(const std::string &node, const CEnvironmentInformation &ei)
CSBMLunitInformation mRootUnit
if mRootObject is "" then this member specifies the units of the root of the expression ...
CSBMLunitInformation * getMappedUnitFromNumberNode(const ASTNode *node)
void push(std::map< std::string, ASTNode * > &frame)
push a mapping for function variables onto the frame stack
CSBMLunitInformation * mpSBMLAmountUnit
void initializeDefaultUnits()
initializes the base units with the defaults (moles, seconds, l, m, m^2)
static bool isOperator(const ASTNode *node)
void initializeFromSBMLModel(bool unitsFromModel)
const ASTNode * mpExpression
root node of an expression in the sbml model
long int flag
Definition: f2c.h:52
CSBMLunitInformation * mpSBMLVolumeUnit
UnitDefinition & getSBMLUnitDefinition()
Definition: CSBMLunit.h:69
std::vector< unsigned int > numbers
std::string mReactionId
id of a reaction. If this is set the local parameters of the reaction are in scope for the expression...
void setInfo(INFO info)
set the status information
Definition: CSBMLunit.h:127
std::map< std::string, std::map< std::string, CSBMLunitInformation > > mSBMLLocalParametersMap
static void outputStatistics(const Statistics &stat, bool flag)
static bool isObject(const ASTNode *node)
static bool isFunctionCall(const ASTNode *node)
determines if the node contains a call to a function definition
the units was determined by reasoning
Definition: CSBMLunit.h:114
std::map< const ASTNode *, CSBMLunitInformation > mSBMLNumbersMap
std::vector< unsigned int > all
static double getValueFromNumberNode(const ASTNode *node)
return the value of a number node
static bool isBuiltInFunctionCall(const ASTNode *node)
determines if the node contains a call to a built in function
bool multiply(const CSBMLunit &unit)
Definition: CSBMLunit.cpp:60
if(!yymsg) yymsg
CSBMLunitInformation recursionDivide(const ASTNode *node, const CSBMLunitInformation &ui, const CEnvironmentInformation &ei)
void drop()
drop the top level of the function variables frame stack
std::set< const ASTNode * > mConflictingNodes
contains the (terminal) nodes where a conflict appeared
void invert()
Definition: CSBMLunit.cpp:88
CSBMLunitInformation * mpSBMLTimeUnit
CSBMLunitInformation recursionPower(const ASTNode *node, const CSBMLunitInformation &ui, const CEnvironmentInformation &ei)
bool isConflict() const
retrieve the conflict flag
Definition: CSBMLunit.h:136
const std::set< const ASTNode * > & getListOfConflictingNodes() const
std::string getMessageAboutUnknownUnits() const
std::vector< std::pair< std::string, std::string > > getListOfLocalParametersWithGivenUnitStatus(int status) const
CSBMLunitInformation recursionEqual(const ASTNode *node, const CSBMLunitInformation &ui, const CEnvironmentInformation &ei)