COPASI API  4.16.103
SBMLImporter.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) 2004 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #ifdef WIN32
16 # pragma warning (disable: 4786)
17 # pragma warning (disable: 4243)
18 // warning C4355: 'this' : used in base member initializer list
19 # pragma warning (disable: 4355)
20 #endif // WIN32
21 
22 #define USE_LAYOUT 1
23 
24 #include <iostream>
25 #include <vector>
26 #include <sstream>
27 #include <map>
28 #include <limits>
29 #include <cmath>
30 
31 #include <sbml/SBMLReader.h>
32 #include <sbml/SBMLDocument.h>
33 #include <sbml/Compartment.h>
34 #include <sbml/Species.h>
35 #include <sbml/SpeciesReference.h>
36 #include <sbml/Reaction.h>
37 #include <sbml/LocalParameter.h>
38 
39 #if LIBSBML_VERSION >= 50400
40 #include <sbml/SBMLTransforms.h>
41 #include "IdList.h" // this file is missing from the libSBML distribution!
42 #include <sbml/conversion/ConversionProperties.h>
43 #include <sbml/packages/layout/extension/LayoutModelPlugin.h>
44 #define INIT_DEFAULTS(element) \
45  {\
46  element.initDefaults();\
47  }
48 #else
49 #define INIT_DEFAULTS(element) \
50  {\
51  }
52 #endif // LIBSBML_VERSION
53 
54 #include <sbml/KineticLaw.h>
55 #include <sbml/math/FormulaFormatter.h>
56 #include <sbml/Model.h>
57 #include <sbml/UnitKind.h>
58 #include <sbml/Unit.h>
59 #include <sbml/Parameter.h>
60 #include <sbml/InitialAssignment.h>
61 #include <sbml/Rule.h>
62 #include <sbml/FunctionDefinition.h>
63 #include <sbml/UnitDefinition.h>
64 
65 #include "copasi.h"
66 
67 #include "report/CKeyFactory.h"
68 #include "model/CModel.h"
69 #include "model/CCompartment.h"
70 #include "model/CMetab.h"
71 #include "model/CReaction.h"
72 #include "model/CModelValue.h"
73 #include "model/CEvent.h"
74 #include "function/CNodeK.h"
75 #include "function/CFunctionDB.h"
77 #include "function/CExpression.h"
80 #include "utilities/CCopasiTree.h"
87 #include "commandline/COptions.h"
88 
89 #include "SBMLImporter.h"
90 #include "SBMLUtils.h"
91 #include "ConverterASTNode.h"
94 
96 #include "layout/CListOfLayouts.h"
97 
99 
100 #if LIBSBML_VERSION >= 50903 && LIBSBML_HAS_PACKAGE_COMP
101 
102 #include <sbml/util/PrefixTransformer.h>
103 #include <sbml/packages/comp/extension/CompModelPlugin.h>
104 
105 class CPrefixNameTransformer : public PrefixTransformer
106 {
107 public:
108  CPrefixNameTransformer() {}
109 
110  static void replaceStringInPlace(std::string& subject, const std::string& search,
111  const std::string& replace)
112  {
113  size_t pos = 0;
114 
115  while ((pos = subject.find(search, pos)) != std::string::npos)
116  {
117  subject.replace(pos, search.length(), replace);
118  pos += replace.length();
119  }
120  }
121 
122  static inline std::string &rtrim(std::string &str)
123  {
124  size_t endpos = str.find_last_not_of(" \t");
125 
126  if (std::string::npos != endpos)
127  {
128  str = str.substr(0, endpos + 1);
129  }
130 
131  return str;
132  }
133 
134  const std::string& cleanName(std::string& prefix)
135  {
136  std::replace(prefix.begin(), prefix.end(), '_', ' ');
137  replaceStringInPlace(prefix, " ", " ");
138  rtrim(prefix);
139  return prefix;
140  }
141 
142  virtual int transform(SBase* element)
143  {
144  if (element == NULL || getPrefix().empty())
145  return LIBSBML_OPERATION_SUCCESS;
146 
147  // set up ids
149 
150  // skip local parameters, as they are not renamed
151  if (element->getTypeCode() == SBML_LOCAL_PARAMETER)
152  return LIBSBML_OPERATION_SUCCESS;
153 
154  // setup names
155  if (element->isSetName())
156  {
157  std::stringstream newName;
158  std::string prefix = getPrefix();
159  newName << element->getName() << " (" << cleanName(prefix) << ")";
160  element->setName(newName.str());
161  }
162 
163  return LIBSBML_OPERATION_SUCCESS;
164  }
165 };
166 
167 #endif
168 
169 // static
171 {
172  return
173  x < 0.0 ? -floor(-x + 0.5) : floor(x + 0.5);
174 }
175 
176 // static
177 bool SBMLImporter::areApproximatelyEqual(const double & x, const double & y, const double & t)
178 {
179  double Scale =
180  (fabs(x) + fabs(y)) * t;
181 
182  // Avoid underflow
183  if (Scale < 100.0 * std::numeric_limits< C_FLOAT64 >::min())
184  return true;
185 
186  return 2 * fabs(x - y) < Scale;
187 }
188 
189 std::string getOriginalSBMLId(Parameter* parameter)
190 {
191  if (parameter == NULL) return "";
192 
193  if (!parameter->isSetAnnotation()) return "";
194 
195  XMLNode* node = parameter->getAnnotation();
196 
197  if (node->getNumChildren() < 1) return "";
198 
199  for (unsigned int i = 0; i < node->getNumChildren(); ++i)
200  {
201  const XMLNode& current = node->getChild(i);
202 
203  if (current.getNamespaces().containsUri("http://copasi.org/initialValue"))
204  {
205  return current.getAttrValue("parent");
206  }
207  }
208 
209  return "";
210 }
211 
212 std::string getInitialCNForSBase(SBase* sbase, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
213 {
214  std::map<CCopasiObject*, SBase*>::const_iterator it;
215 
216  for (it = copasi2sbmlmap.begin(); it != copasi2sbmlmap.end(); ++it)
217  {
218  if (it->second != sbase)
219  continue;
220 
221  CMetab *metab = dynamic_cast<CMetab*>(it->first);
222 
223  if (metab != NULL)
224  return metab->getInitialConcentrationReference()->getCN();
225 
226  CCompartment *comp = dynamic_cast<CCompartment*>(it->first);
227 
228  if (comp != NULL)
229  return comp->getInitialValueReference()->getCN();
230 
231  CModelValue *param = dynamic_cast<CModelValue*>(it->first);
232 
233  if (param != NULL)
234  return param->getInitialValueReference()->getCN();
235  }
236 
237  return "";
238 }
239 
240 /**
241  * Creates and returns a Copasi CModel from the SBMLDocument given as argument.
242  */
243 CModel* SBMLImporter::createCModelFromSBMLDocument(SBMLDocument* sbmlDocument, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
244 {
245  Model* sbmlModel = sbmlDocument->getModel();
246 
247  /* Create an empty model and set the title. */
248  this->mpCopasiModel = new CModel(mpDataModel);
249  copasi2sbmlmap[this->mpCopasiModel] = sbmlModel;
255  this->mpCopasiModel->setSBMLId(sbmlModel->getId());
256 
257  unsigned C_INT32 step = 0, totalSteps = 0;
258  size_t hStep = C_INVALID_INDEX;
259 
260  mImportStep = 1;
261 
263 
264  SBMLImporter::importMIRIAM(sbmlModel, this->mpCopasiModel);
265  UnitDefinition *pSubstanceUnits = NULL;
266  UnitDefinition *pTimeUnits = NULL;
267  UnitDefinition *pVolumeUnits = NULL;
268  UnitDefinition *pAreaUnits = NULL;
269  UnitDefinition *pLengthUnits = NULL;
270  /* Set standard units to match the standard units of SBML files. */
271 
272  // for SBML L3 files the default units are defined on the model
273  if (this->mLevel > 2)
274  {
275  this->mpCopasiModel->setAvogadro(6.02214179e23);
276  this->mAvogadroSet = true;
277 
278  // we make copies of the unit definitions so that we do not have to remember
279  // if we created them or not
280  std::string units;
281  UnitDefinition* pUDef = NULL;
282  Unit unit(sbmlModel->getLevel(), sbmlModel->getVersion());
283  INIT_DEFAULTS(unit);
284 
285  if (sbmlModel->isSetSubstanceUnits())
286  {
287  units = sbmlModel->getSubstanceUnits();
288  assert(units != "");
289  pUDef = sbmlModel->getUnitDefinition(units);
290 
291  if (pUDef != NULL)
292  {
293  pSubstanceUnits = new UnitDefinition(*pUDef);
294  }
295  else
296  {
297  if (units == "mole")
298  {
299  unit.setKind(UNIT_KIND_MOLE);
300  unit.setExponent(1);
301  unit.setMultiplier(1.0);
302  unit.setScale(0);
303  pSubstanceUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
304  pSubstanceUnits->addUnit(&unit);
305  }
306  else if (units == "item")
307  {
308  unit.setKind(UNIT_KIND_ITEM);
309  unit.setExponent(1);
310  unit.setMultiplier(1.0);
311  unit.setScale(0);
312  pSubstanceUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
313  pSubstanceUnits->addUnit(&unit);
314  }
315  else if (units == "dimensionless")
316  {
317  unit.setKind(UNIT_KIND_DIMENSIONLESS);
318  unit.setExponent(1);
319  unit.setMultiplier(1.0);
320  unit.setScale(0);
321  pSubstanceUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
322  pSubstanceUnits->addUnit(&unit);
323  }
324  else
325  {
326  std::string message = "COPASI can't handle substance unit \"" + units + "\". Setting unit for substances to dimensionless.";
327  CCopasiMessage(CCopasiMessage::WARNING, message.c_str());
328  unit.setKind(UNIT_KIND_DIMENSIONLESS);
329  unit.setExponent(1);
330  unit.setMultiplier(1.0);
331  unit.setScale(0);
332  pSubstanceUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
333  pSubstanceUnits->addUnit(&unit);
334  }
335  }
336  }
337 
338  if (sbmlModel->isSetTimeUnits())
339  {
340  units = sbmlModel->getTimeUnits();
341  assert(units != "");
342  pUDef = sbmlModel->getUnitDefinition(units);
343 
344  if (pUDef != NULL)
345  {
346  pTimeUnits = new UnitDefinition(*pUDef);
347  }
348  else
349  {
350  if (units == "second")
351  {
352  unit.setKind(UNIT_KIND_SECOND);
353  unit.setExponent(1);
354  unit.setMultiplier(1.0);
355  unit.setScale(0);
356  pTimeUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
357  pTimeUnits->addUnit(&unit);
358  }
359  else if (units == "dimensionless")
360  {
361  unit.setKind(UNIT_KIND_DIMENSIONLESS);
362  unit.setExponent(1);
363  unit.setMultiplier(1.0);
364  unit.setScale(0);
365  pTimeUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
366  pTimeUnits->addUnit(&unit);
367  }
368  else
369  {
370  std::string message = "COPASI can't handle time unit \"" + units + "\". Setting unit for time to dimensionless.";
371  CCopasiMessage(CCopasiMessage::WARNING, message.c_str());
372  unit.setKind(UNIT_KIND_DIMENSIONLESS);
373  unit.setExponent(1);
374  unit.setMultiplier(1.0);
375  unit.setScale(0);
376  pTimeUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
377  pTimeUnits->addUnit(&unit);
378  }
379  }
380  }
381 
382  if (sbmlModel->isSetVolumeUnits())
383  {
384  units = sbmlModel->getVolumeUnits();
385  assert(units != "");
386  pUDef = sbmlModel->getUnitDefinition(units);
387 
388  if (pUDef != NULL)
389  {
390  pVolumeUnits = new UnitDefinition(*pUDef);
391  }
392  else
393  {
394  if (units == "litre")
395  {
396  unit.setKind(UNIT_KIND_LITRE);
397  unit.setExponent(1);
398  unit.setMultiplier(1.0);
399  unit.setScale(0);
400  pVolumeUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
401  pVolumeUnits->addUnit(&unit);
402  }
403  else if (units == "dimensionless")
404  {
405  unit.setKind(UNIT_KIND_DIMENSIONLESS);
406  unit.setExponent(1);
407  unit.setMultiplier(1.0);
408  unit.setScale(0);
409  pVolumeUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
410  pVolumeUnits->addUnit(&unit);
411  }
412  else
413  {
414  std::string message = "COPASI can't handle volume unit \"" + units + "\". Setting unit for volume to dimensionless.";
415  CCopasiMessage(CCopasiMessage::WARNING, message.c_str());
416  unit.setKind(UNIT_KIND_DIMENSIONLESS);
417  unit.setExponent(1);
418  unit.setMultiplier(1.0);
419  unit.setScale(0);
420  pVolumeUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
421  pVolumeUnits->addUnit(&unit);
422  }
423  }
424  }
425 
426  if (sbmlModel->isSetAreaUnits())
427  {
428  units = sbmlModel->getAreaUnits();
429  assert(units != "");
430  pUDef = sbmlModel->getUnitDefinition(units);
431 
432  if (pUDef != NULL)
433  {
434  pAreaUnits = new UnitDefinition(*pUDef);
435  }
436  else
437  {
438  if (units == "dimensionless")
439  {
440  unit.setKind(UNIT_KIND_DIMENSIONLESS);
441  unit.setExponent(1);
442  unit.setMultiplier(1.0);
443  unit.setScale(0);
444  pAreaUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
445  pAreaUnits->addUnit(&unit);
446  }
447  else
448  {
449  std::string message = "COPASI can't handle area unit \"" + units + "\". Setting unit for area to dimensionless.";
450  CCopasiMessage(CCopasiMessage::WARNING, message.c_str());
451  unit.setKind(UNIT_KIND_DIMENSIONLESS);
452  unit.setExponent(1);
453  unit.setMultiplier(1.0);
454  unit.setScale(0);
455  pAreaUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
456  pAreaUnits->addUnit(&unit);
457  }
458  }
459  }
460 
461  if (sbmlModel->isSetLengthUnits())
462  {
463  units = sbmlModel->getLengthUnits();
464  assert(units != "");
465  pUDef = sbmlModel->getUnitDefinition(units);
466 
467  if (pUDef != NULL)
468  {
469  pLengthUnits = new UnitDefinition(*pUDef);
470  }
471  else
472  {
473  if (units == "litre")
474  {
475  unit.setKind(UNIT_KIND_LITRE);
476  unit.setExponent(1);
477  unit.setMultiplier(1.0);
478  unit.setScale(0);
479  pLengthUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
480  pLengthUnits->addUnit(&unit);
481  }
482  else if (units == "metre")
483  {
484  unit.setKind(UNIT_KIND_METRE);
485  unit.setExponent(1);
486  unit.setMultiplier(1.0);
487  unit.setScale(0);
488  pLengthUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
489  pLengthUnits->addUnit(&unit);
490  }
491  else if (units == "dimensionless")
492  {
493  unit.setKind(UNIT_KIND_DIMENSIONLESS);
494  unit.setExponent(1);
495  unit.setMultiplier(1.0);
496  unit.setScale(0);
497  pLengthUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
498  pLengthUnits->addUnit(&unit);
499  }
500  else
501  {
502  std::string message = "COPASI can't handle length unit \"" + units + "\". Setting unit for length to dimensionless.";
503  CCopasiMessage(CCopasiMessage::WARNING, message.c_str());
504  unit.setKind(UNIT_KIND_DIMENSIONLESS);
505  unit.setExponent(1);
506  unit.setMultiplier(1.0);
507  unit.setScale(0);
508  pLengthUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
509  pLengthUnits->addUnit(&unit);
510  }
511  }
512  }
513  }
514  else
515  {
516  if (sbmlModel->getNumUnitDefinitions() != 0)
517  {
518  unsigned int counter;
519 
520  // we make copies o the unit definitions so that we do not have to remember
521  // if we created them or not
522  for (counter = 0; counter < sbmlModel->getNumUnitDefinitions(); counter++)
523  {
524  UnitDefinition* uDef = sbmlModel->getUnitDefinition(counter);
525  std::string unitId = uDef->getId();
526 
527  if (unitId == "substance")
528  {
529  pSubstanceUnits = new UnitDefinition(*uDef);
530  }
531  else if (unitId == "time")
532  {
533  pTimeUnits = new UnitDefinition(*uDef);
534  }
535  else if (unitId == "volume")
536  {
537  pVolumeUnits = new UnitDefinition(*uDef);
538  }
539  else if ((unitId == "area"))
540  {
541  pAreaUnits = new UnitDefinition(*uDef);
542  }
543  else if ((unitId == "length"))
544  {
545  pLengthUnits = new UnitDefinition(*uDef);
546  }
547  }
548  }
549  }
550 
551  // create the default units if some unit has not been specified
552  if (pSubstanceUnits == NULL)
553  {
554  if (this->mLevel > 2)
555  {
556  // issue a warning
557  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 91, "substance", "mole" , "substance");
558  }
559 
560  // create the default units
561  pSubstanceUnits = new UnitDefinition(this->mLevel, this->mVersion);
562  pSubstanceUnits->setId("dummy_substance");
563  pSubstanceUnits->setName("dummy_substance");
564  Unit* pUnit = pSubstanceUnits->createUnit();
565  pUnit->setKind(UNIT_KIND_MOLE);
566  pUnit->setExponent(1);
567  pUnit->setMultiplier(1.0);
568  pUnit->setScale(0);
569  }
570 
571  if (pTimeUnits == NULL)
572  {
573  if (this->mLevel > 2)
574  {
575  // issue a warning
576  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 91, "time", "second" , "time");
577  }
578 
579  // create the default units
580  pTimeUnits = new UnitDefinition(this->mLevel, this->mVersion);
581  pTimeUnits->setId("dummy_time");
582  pTimeUnits->setName("dummy_time");
583  Unit* pUnit = pTimeUnits->createUnit();
584  pUnit->setKind(UNIT_KIND_SECOND);
585  pUnit->setExponent(1);
586  pUnit->setMultiplier(1.0);
587  pUnit->setScale(0);
588  }
589 
590  if (pVolumeUnits == NULL)
591  {
592  if (this->mLevel > 2)
593  {
594  // issue a warning
595  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 91, "volume", "litre" , "volume");
596  }
597 
598  // create the default units
599  pVolumeUnits = new UnitDefinition(this->mLevel, this->mVersion);
600  pVolumeUnits->setId("dummy_volume");
601  pVolumeUnits->setName("dummy_volume");
602  Unit* pUnit = pVolumeUnits->createUnit();
603  pUnit->setKind(UNIT_KIND_LITRE);
604  pUnit->setExponent(1);
605  pUnit->setMultiplier(1.0);
606  pUnit->setScale(0);
607  }
608 
609  if (pAreaUnits == NULL)
610  {
611  if (this->mLevel > 2)
612  {
613  // issue a warning
614  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 91, "area", "m^2" , "area");
615  }
616 
617  // create the default units
618  pAreaUnits = new UnitDefinition(this->mLevel, this->mVersion);
619  pAreaUnits->setId("dummy_area");
620  pAreaUnits->setName("dummy_area");
621  Unit* pUnit = pAreaUnits->createUnit();
622  pUnit->setKind(UNIT_KIND_METRE);
623  pUnit->setExponent(2);
624  pUnit->setMultiplier(1.0);
625  pUnit->setScale(0);
626  }
627 
628  if (pLengthUnits == NULL)
629  {
630  if (this->mLevel > 2)
631  {
632  // issue a warning
633  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 91, "length", "m" , "length");
634  }
635 
636  // create the default units
637  pLengthUnits = new UnitDefinition(this->mLevel, this->mVersion);
638  pLengthUnits->setId("dummy_length");
639  pLengthUnits->setName("dummy_length");
640  Unit* pUnit = pLengthUnits->createUnit();
641  pUnit->setKind(UNIT_KIND_METRE);
642  pUnit->setExponent(1);
643  pUnit->setMultiplier(1.0);
644  pUnit->setScale(0);
645  }
646 
647  // now we have some common code to actually import the units
648 
649  // handle the substance units
650  assert(pSubstanceUnits != NULL);
651 
652  if (pSubstanceUnits != NULL)
653  {
654  std::pair<CModel::QuantityUnit, bool> qUnit;
655 
656  try
657  {
658  qUnit = this->handleSubstanceUnit(pSubstanceUnits);
659  }
660  catch (...)
661  {
662  std::ostringstream os;
663  os << "Error while importing substance units.";
664 
665  // check if the last message on the stack is an exception
666  // and if so, add the message text to the current exception
668  {
669  // we only want the message, not the timestamp line
670  std::string text = CCopasiMessage::peekLastMessage().getText();
671  os << text.substr(text.find("\n"));
672  }
673 
674  CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
675  }
676 
677  if (qUnit.second == false)
678  {
679  // the unit could not be handled, give an error message and
680  // set the units to mole
681  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 66, "substance", "Mole");
683  }
684  else
685  {
686  this->mpCopasiModel->setQuantityUnit(qUnit.first);
687  }
688 
689  // check if the extends units are set and if they are equal to the substance units
690  // otherwise issue a warning
691  if (this->mLevel > 2)
692  {
693  if (!sbmlModel->isSetExtentUnits())
694  {
696  }
697  else
698  {
699  const UnitDefinition* pExtendsUnits = sbmlModel->getUnitDefinition(sbmlModel->getExtentUnits());
700 
701  if (pExtendsUnits != NULL)
702  {
703  if (!areSBMLUnitDefinitionsIdentical(pSubstanceUnits, pExtendsUnits))
704  {
706  }
707  }
708  else
709  {
710  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 66, "extends", "the same units as the substances");
711  }
712  }
713  }
714 
715  delete pSubstanceUnits;
716  pSubstanceUnits = NULL;
717  }
718 
719  // handle the time units
720  assert(pTimeUnits != NULL);
721 
722  if (pTimeUnits != NULL)
723  {
724  std::pair<CModel::TimeUnit, bool> tUnit;
725 
726  try
727  {
728  tUnit = this->handleTimeUnit(pTimeUnits);
729  }
730  catch (...)
731  {
732  std::ostringstream os;
733  os << "Error while importing time units.";
734 
735  // check if the last message on the stack is an exception
736  // and if so, add the message text to the current exception
738  {
739  // we only want the message, not the timestamp line
740  std::string text = CCopasiMessage::peekLastMessage().getText();
741  os << text.substr(text.find("\n"));
742  }
743 
744  CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
745  }
746 
747  if (tUnit.second == false)
748  {
749  // the unit could not be handled, give an error message and
750  // set the units to second
751  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 66, "time", "second");
753  }
754  else
755  {
756  this->mpCopasiModel->setTimeUnit(tUnit.first);
757  }
758 
759  delete pTimeUnits;
760  pTimeUnits = NULL;
761  }
762 
763  // handle the volume units
764  assert(pVolumeUnits != NULL);
765 
766  if (pVolumeUnits != NULL)
767  {
768  std::pair<CModel::VolumeUnit, bool> vUnit;
769 
770  try
771  {
772  vUnit = this->handleVolumeUnit(pVolumeUnits);
773  }
774  catch (...)
775  {
776  std::ostringstream os;
777  os << "Error while importing volume units.";
778 
779  // check if the last message on the stack is an exception
780  // and if so, add the message text to the current exception
782  {
783  // we only want the message, not the timestamp line
784  std::string text = CCopasiMessage::peekLastMessage().getText();
785  os << text.substr(text.find("\n"));
786  }
787 
788  CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
789  }
790 
791  if (vUnit.second == false)
792  {
793  // the unit could not be handled, give an error message and
794  // set the units to litre
795  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 66, "volume", "litre");
797  }
798  else
799  {
800  this->mpCopasiModel->setVolumeUnit(vUnit.first);
801  }
802 
803  delete pVolumeUnits;
804  pVolumeUnits = NULL;
805  }
806 
807  // handle the area units
808  assert(pAreaUnits != NULL);
809 
810  if (pAreaUnits != NULL)
811  {
812  std::pair<CModel::AreaUnit, bool> vUnit;
813 
814  try
815  {
816  vUnit = this->handleAreaUnit(pAreaUnits);
817  }
818  catch (...)
819  {
820  std::ostringstream os;
821  os << "Error while importing area units.";
822 
823  // check if the last message on the stack is an exception
824  // and if so, add the message text to the current exception
826  {
827  // we only want the message, not the timestamp line
828  std::string text = CCopasiMessage::peekLastMessage().getText();
829  os << text.substr(text.find("\n"));
830  }
831 
832  CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
833  }
834 
835  if (vUnit.second == false)
836  {
837  // the unit could not be handled, give an error message and
838  // set the units to litre
839  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 66, "area", "square meter");
841  }
842  else
843  {
844  this->mpCopasiModel->setAreaUnit(vUnit.first);
845  }
846 
847  delete pAreaUnits;
848  pAreaUnits = NULL;
849  }
850 
851  // handle the length units
852  assert(pLengthUnits != NULL);
853 
854  if (pLengthUnits != NULL)
855  {
856  std::pair<CModel::LengthUnit, bool> vUnit;
857 
858  try
859  {
860  vUnit = this->handleLengthUnit(pLengthUnits);
861  }
862  catch (...)
863  {
864  std::ostringstream os;
865  os << "Error while importing length units.";
866 
867  // check if the last message on the stack is an exception
868  // and if so, add the message text to the current exception
870  {
871  // we only want the message, not the timestamp line
872  std::string text = CCopasiMessage::peekLastMessage().getText();
873  os << text.substr(text.find("\n"));
874  }
875 
876  CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
877  }
878 
879  if (vUnit.second == false)
880  {
881  // the unit could not be handled, give an error message and
882  // set the units to litre
883  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 66, "length", "meter");
885  }
886  else
887  {
888  this->mpCopasiModel->setLengthUnit(vUnit.first);
889  }
890 
891  delete pLengthUnits;
892  pLengthUnits = NULL;
893  }
894 
895  // go through all compartments and species and check if the units are
896  // consistent
897  checkElementUnits(sbmlModel, this->mpCopasiModel, this->mLevel, this->mVersion);
898  std::string title;
899 
900  if (this->isStochasticModel(sbmlModel))
901  {
903  }
904  else
905  {
907  }
908 
909  SBMLImporter::importNotes(this->mpCopasiModel, sbmlModel);
910 
911  title = sbmlModel->getName();
912 
913  if (title == "")
914  {
915  title = "NoName";
916  }
917 
918  size_t idCount = 0;
919 
920  while (mpDataModel->getObject("Model=" + title) != NULL)
921  {
922  std::stringstream str; str << sbmlModel->getName() << "_" << ++idCount;
923  title = str.str();
924  }
925 
926  this->mpCopasiModel->setObjectName(title.c_str());
927  // fill the set of SBML species reference ids because
928  // we need this to check for references to species references in all expressions
929  // as long as we do not support these references
930  SBMLImporter::updateSBMLSpeciesReferenceIds(sbmlModel, this->mSBMLSpeciesReferenceIds);
931 
932  /* import the functions */
933  unsigned int counter;
934  CCopasiVectorN< CFunction >* functions = &(this->functionDB->loadedFunctions());
935 
936  size_t num = (*functions).size();
937 
938  if (mpImportHandler)
939  {
940  mImportStep = 2;
941 
942  if (!mpImportHandler->progressItem(mhImportStep)) return NULL;
943 
944  totalSteps = (unsigned C_INT32) num;
945  hStep = mpImportHandler->addItem("Importing function definitions",
946  step,
947  &totalSteps);
948  }
949 
950  this->sbmlIdMap.clear();
951 
952  for (counter = 0; counter < num; ++counter)
953  {
954  CFunction* tree = (*functions)[counter];
955 
956  if (!tree->getSBMLId().empty())
957  {
958  this->sbmlIdMap[tree] = tree->getSBMLId();
959  tree->setSBMLId("");
960  }
961  }
962 
963  CFunctionDB* pTmpFunctionDB = this->importFunctionDefinitions(sbmlModel, copasi2sbmlmap);
964  // try to find global parameters that represent avogadros number
966 
967  std::map<std::string, CCompartment*> compartmentMap;
968 
969  /* Create the compartments */
970  num = sbmlModel->getNumCompartments();
971 
972  if (mpImportHandler)
973  {
974  mpImportHandler->finishItem(hStep);
975  mImportStep = 3;
976 
977  if (!mpImportHandler->progressItem(mhImportStep)) return NULL;
978 
979  step = 0;
980  totalSteps = (unsigned C_INT32) num;
981  hStep = mpImportHandler->addItem("Importing compartments...",
982  step,
983  &totalSteps);
984  }
985 
986  for (counter = 0; counter < num; counter++)
987  {
988  Compartment* sbmlCompartment = sbmlModel->getCompartment(counter);
989 
990  if (sbmlCompartment == NULL)
991  {
992  fatalError();
993  }
994 
995  CCompartment* pCopasiCompartment = NULL;
996 
997  try
998  {
999  pCopasiCompartment = this->createCCompartmentFromCompartment(sbmlCompartment, this->mpCopasiModel, copasi2sbmlmap, sbmlModel);
1000  }
1001  catch (...)
1002  {
1003  std::ostringstream os;
1004  os << "Error while importing compartment \"";
1005  os << sbmlCompartment->getId() << "\".";
1006 
1007  // check if the last message on the stack is an exception
1008  // and if so, add the message text to the current exception
1010  {
1011  // we only want the message, not the timestamp line
1012  std::string text = CCopasiMessage::peekLastMessage().getText();
1013  os << text.substr(text.find("\n"));
1014  }
1015 
1016  CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
1017  }
1018 
1019  std::string key = sbmlCompartment->getId();
1020 
1021  if (mLevel == 1)
1022  {
1023  key = sbmlCompartment->getName();
1024  }
1025 
1026  compartmentMap[key] = pCopasiCompartment;
1027  ++step;
1028 
1029  if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return NULL;
1030  }
1031 
1032  /* Create all species */
1033  num = sbmlModel->getNumSpecies();
1034 
1035  if (mpImportHandler)
1036  {
1037  mpImportHandler->finishItem(hStep);
1038  mImportStep = 4;
1039 
1040  if (!mpImportHandler->progressItem(mhImportStep)) return NULL;
1041 
1042  step = 0;
1043  totalSteps = (unsigned C_INT32) num;
1044  hStep = mpImportHandler->addItem("Importing species...",
1045  step,
1046  &totalSteps);
1047  }
1048 
1049  for (counter = 0; counter < num; ++counter)
1050  {
1051  Species* sbmlSpecies = sbmlModel->getSpecies(counter);
1052 
1053  if (sbmlSpecies == NULL)
1054  {
1055  fatalError();
1056  }
1057 
1058  CCompartment* pCopasiCompartment = compartmentMap[sbmlSpecies->getCompartment()];
1059 
1060  if (pCopasiCompartment != NULL)
1061  {
1062  CMetab* pCopasiMetabolite = NULL;
1063 
1064  try
1065  {
1066  pCopasiMetabolite = this->createCMetabFromSpecies(sbmlSpecies, this->mpCopasiModel, pCopasiCompartment, copasi2sbmlmap, sbmlModel);
1067  }
1068  catch (...)
1069  {
1070  std::ostringstream os;
1071  os << "Error while importing species \"";
1072  os << sbmlSpecies->getId() << "\".";
1073 
1074  // check if the last message on the stack is an exception
1075  // and if so, add the message text to the current exception
1077  {
1078  // we only want the message, not the timestamp line
1079  std::string text = CCopasiMessage::peekLastMessage().getText();
1080  os << text.substr(text.find("\n"));
1081  }
1082 
1083  CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
1084  }
1085 
1086  std::string key;
1087 
1088  if (this->mLevel == 1)
1089  {
1090  key = sbmlSpecies->getName();
1091  }
1092  else
1093  {
1094  key = sbmlSpecies->getId();
1095  }
1096 
1097  this->speciesMap[key] = pCopasiMetabolite;
1098  }
1099  else
1100  {
1101  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 5 , sbmlSpecies->getCompartment().c_str(), sbmlSpecies->getId().c_str());
1102  }
1103 
1104  ++step;
1105 
1106  if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return NULL;
1107  }
1108 
1109  /* Create the global Parameters */
1110  num = sbmlModel->getNumParameters();
1111 
1112  if (mpImportHandler)
1113  {
1114  mpImportHandler->finishItem(hStep);
1115  mImportStep = 5;
1116 
1117  if (!mpImportHandler->progressItem(mhImportStep)) return NULL;
1118 
1119  step = 0;
1120  totalSteps = (unsigned C_INT32) num;
1121  hStep = mpImportHandler->addItem("Importing global parameters...",
1122  step,
1123  &totalSteps);
1124  }
1125 
1126  for (counter = 0; counter < num; ++counter)
1127  {
1128  Parameter* sbmlParameter = sbmlModel->getParameter(counter);
1129 
1130  if (sbmlParameter == NULL)
1131  {
1132  fatalError();
1133  }
1134 
1135  try
1136  {
1137  std::string sbmlId = getOriginalSBMLId(sbmlParameter);
1138 
1139  // don't import parameter if it is one of our special ones instead get the cn to the target
1140  if (!sbmlId.empty())
1141  {
1142  Species* species = sbmlModel->getSpecies(sbmlId);
1143 
1144  if (species != NULL)
1145  {
1146  mKnownInitalValues[sbmlParameter->getId()] = getInitialCNForSBase(species, copasi2sbmlmap);
1147  continue;
1148  }
1149 
1150  Compartment *comp = sbmlModel->getCompartment(sbmlId);
1151 
1152  if (comp != NULL)
1153  {
1154  mKnownInitalValues[sbmlParameter->getId()] = getInitialCNForSBase(comp, copasi2sbmlmap);
1155  continue;
1156  }
1157 
1158  Parameter* param = sbmlModel->getParameter(sbmlId);
1159 
1160  if (param != NULL)
1161  {
1162  mKnownInitalValues[sbmlParameter->getId()] = getInitialCNForSBase(param, copasi2sbmlmap);
1163  continue;
1164  }
1165  }
1166 
1167  this->createCModelValueFromParameter(sbmlParameter, this->mpCopasiModel, copasi2sbmlmap);
1168  }
1169  catch (...)
1170  {
1171  std::ostringstream os;
1172  os << "Error while importing parameter \"";
1173  os << sbmlParameter->getId() << "\".";
1174 
1175  // check if the last message on the stack is an exception
1176  // and if so, add the message text to the current exception
1178  {
1179  // we only want the message, not the timestamp line
1180  std::string text = CCopasiMessage::peekLastMessage().getText();
1181  os << text.substr(text.find("\n"));
1182  }
1183 
1184  CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
1185  }
1186 
1187  ++step;
1188 
1189  if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return NULL;
1190  }
1191 
1192  // the information that is collected here is used in the creation of the reactions to
1193  // adjust the stoichiometry of substrates and products with the conversion factos from
1194  // the SBML model
1195 
1196  // now that all parameters have been imported, we can check if the model
1197  // defines a global conversion factor
1198  if (this->mLevel > 2 && sbmlModel->isSetConversionFactor())
1199  {
1200  this->mConversionFactorFound = true;
1201  std::string id = sbmlModel->getConversionFactor();
1202  assert(id != "");
1203  std::map<std::string, const CModelValue*>::const_iterator pos = this->mSBMLIdModelValueMap.find(id);
1204  assert(pos != this->mSBMLIdModelValueMap.end());
1205 
1206  if (pos != this->mSBMLIdModelValueMap.end())
1207  {
1208  this->mpModelConversionFactor = pos->second;
1209  }
1210  else
1211  {
1212  fatalError();
1213  }
1214  }
1215 
1216  // we need to go through all species again and find conversion factors
1217  // right now I don't want to change the order of importing model values
1218  // and species since I am not sure what implications that would have
1219  if (this->mLevel > 2)
1220  {
1221  num = sbmlModel->getNumSpecies();
1222 
1223  for (counter = 0; counter < num; ++counter)
1224  {
1225  Species* sbmlSpecies = sbmlModel->getSpecies(counter);
1226 
1227  if (sbmlSpecies->isSetConversionFactor())
1228  {
1229  this->mConversionFactorFound = true;
1230  std::map<std::string, const CModelValue*>::const_iterator pos = this->mSBMLIdModelValueMap.find(sbmlSpecies->getConversionFactor());
1231 
1232  if (pos != this->mSBMLIdModelValueMap.end())
1233  {
1234  this->mSpeciesConversionParameterMap[sbmlSpecies->getId()] = pos->second;
1235  }
1236  }
1237  }
1238  }
1239 
1240  if (!this->mIgnoredParameterUnits.empty())
1241  {
1242  std::ostringstream os;
1243  std::vector<std::string>::iterator errorIt = this->mIgnoredParameterUnits.begin();
1244 
1245  while (errorIt != this->mIgnoredParameterUnits.end())
1246  {
1247  os << *errorIt << ", ";
1248  ++errorIt;
1249  }
1250 
1251  std::string s = os.str();
1252  CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 26, s.substr(0, s.size() - 2).c_str());
1253  }
1254 
1255  /* Create all reactions */
1256  num = sbmlModel->getNumReactions();
1257 
1258  if (mpImportHandler)
1259  {
1260  mpImportHandler->finishItem(hStep);
1261  mImportStep = 6;
1262 
1263  if (!mpImportHandler->progressItem(mhImportStep)) return NULL;
1264 
1265  step = 0;
1266  totalSteps = (unsigned C_INT32) num;
1267  hStep = mpImportHandler->addItem("Importing reactions...",
1268  step,
1269  &totalSteps);
1270  }
1271 
1272  this->mDivisionByCompartmentReactions.clear();
1273  this->mFastReactions.clear();
1275 
1276  for (counter = 0; counter < num; counter++)
1277  {
1278  try
1279  {
1280  this->createCReactionFromReaction(sbmlModel->getReaction(counter), sbmlModel, this->mpCopasiModel, copasi2sbmlmap, pTmpFunctionDB);
1281  }
1282  catch (...)
1283  {
1284  std::ostringstream os;
1285  os << "Error while importing reaction \"";
1286  os << sbmlModel->getReaction(counter)->getId() << "\".";
1287 
1288  // check if the last message on the stack is an exception
1289  // and if so, add the message text to the current exception
1291  {
1292  // we only want the message, not the timestamp line
1293  std::string text = CCopasiMessage::peekLastMessage().getText();
1294  os << text.substr(text.find("\n"));
1295  }
1296 
1297  CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
1298  }
1299 
1300  ++step;
1301 
1302  if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return NULL;
1303  }
1304 
1305  if (!this->mDivisionByCompartmentReactions.empty())
1306  {
1307  // create the error message
1308  std::string idList;
1309  std::set<std::string>::const_iterator it = this->mDivisionByCompartmentReactions.begin();
1310  std::set<std::string>::const_iterator endit = this->mDivisionByCompartmentReactions.end();
1311 
1312  while (it != endit)
1313  {
1314  idList += (*it);
1315  idList += ", ";
1316  ++it;
1317  }
1318 
1319  idList = idList.substr(0, idList.length() - 2);
1320  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 17, idList.c_str());
1321  }
1322 
1323  if (!this->mFastReactions.empty())
1324  {
1325  // create the error message
1326  std::string idList;
1327  std::set<std::string>::const_iterator it = this->mFastReactions.begin();
1328  std::set<std::string>::const_iterator endit = this->mFastReactions.end();
1329 
1330  while (it != endit)
1331  {
1332  idList += (*it);
1333  idList += ", ";
1334  ++it;
1335  }
1336 
1337  idList = idList.substr(0, idList.length() - 2);
1338  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 29, idList.c_str());
1339  }
1340 
1341  if (!this->mReactionsWithReplacedLocalParameters.empty())
1342  {
1343  // create the error message
1344  std::string idList;
1345  std::set<std::string>::const_iterator it = this->mReactionsWithReplacedLocalParameters.begin();
1346  std::set<std::string>::const_iterator endit = this->mReactionsWithReplacedLocalParameters.end();
1347 
1348  while (it != endit)
1349  {
1350  idList += (*it);
1351  idList += ", ";
1352  ++it;
1353  }
1354 
1355  idList = idList.substr(0, idList.length() - 2);
1356  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 87, idList.c_str());
1357  }
1358 
1359  // import the initial assignments
1360  // we do this after the reactions since intial assignments can reference reaction ids.
1361  if (mpImportHandler)
1362  {
1363  mpImportHandler->finishItem(hStep);
1364  mImportStep = 7;
1365 
1366  if (!mpImportHandler->progressItem(mhImportStep)) return NULL;
1367  }
1368 
1369  importInitialAssignments(sbmlModel, copasi2sbmlmap, this->mpCopasiModel);
1370 
1371  // import all rules
1372  // we have to import them after the reactions since they can reference a reaction
1373  // id which has to be known when importing the rule
1374 
1375  /* Create the rules */
1376  // rules should be imported after reactions because we use the mSBMLSpeciesReferenceIds to determine if a
1377  // rule changes a species reference (stoichiometry
1378  // Since COPASI does not support this, we need to ignore the rule
1379  this->areRulesUnique(sbmlModel);
1380  num = sbmlModel->getNumRules();
1381 
1382  if (mpImportHandler)
1383  {
1384  mpImportHandler->finishItem(hStep);
1385  mImportStep = 8;
1386 
1387  if (!mpImportHandler->progressItem(mhImportStep)) return NULL;
1388 
1389  step = 0;
1390  totalSteps = (unsigned C_INT32) num;
1391  hStep = mpImportHandler->addItem("Importing global parameters...",
1392  step,
1393  &totalSteps);
1394  }
1395 
1396  for (counter = 0; counter < num; ++counter)
1397  {
1398  Rule* sbmlRule = sbmlModel->getRule(counter);
1399 
1400  if (sbmlRule == NULL)
1401  {
1402  fatalError();
1403  }
1404 
1405  // improve the error message a bit.
1406  try
1407  {
1408  this->importSBMLRule(sbmlRule, copasi2sbmlmap, sbmlModel);
1409  }
1410  catch (...)
1411  {
1412  std::ostringstream os;
1413  os << "Error while importing rule for variable \"";
1414  os << sbmlRule->getVariable() << "\".";
1415 
1416  // check if the last message on the stack is an exception
1417  // and if so, add the message text to the current exception
1419  {
1420  // we only want the message, not the time stamp line
1421  std::string text = CCopasiMessage::peekLastMessage().getText();
1422  os << text.substr(text.find("\n"));
1423  }
1424 
1425  CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
1426  }
1427 
1428  ++step;
1429 
1430  if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return NULL;
1431  }
1432 
1433  if (sbmlModel->getNumConstraints() > 0)
1434  {
1436  }
1437 
1438  // TODO Create all constraints
1439  // TODO Since we don't have constraints yet, there is no code here.
1440  // TODO When implementing import of constraints, don't forget to replace calls to
1441  // TODO explicitly time dependent functions in the constraints math expression.
1442 
1443  // import all events
1444  // events should be imported after reactions because we use the mSBMLSpeciesReferenceIds to determine if an
1445  // event assignment changes a species reference (stoichiometry
1446  // Since COPASI does not support this, we need to ignore the event assignment
1447  if (mpImportHandler)
1448  {
1449  mpImportHandler->finishItem(hStep);
1450  mImportStep = 9;
1451 
1452  if (!mpImportHandler->progressItem(mhImportStep)) return NULL;
1453  }
1454 
1455  this->importEvents(sbmlModel, this->mpCopasiModel, copasi2sbmlmap);
1456 
1457  this->mpCopasiModel->setCompileFlag();
1458 
1459  if (this->mUnsupportedRuleFound)
1460  {
1462  }
1463 
1464  if (this->mRateRuleForSpeciesReferenceIgnored == true)
1465  {
1466  CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 94 , "Rate rule" , "Rate rule");
1467  }
1468 
1469  if (this->mEventAssignmentForSpeciesReferenceIgnored == true)
1470  {
1471  CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 94 , "Event assignment", "event assignment");
1472  }
1473 
1474  if (this->mConversionFactorFound == true)
1475  {
1477  }
1478 
1479  if (this->mEventPrioritiesIgnored == true)
1480  {
1482  }
1483 
1484  if (this->mInitialTriggerValues == true)
1485  {
1487  }
1488 
1489  if (this->mNonPersistentTriggerFound == true)
1490  {
1492  }
1493 
1494  if (mpImportHandler)
1495  {
1496  mpImportHandler->finishItem(hStep);
1497  mImportStep = 10;
1498 
1499  if (!mpImportHandler->progressItem(mhImportStep)) return NULL;
1500  }
1501 
1502  // unset the hasOnlySubstanceUnits flag on all such species
1503  std::map<Species*, Compartment*>::iterator it = this->mSubstanceOnlySpecies.begin();
1504  std::map<Species*, Compartment*>::iterator endIt = this->mSubstanceOnlySpecies.end();
1505 
1506  while (it != endIt)
1507  {
1508  it->first->setHasOnlySubstanceUnits(false);
1509  ++it;
1510  }
1511 
1512  setInitialValues(this->mpCopasiModel, copasi2sbmlmap);
1513  // evaluate and apply the initial expressions
1514  this->applyStoichiometricExpressions(copasi2sbmlmap, sbmlModel);
1515 
1516  // now we apply the conversion factors
1517  if (this->mLevel > 2)
1518  {
1519  this->applyConversionFactors();
1520  }
1521 
1522  this->removeUnusedFunctions(pTmpFunctionDB, copasi2sbmlmap);
1523 
1524  // remove the temporary avogadro parameter if one was created
1525  if (this->mAvogadroCreated == true)
1526  {
1527  const Parameter* pParameter = *this->mPotentialAvogadroNumbers.begin();
1528  ListOf* pList = sbmlModel->getListOfParameters();
1529  unsigned i, iMax = pList->size();
1530 
1531  for (i = 0; i < iMax; ++i)
1532  {
1533  if (pList->get(i)->getId() == pParameter->getId())
1534  {
1535  pList->remove(i);
1536  break;
1537  }
1538  }
1539  }
1540 
1541  delete pTmpFunctionDB;
1542 
1543  // create a warning if the delay function is used in the model
1544  if (this->mDelayFound)
1545  {
1547  }
1548 
1549  // give a warning that units on pure number as allowed in SBML L3 and above
1550  // are ignored by COPASI
1551  if (this->mLevel > 2 && this->mUnitOnNumberFound)
1552  {
1554  }
1555 
1557  return this->mpCopasiModel;
1558 }
1559 
1560 std::string isKnownCustomFunctionDefinition(const FunctionDefinition* sbmlFunction,
1561  const std::string& sNamespace,
1562  const std::string& elementName,
1563  const std::string& definition)
1564 {
1565  FunctionDefinition* current = const_cast<FunctionDefinition*>(sbmlFunction);
1566 
1567  if (current == NULL) return "";
1568 
1569  if (!current->isSetAnnotation()) return "";
1570 
1571  const XMLNode* element = current->getAnnotation();
1572 
1573  if (element == NULL) return "";
1574 
1575  for (unsigned int i = 0 ; i < element->getNumChildren(); ++i)
1576  {
1577  const XMLNode& annot = element->getChild(i);
1578 
1579  if (annot.getURI() == sNamespace &&
1580  annot.getName() == elementName &&
1581  annot.getAttrValue("definition") == definition)
1582  {
1583  return current->getId();
1584  }
1585  }
1586 
1587  return "";
1588 }
1589 
1590 bool addToKnownFunctionToMap(std::map<std::string, std::string>& map, const FunctionDefinition* sbmlFunction)
1591 {
1592  if (!sbmlFunction->isSetAnnotation())
1593  return false;
1594 
1595  std::string id = isKnownCustomFunctionDefinition(sbmlFunction,
1596  "http://sbml.org/annotations/symbols",
1597  "symbols",
1598  "http://en.wikipedia.org/wiki/Derivative");
1599 
1600  if (!id.empty())
1601  {
1602  map[id] = "RATE";
1603  return true;
1604  }
1605 
1606  id = isKnownCustomFunctionDefinition(sbmlFunction,
1607  "http://sbml.org/annotations/distribution",
1608  "distribution",
1609  "http://www.uncertml.org/distributions/normal");
1610 
1611  if (!id.empty())
1612  {
1613  map[id] = "RNORMAL";
1614  return true;
1615  }
1616 
1617  id = isKnownCustomFunctionDefinition(sbmlFunction,
1618  "http://sbml.org/annotations/distribution",
1619  "distribution",
1620  "http://www.uncertml.org/distributions/uniform");
1621 
1622  if (!id.empty())
1623  {
1624  map[id] = "RUNIFORM";
1625  return true;
1626  }
1627 
1628  id = isKnownCustomFunctionDefinition(sbmlFunction,
1629  "http://sbml.org/annotations/distribution",
1630  "distribution",
1631  "http://www.uncertml.org/distributions/gamma");
1632 
1633  if (!id.empty())
1634  {
1635  map[id] = "RGAMMA";
1636  return true;
1637  }
1638 
1639  id = isKnownCustomFunctionDefinition(sbmlFunction,
1640  "http://sbml.org/annotations/distribution",
1641  "distribution",
1642  "http://www.uncertml.org/distributions/poisson");
1643 
1644  if (!id.empty())
1645  {
1646  map[id] = "RPOISSON";
1647  return true;
1648  }
1649 
1650  return false;
1651 }
1652 
1653 int
1654 AstStrCmp(const void *s1, const void *s2)
1655 {
1656  const char* a = static_cast<const ASTNode *>(s1)->getName();
1657  const char* b = static_cast<const ASTNode *>(s2)->getName();
1658 
1659  if (a == NULL && b == NULL) return 0;
1660 
1661  if (a == NULL && b != NULL) return -1;
1662 
1663  if (a != NULL && b == NULL) return 1;
1664 
1665  return strcmp(
1666  a,
1667  b);
1668 }
1669 
1670 /**
1671  * This function checks the function definition for unused arguments (that would not
1672  * be properly displayed in the CopasiUI). If found, the function body will be replaced
1673  * with one including all arguments.
1674  */
1675 void ensureAllArgsAreBeingUsedInFunctionDefinition(const FunctionDefinition* sbmlFunction)
1676 {
1677  if (sbmlFunction == NULL || sbmlFunction->getNumArguments() == 0 || sbmlFunction->getBody() == NULL) return;
1678 
1679  // get all variables
1680  List *variables = sbmlFunction->getBody()->getListOfNodes(ASTNode_isName);
1681 
1682  // find unused ones
1683  std::vector<std::string> unused;
1684 
1685  for (unsigned int i = 0; i < sbmlFunction->getNumArguments(); ++i)
1686  {
1687  const ASTNode *arg = sbmlFunction->getArgument(i);
1688 
1689  if (variables->find(arg, AstStrCmp) == NULL)
1690  {
1691  if (arg->getName() != NULL)
1692  unused.push_back(arg->getName());
1693  }
1694  }
1695 
1696  // get rid of the list
1697  delete variables;
1698 
1699  // let us hope this is empty
1700  if (unused.size() == 0)
1701  return;
1702 
1703  // it is not, so modify the function definition to include all of them
1704  std::stringstream str;
1705  str << "lambda(";
1706 
1707  for (unsigned int i = 0; i < sbmlFunction->getNumArguments(); ++i)
1708  str << sbmlFunction->getArgument(i)->getName() << ", ";
1709 
1710  char* formula = SBML_formulaToString(sbmlFunction->getBody());
1711  str << formula;
1712 
1713  std::vector<std::string>::iterator it;
1714 
1715  for (it = unused.begin(); it != unused.end(); ++it)
1716  str << " + 0*" << *it;
1717 
1718  str << ")";
1719 
1720  // update the function definition
1721  const_cast<FunctionDefinition*>(sbmlFunction)->setMath(SBML_parseFormula(str.str().c_str()));
1722 
1723  // free the formula
1724  free(formula);
1725 }
1726 
1727 CFunction* SBMLImporter::createCFunctionFromFunctionDefinition(const FunctionDefinition* sbmlFunction, CFunctionDB* pTmpFunctionDB, Model* pSBMLModel, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
1728 {
1729 
1731 
1733 
1734  CFunction* pTmpFunction = this->createCFunctionFromFunctionTree(sbmlFunction, pSBMLModel, copasi2sbmlmap);
1735 
1736  if (pTmpFunction == NULL)
1737  {
1738  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 14, sbmlFunction->getId().c_str());
1739  return NULL;
1740  }
1741 
1742  std::string sbmlId = sbmlFunction->getId();
1743  pTmpFunction->setSBMLId(sbmlId);
1744  // check if the id is already taken by another function definition, maybe
1745  // from an earlier import, if this is the case, delete the id on the old
1746  // function definition
1747  // if we don't do this, two functions might have the same SBML id during
1748  // export which makes the exporter code so much more difficult
1749  size_t i, iMax = this->functionDB->loadedFunctions().size();
1750 
1751  for (i = 0; i < iMax; ++i)
1752  {
1753  CFunction* pFun = this->functionDB->loadedFunctions()[i];
1754 
1755  if (pFun->getSBMLId() == sbmlId)
1756  {
1757  pFun->setSBMLId("");
1758  }
1759  }
1760 
1761  std::string functionName = sbmlFunction->getName();
1762 
1763  if (functionName == "")
1764  {
1765  functionName = sbmlFunction->getId();
1766  }
1767 
1768  unsigned int counter = 1;
1769  std::ostringstream numberStream;
1770  std::string appendix = "";
1771  CFunction * pExistingFunction = NULL;
1772 
1773  while ((pExistingFunction = functionDB->findFunction(functionName + appendix)))
1774  {
1775  if (areEqualFunctions(pExistingFunction, pTmpFunction))
1776  {
1777  pdelete(pTmpFunction);
1778  pTmpFunction = pExistingFunction;
1779 
1780  break;
1781  }
1782 
1783  // We need to check whether the functions are identical.
1784  numberStream.str("");
1785  numberStream << "_" << counter;
1786  counter++;
1787  appendix = numberStream.str();
1788  }
1789 
1790  if (pTmpFunction != pExistingFunction)
1791  {
1792  pTmpFunction->setObjectName(functionName + appendix);
1793  functionDB->add(pTmpFunction, true);
1794  pTmpFunctionDB->add(pTmpFunction, false);
1795  }
1796 
1797  if (pTmpFunction->getType() == CEvaluationTree::UserDefined)
1798  {
1799  SBMLImporter::importMIRIAM(sbmlFunction, pTmpFunction);
1800  SBMLImporter::importNotes(pTmpFunction, sbmlFunction);
1801  }
1802 
1803  return pTmpFunction;
1804 }
1805 
1806 CFunction* SBMLImporter::createCFunctionFromFunctionTree(const FunctionDefinition* pSBMLFunction, Model* pSBMLModel, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
1807 {
1808  CFunction* pFun = NULL;
1809 
1810  if (!pSBMLFunction->isSetMath())
1811  return NULL;
1812 
1813  ConverterASTNode root(*pSBMLFunction->getMath());
1814 
1816  {
1817  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 85, pSBMLFunction->getId().c_str());
1818  }
1819 
1820  this->preprocessNode(&root, pSBMLModel, copasi2sbmlmap);
1821 
1822  if (root.getType() != AST_LAMBDA)
1823  {
1824  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 11, pSBMLFunction->getId().c_str());
1825  return NULL;
1826  }
1827 
1828  // get the number of children.
1829  // the first n-1 children are the parameters for the function
1830  // the last child is the actual function
1831  pFun = new CKinFunction();
1832  unsigned int i, iMax = root.getNumChildren() - 1;
1833  std::set<std::string> variableNames;
1834 
1835  for (i = 0; i < iMax; ++i)
1836  {
1837  ASTNode* pVarNode = root.getChild(i);
1838 
1839  if (pVarNode->getType() != AST_NAME)
1840  {
1841  delete pFun;
1842  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 12, pSBMLFunction->getId().c_str());
1843  }
1844 
1845  pFun->addVariable(pVarNode->getName());
1846  variableNames.insert(pVarNode->getName());
1847  }
1848 
1849  // now we check if the AST tree has a node that represents the time
1850  // object
1851  // find a unique name for the time variable
1852  std::ostringstream sstream;
1853  std::string timeVariableName = "time";
1854  unsigned int postfix = 1;
1855 
1856  while (variableNames.find(timeVariableName) != variableNames.end())
1857  {
1858  sstream.str("");
1859  sstream << "time_" << postfix;
1860  timeVariableName = sstream.str();
1861  ++postfix;
1862  }
1863 
1864  if (this->replaceTimeNodesInFunctionDefinition(root.getChild(iMax), timeVariableName))
1865  {
1866  // add another variable to the function
1867  ASTNode* pVarNode = new ASTNode(AST_NAME);
1868  pVarNode->setName(timeVariableName.c_str());
1869  ASTNode* pTmpNode = root.removeChild(iMax);
1870  root.addChild(pVarNode);
1871  root.addChild(pTmpNode);
1872  // increase iMax since we now have one more child
1873  ++iMax;
1874  pFun->addVariable(timeVariableName);
1875  this->mExplicitelyTimeDependentFunctionDefinitions.insert(pSBMLFunction->getId());
1876  }
1877 
1878  pFun->setTree(*root.getChild(iMax));
1880 
1881  // if the root node already is an object node, this has to be dealt with separately
1882  if (dynamic_cast<CEvaluationNodeObject*>(&(*treeIt)))
1883  {
1884  CEvaluationNodeVariable* pVariableNode = new CEvaluationNodeVariable(CEvaluationNodeVariable::ANY, (*treeIt).getData().substr(1, (*treeIt).getData().length() - 2));
1885  pFun->setRoot(pVariableNode);
1886  }
1887  else
1888  {
1889  while (treeIt != NULL)
1890  {
1891  if (dynamic_cast<CEvaluationNodeObject*>(&(*treeIt)))
1892  {
1893  CEvaluationNodeVariable* pVariableNode = new CEvaluationNodeVariable(CEvaluationNodeVariable::ANY, (*treeIt).getData().substr(1, (*treeIt).getData().length() - 2));
1894 
1895  if ((*treeIt).getParent())
1896  {
1897  (*treeIt).getParent()->addChild(pVariableNode, &(*treeIt));
1898  (*treeIt).getParent()->removeChild(&(*treeIt));
1899  }
1900 
1901  delete &(*treeIt);
1902  treeIt = pVariableNode;
1903  }
1904 
1905  ++treeIt;
1906  }
1907  }
1908 
1909  pFun->updateTree();
1910 
1911  if (!pFun->compile())
1912  {
1913  delete pFun;
1914  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 28, pSBMLFunction->getId().c_str());
1915  }
1916 
1917  if (pFun->getRoot() == NULL)
1918  {
1919  delete pFun;
1920  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 13, pSBMLFunction->getId().c_str());
1921  }
1922 
1923  return pFun;
1924 }
1925 
1926 CFunction* getFunctionForKey(CCopasiVectorN<CFunction> &functionDb, const std::string& key)
1927 {
1928  CFunction* pFunc = NULL;
1929  CCopasiVectorN<CFunction>::iterator it = functionDb.begin(), endit = functionDb.end();
1930 
1931  while (it != endit)
1932  {
1933  if ((*it)->getKey() == key)
1934  {
1935  pFunc = dynamic_cast<CFunction*>(*it);
1936  break;
1937  }
1938 
1939  ++it;
1940  }
1941 
1942  return pFunc;
1943 }
1944 
1945 /**
1946  * Creates and returns a COPASI CCompartment from the SBML Compartment
1947  * given as argument.
1948  */
1949 CCompartment*
1950 SBMLImporter::createCCompartmentFromCompartment(const Compartment* sbmlCompartment, CModel* copasiModel, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap, const Model* /*pSBMLModel*/)
1951 {
1952  if (sbmlCompartment->isSetUnits())
1953  {
1954  std::string cU = sbmlCompartment->getUnits();
1955  }
1956 
1957  unsigned int dimensionality = 3;
1958 
1959  if (sbmlCompartment->isSetSpatialDimensions())
1960  {
1961  dimensionality = sbmlCompartment->getSpatialDimensions();
1962 
1963  // starting with SBML Level 3, the spatial dimensions of a compartment can be
1964  // any rational number
1965  double dDimensionality = sbmlCompartment->getSpatialDimensions();
1966 
1967  if (this->mLevel > 2)
1968  {
1969  dDimensionality = sbmlCompartment->getSpatialDimensionsAsDouble();
1970  }
1971 
1972  // check if the unsigned int dimensionality corresponds to the double
1973  // dimensionality, otherwise give an error message
1974  // Actually if we wanted to be correct, we would have to check if the
1975  // part before the komma also matches and maybe have a separate error message if not
1976  dDimensionality -= dimensionality;
1977 
1978  if (fabs(dDimensionality) > 1e-9)
1979  {
1980  CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 89, sbmlCompartment->getId().c_str());
1981  dimensionality = 3;
1982  }
1983  }
1984  else
1985  {
1986  CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 90, sbmlCompartment->getId().c_str());
1987  dimensionality = 3;
1988  }
1989 
1990  if (dimensionality > 3)
1991  {
1993  "Reaction with id \"%s\" has dimensions of %d, this is not supported by COPASI. COPASI will assume that the compartment is three dimensional."
1994  , sbmlCompartment->getId().c_str(), dimensionality);
1995  dimensionality = 3;
1996  //fatalError();
1997  }
1998 
1999  std::string name = sbmlCompartment->getName();
2000 
2001  if (name == "")
2002  {
2003  name = sbmlCompartment->getId();
2004  }
2005 
2006  std::string appendix = "";
2007  unsigned int counter = 2;
2008  std::ostringstream numberStream;
2009 
2010  while (copasiModel->getCompartments().getIndex(name + appendix) != C_INVALID_INDEX)
2011  {
2012  numberStream.str("");
2013  numberStream << "_" << counter;
2014  counter++;
2015  appendix = numberStream.str();
2016  }
2017 
2018  double value;
2019  CCompartment* copasiCompartment = copasiModel->createCompartment(name + appendix, value);
2020 
2021  if (this->mLevel == 1)
2022  {
2023  copasiCompartment->setSBMLId(sbmlCompartment->getName());
2024  }
2025  else
2026  {
2027  copasiCompartment->setSBMLId(sbmlCompartment->getId());
2028  }
2029 
2030  // set the dimension of the compartment
2031  copasiCompartment->setDimensionality(dimensionality);
2032 
2033  //DebugFile << "Created Compartment: " << copasiCompartment->getObjectName() << std::endl;
2034  SBMLImporter::importMIRIAM(sbmlCompartment, copasiCompartment);
2035  SBMLImporter::importNotes(copasiCompartment, sbmlCompartment);
2036  copasi2sbmlmap[copasiCompartment] = const_cast<Compartment*>(sbmlCompartment);
2037  return copasiCompartment;
2038 }
2039 
2040 /**
2041  * Creates and returns a Copasi CMetab from the given SBML Species object.
2042  */
2043 CMetab*
2044 SBMLImporter::createCMetabFromSpecies(const Species* sbmlSpecies, CModel* copasiModel, CCompartment* copasiCompartment, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap, const Model* /*pSBMLModel*/)
2045 {
2046  if (sbmlSpecies->isSetSubstanceUnits())
2047  {
2048  std::string cU = sbmlSpecies->getSubstanceUnits();
2049  }
2050 
2051  std::map<CCopasiObject*, SBase*>::iterator it = copasi2sbmlmap.find(copasiCompartment);
2052 
2053  if (it == copasi2sbmlmap.end())
2054  {
2055  fatalError();
2056  }
2057 
2058  Compartment* pSBMLCompartment = (Compartment*)it->second;
2059 
2060  if (sbmlSpecies->getHasOnlySubstanceUnits() == true)
2061  {
2062  this->mSubstanceOnlySpecies.insert(std::make_pair(const_cast<Species*>(sbmlSpecies), pSBMLCompartment));
2063  }
2064 
2065  std::string name = sbmlSpecies->getName();
2066 
2067  if (name == "")
2068  {
2069  name = sbmlSpecies->getId();
2070  }
2071 
2072  std::string appendix = "";
2073  unsigned int counter = 2;
2074  std::ostringstream numberStream;
2075 
2076  while (copasiCompartment->getMetabolites().getIndex(name + appendix) != C_INVALID_INDEX)
2077  {
2078  numberStream.str("");
2079  numberStream << "_" << counter;
2080  counter++;
2081  appendix = numberStream.str();
2082  }
2083 
2084  CMetab* copasiMetabolite = copasiModel->createMetabolite(name + appendix, copasiCompartment->getObjectName());
2085 
2086  if (copasiMetabolite == NULL)
2087  {
2088  //DebugFile << "Could not create Copasi species." << std::endl;
2089  fatalError();
2090  }
2091 
2092  if (sbmlSpecies->getConstant() || sbmlSpecies->getBoundaryCondition())
2093  {
2094  copasiMetabolite->setStatus(CModelEntity::FIXED);
2095  }
2096  else
2097  {
2098  copasiMetabolite->setStatus(CModelEntity::REACTIONS);
2099  }
2100 
2101  // also check if the compartment has a spatialSize of 0 because this also implies hasOnlySubstanceUnits for the species in this compartment
2102  // In Level 3 files this no longer seems to be the case however, so we only
2103  // do this for L1 and L2 files
2104  if (pSBMLCompartment->isSetSpatialDimensions() && pSBMLCompartment->getSpatialDimensions() == 0 && this->mLevel < 3)
2105  {
2106  this->mSubstanceOnlySpecies.insert(std::make_pair(const_cast<Species*>(sbmlSpecies), pSBMLCompartment));
2107  }
2108 
2109  //DebugFile << "Created species: " << copasiMetabolite->getObjectName() << std::endl;
2110  copasi2sbmlmap[copasiMetabolite] = const_cast<Species*>(sbmlSpecies);
2111 
2112  if (this->mLevel == 1)
2113  {
2114  copasiMetabolite->setSBMLId(sbmlSpecies->getName());
2115  }
2116  else
2117  {
2118  copasiMetabolite->setSBMLId(sbmlSpecies->getId());
2119  }
2120 
2121  SBMLImporter::importMIRIAM(sbmlSpecies, copasiMetabolite);
2122  SBMLImporter::importNotes(copasiMetabolite, sbmlSpecies);
2123  return copasiMetabolite;
2124 }
2125 
2127 {
2128  size_t i, iMax = db.size();
2129 
2130  for (i = 0; i < iMax; ++i)
2131  {
2132  CFunction* pFun = (db[i]);
2133 
2134  if (*pFun == *func)
2135  return pFun;
2136  }
2137 
2138  return NULL;
2139 }
2140 
2141 /**
2142  * Creates and returns a COPASI CReaction object from the given SBML
2143  * Reaction object.
2144  */
2145 CReaction*
2146 SBMLImporter::createCReactionFromReaction(Reaction* sbmlReaction, Model* pSBMLModel, CModel* copasiModel, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap, CFunctionDB* pTmpFunctionDB)
2147 {
2148  if (sbmlReaction == NULL)
2149  {
2150  fatalError();
2151  }
2152 
2153  std::string name = sbmlReaction->getName();
2154 
2155  if (name == "")
2156  {
2157  name = sbmlReaction->getId();
2158  }
2159 
2160  /* Check if the name of the reaction is unique. */
2161  std::string appendix = "";
2162  unsigned int counter = 2;
2163  std::ostringstream numberStream;
2164 
2165  while (copasiModel->getReactions().getIndex(name + appendix) != C_INVALID_INDEX)
2166  {
2167  numberStream.str("");
2168  numberStream << "_" << counter;
2169  counter++;
2170  appendix = numberStream.str();
2171  }
2172 
2173  /* create a new reaction with the unique name */
2174  CReaction* copasiReaction = copasiModel->createReaction(name + appendix);
2175 
2176  if (copasiReaction == NULL)
2177  {
2178  fatalError();
2179  }
2180 
2181  copasiReaction->setReversible(sbmlReaction->getReversible());
2182 
2183  if (this->mLevel == 1)
2184  {
2185  copasiReaction->setSBMLId(sbmlReaction->getName());
2186  }
2187  else
2188  {
2189  copasiReaction->setSBMLId(sbmlReaction->getId());
2190  }
2191 
2192  copasi2sbmlmap[copasiReaction] = const_cast<Reaction*>(sbmlReaction);
2193 
2194  if (sbmlReaction->isSetFast() && sbmlReaction->getFast() == true)
2195  {
2196  const_cast<Reaction*>(sbmlReaction)->setFast(false);
2197  this->mFastReactions.insert(sbmlReaction->getId());
2198  }
2199 
2200  // store if the reaction involves species that need a conversion factor
2201  bool mConversionFactorNeeded = false;
2202 
2203  /* Add all substrates to the reaction */
2204  unsigned int num = sbmlReaction->getNumReactants();
2205  bool singleCompartment = true;
2206  const CCompartment* compartment = NULL;
2207  bool hasOnlySubstanceUnitPresent = false;
2208  bool ignoreMassAction = false;
2209 
2210  for (counter = 0; counter < num; counter++)
2211  {
2212  const SpeciesReference* sr = sbmlReaction->getReactant(counter);
2213 
2214  if (sr == NULL)
2215  {
2216  delete copasiReaction;
2217  fatalError();
2218  }
2219 
2220  C_FLOAT64 stoi = 1.0;
2221 
2222  if (this->mLevel < 3)
2223  {
2224  // We may not use Mass Action kinetics if we do not know the actual stoichiometry.
2225  if (sr->isSetStoichiometryMath())
2226  {
2227  ignoreMassAction = true;
2228  }
2229  else
2230  {
2231  stoi = sr->getStoichiometry() / sr->getDenominator();
2232  }
2233  }
2234  else
2235  {
2236  // We may not use Mass Action kinetics if we do not know the actual stoichiometry.
2237  if (sr->isSetId() &&
2238  (pSBMLModel->getInitialAssignment(sr->getId()) ||
2239  pSBMLModel->getRule(sr->getId())))
2240  {
2241  ignoreMassAction = true;
2242  }
2243 
2244  if (sr->isSetStoichiometry())
2245  {
2246  stoi = sr->getStoichiometry();
2247  }
2248  }
2249 
2250  std::map<std::string, CMetab*>::iterator pos;
2251  pos = this->speciesMap.find(sr->getSpecies());
2252 
2253  // check if we have a conversion factor
2254  if (pos == this->speciesMap.end())
2255  {
2256  delete copasiReaction;
2257  fatalError();
2258  }
2259  else
2260  {
2261  // check if there is a conversion factor on the species
2262  if (this->mpModelConversionFactor != NULL || this->mSpeciesConversionParameterMap.find(pos->second->getSBMLId()) != this->mSpeciesConversionParameterMap.end())
2263  {
2264  mConversionFactorNeeded = true;
2265  }
2266  }
2267 
2268  std::map<CCopasiObject*, SBase*>::const_iterator spos = copasi2sbmlmap.find(pos->second);
2269  assert(spos != copasi2sbmlmap.end());
2270  Species* pSBMLSpecies = dynamic_cast<Species*>(spos->second);
2271  assert(pSBMLSpecies != NULL);
2272  hasOnlySubstanceUnitPresent = (hasOnlySubstanceUnitPresent | (pSBMLSpecies->getHasOnlySubstanceUnits() == true));
2273 
2274  if (compartment == NULL)
2275  {
2276  compartment = pos->second->getCompartment();
2277  }
2278  else
2279  {
2280  if (singleCompartment && compartment != pos->second->getCompartment())
2281  {
2282  singleCompartment = false;
2283  }
2284  }
2285 
2286  copasiReaction->addSubstrate(pos->second->getKey(), stoi);
2287 
2288  // we need to store the id of the species reference if it is set because SBML Level 3 allows
2289  // references to species references and if we want to support his, we need the id to import
2290  // expressions that reference a species reference
2291  if (this->mLevel > 2)
2292  {
2293  CCopasiVector<CChemEqElement>::const_iterator it = copasiReaction->getChemEq().getSubstrates().begin(), endit = copasiReaction->getChemEq().getSubstrates().end();
2294  CChemEqElement* pElement = NULL;
2295 
2296  while (it != endit)
2297  {
2298  if ((*it)->getMetaboliteKey() == pos->second->getKey())
2299  {
2300  pElement = const_cast<CChemEqElement*>(*it);
2301  break;
2302  }
2303 
2304  ++it;
2305  }
2306 
2307  assert(pElement != NULL);
2308 
2309  if (pElement != NULL)
2310  {
2311  copasi2sbmlmap[pElement] = const_cast<SpeciesReference*>(sr);
2312 
2313  this->mChemEqElementSpeciesIdMap[pElement] = std::pair<std::string, CChemEq::MetaboliteRole>(sr->getSpecies(), CChemEq::SUBSTRATE);
2314  }
2315  }
2316 
2317  // find the CChemEqElement that belongs to the added substrate
2318  if (this->mLevel < 3 && sr->isSetStoichiometryMath())
2319  {
2320  CChemEq& chemEq = copasiReaction->getChemEq();
2323  CChemEqElement* pChemEqElement = NULL;
2324 
2325  while (it != end)
2326  {
2327  if ((*it)->getMetabolite() == pos->second)
2328  {
2329  pChemEqElement = (*it);
2330  break;
2331  }
2332 
2333  ++it;
2334  }
2335 
2336  assert(pChemEqElement != NULL);
2337  mStoichiometricExpressionMap.insert(std::make_pair(sr->getStoichiometryMath()->getMath(), pChemEqElement));
2338  }
2339  }
2340 
2341  /* Add all products to the reaction */
2342  num = sbmlReaction->getNumProducts();
2343 
2344  for (counter = 0; counter < num; counter++)
2345  {
2346  const SpeciesReference* sr = sbmlReaction->getProduct(counter);
2347 
2348  if (sr == NULL)
2349  {
2350  delete copasiReaction;
2351  fatalError();
2352  }
2353 
2354  C_FLOAT64 stoi = 1.0;
2355 
2356  if (this->mLevel < 3)
2357  {
2358  if (sr->isSetStoichiometryMath())
2359  {
2360  ignoreMassAction = true;
2361  }
2362  else
2363  {
2364  stoi = sr->getStoichiometry() / sr->getDenominator();
2365  }
2366  }
2367  else
2368  {
2369  // We may not use Mass Action kinetics if we do not know the actual stoichiometry.
2370  if (sr->isSetId() &&
2371  (pSBMLModel->getInitialAssignment(sr->getId()) ||
2372  pSBMLModel->getRule(sr->getId())))
2373  {
2374  ignoreMassAction = true;
2375  }
2376 
2377  if (sr->isSetStoichiometry())
2378  {
2379  stoi = sr->getStoichiometry();
2380  }
2381  }
2382 
2383  std::map<std::string, CMetab*>::iterator pos;
2384  pos = this->speciesMap.find(sr->getSpecies());
2385 
2386  if (pos == this->speciesMap.end())
2387  {
2388  delete copasiReaction;
2389  fatalError();
2390  }
2391  else
2392  {
2393  // check if there is a conversion factor on the species
2394  if (this->mpModelConversionFactor != NULL || this->mSpeciesConversionParameterMap.find(pos->second->getSBMLId()) != this->mSpeciesConversionParameterMap.end())
2395  {
2396  mConversionFactorNeeded = true;
2397  }
2398  }
2399 
2400  std::map<CCopasiObject*, SBase*>::const_iterator spos = copasi2sbmlmap.find(pos->second);
2401  assert(spos != copasi2sbmlmap.end());
2402  Species* pSBMLSpecies = dynamic_cast<Species*>(spos->second);
2403  assert(pSBMLSpecies != NULL);
2404  hasOnlySubstanceUnitPresent = (hasOnlySubstanceUnitPresent | (pSBMLSpecies->getHasOnlySubstanceUnits() == true));
2405 
2406  if (compartment == NULL)
2407  {
2408  compartment = pos->second->getCompartment();
2409  }
2410  else
2411  {
2412  if (singleCompartment && compartment != pos->second->getCompartment())
2413  {
2414  singleCompartment = false;
2415  }
2416  }
2417 
2418  copasiReaction->addProduct(pos->second->getKey(), stoi);
2419 
2420  // we need to store the id of the species reference if it is set because SBML Level 3 allows
2421  // references to species references and if we want to support his, we need the id to import
2422  // expressions that reference a species reference
2423  if (this->mLevel > 2)
2424  {
2425  CCopasiVector<CChemEqElement>::const_iterator it = copasiReaction->getChemEq().getProducts().begin(), endit = copasiReaction->getChemEq().getProducts().end();
2426  CChemEqElement* pElement = NULL;
2427 
2428  while (it != endit)
2429  {
2430  if ((*it)->getMetaboliteKey() == pos->second->getKey())
2431  {
2432  pElement = const_cast<CChemEqElement*>(*it);
2433  break;
2434  }
2435 
2436  ++it;
2437  }
2438 
2439  assert(pElement != NULL);
2440 
2441  if (pElement != NULL)
2442  {
2443  copasi2sbmlmap[pElement] = const_cast<SpeciesReference*>(sr);
2444  this->mChemEqElementSpeciesIdMap[pElement] = std::pair<std::string, CChemEq::MetaboliteRole>(sr->getSpecies(), CChemEq::PRODUCT);
2445  }
2446  }
2447 
2448  if (sr->isSetStoichiometryMath())
2449  {
2450  CChemEq& chemEq = copasiReaction->getChemEq();
2453  CChemEqElement* pChemEqElement = NULL;
2454 
2455  while (it != end)
2456  {
2457  if ((*it)->getMetabolite() == pos->second)
2458  {
2459  pChemEqElement = (*it);
2460  break;
2461  }
2462 
2463  ++it;
2464  }
2465 
2466  assert(pChemEqElement != NULL);
2467  mStoichiometricExpressionMap.insert(std::make_pair(sr->getStoichiometryMath()->getMath(), pChemEqElement));
2468  }
2469  }
2470 
2471  /* Add all modifiers to the reaction */
2472  num = sbmlReaction->getNumModifiers();
2473 
2474  for (counter = 0; counter < num; counter++)
2475  {
2476  const ModifierSpeciesReference* sr = sbmlReaction->getModifier(counter);
2477 
2478  if (sr == NULL)
2479  {
2480  delete copasiReaction;
2481  fatalError();
2482  }
2483 
2484  std::map<std::string, CMetab*>::iterator pos;
2485  pos = this->speciesMap.find(sr->getSpecies());
2486 
2487  if (pos == this->speciesMap.end())
2488  {
2489  delete copasiReaction;
2490  fatalError();
2491  }
2492 
2493  std::map<CCopasiObject*, SBase*>::const_iterator spos = copasi2sbmlmap.find(pos->second);
2494  assert(spos != copasi2sbmlmap.end());
2495  Species* pSBMLSpecies = dynamic_cast<Species*>(spos->second);
2496  assert(pSBMLSpecies != NULL);
2497  hasOnlySubstanceUnitPresent = (hasOnlySubstanceUnitPresent | (pSBMLSpecies->getHasOnlySubstanceUnits() == true));
2498  copasiReaction->addModifier(pos->second->getKey());
2499 
2500  // we need to store the id of the species reference if it is set because SBML Level 3 allows
2501  // references to species references and if we want to support his, we need the id to import
2502  // expressions that reference a species reference
2503  if (this->mLevel > 2)
2504  {
2505 
2506  CCopasiVector<CChemEqElement>::const_iterator it = copasiReaction->getChemEq().getModifiers().begin(), endit = copasiReaction->getChemEq().getModifiers().end();
2507  CChemEqElement* pElement = NULL;
2508 
2509  while (it != endit)
2510  {
2511  if ((*it)->getMetaboliteKey() == pos->second->getKey())
2512  {
2513  pElement = const_cast<CChemEqElement*>(*it);
2514  break;
2515  }
2516 
2517  ++it;
2518  }
2519 
2520  assert(pElement != NULL);
2521 
2522  if (pElement != NULL)
2523  {
2524  copasi2sbmlmap[pElement] = const_cast<ModifierSpeciesReference*>(sr);
2525  }
2526  }
2527  }
2528 
2529  /* in the newly created CFunction set the types for all parameters and
2530  * either a mapping or a value
2531  */
2532  const KineticLaw* kLaw = sbmlReaction->getKineticLaw();
2533 
2534  if (kLaw != NULL)
2535  {
2536  const ListOfParameters* pParamList = NULL;
2537 
2538  if (this->mLevel > 2)
2539  {
2540  pParamList = kLaw->getListOfLocalParameters();
2541  }
2542  else
2543  {
2544  pParamList = kLaw->getListOfParameters();
2545  }
2546 
2547  for (counter = 0; counter < pParamList->size(); ++counter)
2548  {
2549  const Parameter* pSBMLParameter = pParamList->get(counter);
2550  std::string id;
2551 
2552  if (this->mLevel == 1)
2553  {
2554  id = pSBMLParameter->getName();
2555  }
2556  else
2557  {
2558  id = pSBMLParameter->getId();
2559  }
2560 
2561  double value;
2562 
2563  if (pSBMLParameter->isSetValue() && pSBMLParameter->getValue() == pSBMLParameter->getValue()) // make sure it is not set to NaN
2564  {
2565  value = pSBMLParameter->getValue();
2566  }
2567  else
2568  {
2569  // Set value to NaN and create a warning if it is the first time
2570  // this happend
2571  value = std::numeric_limits<C_FLOAT64>::quiet_NaN();
2572 
2573  if (!this->mIncompleteModel)
2574  {
2575  this->mIncompleteModel = true;
2576  CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 42, pSBMLParameter->getId().c_str());
2577  }
2578  }
2579 
2580  copasiReaction->getParameters().addParameter(id, CCopasiParameter::DOUBLE, value);
2581  }
2582 
2583  const ASTNode* kLawMath = kLaw->getMath();
2584 
2585  if (kLawMath == NULL || kLawMath->getType() == AST_UNKNOWN)
2586  {
2587  copasiReaction->setFunction(NULL);
2588  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 56, sbmlReaction->getId().c_str());
2589  }
2590  else
2591  {
2592  // check for references to species references in the expression because we don't support them yet
2593  if (!SBMLImporter::findIdInASTTree(kLawMath, this->mSBMLSpeciesReferenceIds).empty())
2594  {
2596  }
2597 
2598  ConverterASTNode* node = new ConverterASTNode(*kLawMath);
2599  this->preprocessNode(node, pSBMLModel, copasi2sbmlmap, sbmlReaction);
2600 
2601  if (node == NULL)
2602  {
2603  delete copasiReaction;
2604  fatalError();
2605  }
2606 
2607  /* if it is a single compartment reaction, we have to divide the whole kinetic
2608  ** equation by the compartment because COPASI expects
2609  ** kinetic laws that specify concentration/time for single compartment
2610  ** reactions.
2611  */
2612  if (singleCompartment)
2613  {
2614 
2615  if (compartment != NULL)
2616  {
2617  // only divide if it is not a 0-dimensional compartment
2618  if (compartment->getDimensionality() != 0)
2619  {
2620  // check if the root node is a multiplication node and it's first child
2621  // is a compartment node, if so, drop those two and make the second child
2622  // new new root node
2623  ConverterASTNode* tmpNode1 = this->isMultipliedByVolume(node, compartment->getSBMLId());
2624 
2625  if (tmpNode1)
2626  {
2627  delete node;
2628  node = tmpNode1;
2629 
2630  if (node->getType() == AST_DIVIDE && node->getNumChildren() != 2)
2631  {
2632  delete tmpNode1;
2633  fatalError();
2634  }
2635  }
2636  else
2637  {
2638  tmpNode1 = new ConverterASTNode();
2639  tmpNode1->setType(AST_DIVIDE);
2640  tmpNode1->addChild(node);
2641  ConverterASTNode* tmpNode2 = new ConverterASTNode();
2642  tmpNode2->setType(AST_NAME);
2643  tmpNode2->setName(compartment->getSBMLId().c_str());
2644  tmpNode1->addChild(tmpNode2);
2645  node = tmpNode1;
2646  std::map<CCopasiObject*, SBase*>::const_iterator pos = copasi2sbmlmap.find(const_cast<CCompartment*>(compartment));
2647  assert(pos != copasi2sbmlmap.end());
2648  Compartment* pSBMLCompartment = dynamic_cast<Compartment*>(pos->second);
2649  assert(pSBMLCompartment != NULL);
2650 
2651  if (!hasOnlySubstanceUnitPresent && ((this->mLevel == 1 && pSBMLCompartment->isSetVolume()) || (this->mLevel >= 2 && pSBMLCompartment->isSetSize())) && pSBMLCompartment->getSize() == 1.0)
2652  {
2653  // we have to check if all species used in the reaction
2654  // have the hasOnlySubstance flag set
2655 
2656  if (node->getChild(0)->getType() == AST_FUNCTION && (!this->containsVolume(node->getChild(0), compartment->getSBMLId())))
2657  {
2658  // add the id of the reaction to the set so that we can create an error message later.
2659  this->mDivisionByCompartmentReactions.insert(sbmlReaction->getId());
2660  }
2661  }
2662  }
2663  }
2664  }
2665  else
2666  {
2667  delete node;
2668  delete copasiReaction;
2669  fatalError();
2670  }
2671  }
2672 
2673  /* Create a new user defined CKinFunction */
2674  if (!sbmlId2CopasiCN(node, copasi2sbmlmap, copasiReaction->getParameters()))
2675  {
2676  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 27, copasiReaction->getObjectName().c_str());
2677  }
2678 
2679  CEvaluationNode* pExpressionTreeRoot = CEvaluationTree::fromAST(node);
2680  delete node;
2681  node = NULL;
2682 
2683  if (pExpressionTreeRoot)
2684  {
2685  CExpression KineticLawExpression;
2686  KineticLawExpression.setRoot(pExpressionTreeRoot);
2687 
2688  // check if the expression is constant flux
2689  CCopasiObject* pParamObject = SBMLImporter::isConstantFlux(pExpressionTreeRoot, copasiModel, pTmpFunctionDB);
2690 
2691  if (pParamObject != NULL)
2692  {
2693  // assert that the object really is a local or global
2694  // parameter
2695  assert(dynamic_cast<const CCopasiParameter*>(pParamObject) || dynamic_cast<const CModelValue*>(pParamObject));
2696  std::string functionName;
2697 
2698  if (copasiReaction->isReversible())
2699  {
2700  // set the function to Constant flux (reversible)
2701  functionName = "Constant flux (reversible)";
2702  }
2703  else
2704  {
2705  // set the function to Constant flux (irreversible)
2706  functionName = "Constant flux (irreversible)";
2707  }
2708 
2709  CFunction* pCFFun = dynamic_cast<CFunction*>(this->functionDB->findFunction(functionName));
2710  assert(pCFFun != NULL);
2711  CEvaluationNodeCall* pCallNode = NULL;
2712 
2713  if (CEvaluationNode::type(pExpressionTreeRoot->getType()) == CEvaluationNode::OBJECT)
2714  {
2715  pCallNode = new CEvaluationNodeCall(CEvaluationNodeCall::EXPRESSION, "dummy_call");
2716  // add the parameter
2717  pCallNode->addChild(pExpressionTreeRoot->copyBranch());
2718  }
2719  else
2720  {
2721  pCallNode = dynamic_cast<CEvaluationNodeCall*>(pExpressionTreeRoot->copyBranch());
2722  assert(pCallNode != NULL);
2723  }
2724 
2725  if (pParamObject->getObjectType() == "Parameter")
2726  {
2727  pParamObject->setObjectName("v");
2728  dynamic_cast<CEvaluationNode*>(pCallNode->getChild())->setData("<" + pParamObject->getCN() + ">");
2729  }
2730 
2731  copasiReaction->setFunction(pCFFun);
2732  // map the parameter
2733  this->doMapping(copasiReaction, pCallNode);
2734  delete pCallNode;
2735  }
2736  else
2737  {
2738  // check if the root node is a simple function call
2739  if (this->isSimpleFunctionCall(pExpressionTreeRoot))
2740  {
2741  // if yes, we check if it corresponds to an already existing function
2742  std::string functionName = pExpressionTreeRoot->getData();
2743  CFunction* pImportedFunction = dynamic_cast<CFunction*>(functionDB->findFunction(functionName));
2744  assert(pImportedFunction);
2745  std::vector<CEvaluationNodeObject*>* v = NULL;
2746 
2747  // only check for mass action if there is no conversion factor involved
2748  // for any of the species involved in the reaction (substrates and products)
2749  if (!mConversionFactorNeeded && !ignoreMassAction)
2750  {
2751  v = this->isMassAction(pImportedFunction, copasiReaction->getChemEq(), static_cast<const CEvaluationNodeCall*>(pExpressionTreeRoot));
2752 
2753  if (!v)
2754  {
2755  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 27, copasiReaction->getObjectName().c_str());
2756  }
2757  }
2758 
2759  if (v && !v->empty())
2760  {
2761  CFunction* pFun = NULL;
2762 
2763  if (copasiReaction->isReversible())
2764  {
2765  pFun = static_cast<CFunction*>(this->functionDB->findFunction("Mass action (reversible)"));
2766  }
2767  else
2768  {
2769  pFun = static_cast<CFunction*>(this->functionDB->findFunction("Mass action (irreversible)"));
2770  }
2771 
2772  if (!pFun)
2773  {
2774  fatalError();
2775  }
2776 
2777  // do the mapping
2779  size_t i, iMax = v->size();
2780 
2781  for (i = 0; i < iMax; ++i)
2782  {
2783  pCallNode->addChild((*v)[i]);
2784  }
2785 
2786  this->renameMassActionParameters(pCallNode);
2787  copasiReaction->setFunction(pFun);
2788  this->doMapping(copasiReaction, pCallNode);
2789  delete pCallNode;
2790  }
2791  else
2792  {
2793  // find corresponding reaction does *not* return the function if it is already in the function db
2794  // we better change this to return this instance.
2795  CFunction* pExistingFunction = findCorrespondingFunction(&KineticLawExpression, copasiReaction);
2796 
2797  // if it does, we set the existing function for this reaction
2798  if (pExistingFunction)
2799  {
2800  copasiReaction->setFunction(pExistingFunction);
2801  // do the mapping
2802  doMapping(copasiReaction, dynamic_cast<const CEvaluationNodeCall*>(pExpressionTreeRoot));
2803  }
2804  // else we take the function from the pTmpFunctionDB, copy it and set the usage correctly
2805  else
2806  {
2807  // replace the variable nodes in pImportedFunction with nodes from
2808  std::map<std::string, std::string > arguments;
2809  const CFunctionParameters& funParams = pImportedFunction->getVariables();
2810  const CEvaluationNode* pTmpNode = static_cast<const CEvaluationNode*>(pExpressionTreeRoot->getChild());
2811  size_t i, iMax = funParams.size();
2812 
2813  for (i = 0; (i < iMax) && pTmpNode; ++i)
2814  {
2815  if (!(pTmpNode->getType() == CEvaluationNode::OBJECT)) fatalError();
2816 
2817  arguments[funParams[i]->getObjectName()] = pTmpNode->getData().substr(1, pTmpNode->getData().length() - 2);
2818  pTmpNode = static_cast<const CEvaluationNode*>(pTmpNode->getSibling());
2819  }
2820 
2821  assert((i == iMax) && pTmpNode == NULL);
2822 
2823  CEvaluationNode* pTmpExpression = variables2objects(pImportedFunction->getRoot()->copyBranch(), arguments);
2824 
2825  CExpression TmpTree2;
2826  TmpTree2.setRoot(pTmpExpression);
2827 
2828  // code to fix bug 1874
2829  // since this is a user defined function, we want to retain the original name
2830  // next setFunctionFromExpressionTree will take this name if it is not 'Expression'
2831  // (the default)
2832  TmpTree2.setObjectName(functionName);
2833 
2834  CFunction * pNewFunction = copasiReaction->setFunctionFromExpressionTree(TmpTree2, copasi2sbmlmap, this->functionDB);
2835 
2836  if (pNewFunction != NULL &&
2837  pNewFunction->getType() == CEvaluationTree::UserDefined)
2838  {
2839  // code to fix Bug 1015
2840  if (!pNewFunction->isSuitable(copasiReaction->getChemEq().getSubstrates().size(),
2841  copasiReaction->getChemEq().getProducts().size(),
2842  copasiReaction->isReversible() ? TriTrue : TriFalse))
2843  {
2844  pNewFunction->setReversible(TriUnspecified);
2845  }
2846 
2847  pTmpFunctionDB->add(pNewFunction, false);
2848  mUsedFunctions.insert(pNewFunction->getObjectName());
2849  }
2850  }
2851  }
2852 
2853  pdelete(v);
2854  }
2855  else
2856  {
2857  std::vector<CEvaluationNodeObject*>* v = NULL;
2858 
2859  // only check for mass action if there is no conversion factor involved
2860  // for any of the species involved in the reaction (substrates and products)
2861  if (!mConversionFactorNeeded && !ignoreMassAction)
2862  {
2863  v = this->isMassAction(&KineticLawExpression, copasiReaction->getChemEq());
2864 
2865  if (!v)
2866  {
2867  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 27, copasiReaction->getObjectName().c_str());
2868  }
2869  }
2870 
2871  if (v && !v->empty())
2872  {
2873  CFunction* pFun = NULL;
2874 
2875  if (copasiReaction->isReversible())
2876  {
2877  pFun = static_cast<CFunction*>(this->functionDB->findFunction("Mass action (reversible)"));
2878  }
2879  else
2880  {
2881  pFun = static_cast<CFunction*>(this->functionDB->findFunction("Mass action (irreversible)"));
2882  }
2883 
2884  if (!pFun)
2885  {
2886  fatalError();
2887  }
2888 
2889  // do the mapping
2891  size_t i, iMax = v->size();
2892 
2893  for (i = 0; i < iMax; ++i)
2894  {
2895  pCallNode->addChild((*v)[i]);
2896  }
2897 
2898  // rename the function parameters to k1 and k2
2899  this->renameMassActionParameters(pCallNode);
2900  copasiReaction->setFunction(pFun);
2901  this->doMapping(copasiReaction, pCallNode);
2902  delete pCallNode;
2903  }
2904  else
2905  {
2906  CFunction* pNonconstFun = copasiReaction->setFunctionFromExpressionTree(KineticLawExpression, copasi2sbmlmap, this->functionDB);
2907 
2908  if (pNonconstFun != NULL &&
2909  pNonconstFun->getType() == CEvaluationTree::UserDefined)
2910  {
2911  // code to fix Bug 1015
2912  if (!pNonconstFun->isSuitable(copasiReaction->getChemEq().getSubstrates().size(),
2913  copasiReaction->getChemEq().getProducts().size(),
2914  copasiReaction->isReversible() ? TriTrue : TriFalse))
2915  {
2916  pNonconstFun->setReversible(TriUnspecified);
2917  }
2918 
2919  pTmpFunctionDB->add(pNonconstFun, false);
2920  mUsedFunctions.insert(pNonconstFun->getObjectName());
2921  }
2922  }
2923 
2924  pdelete(v);
2925  }
2926  }
2927  }
2928  else
2929  {
2930  // error message
2931  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 8, copasiReaction->getObjectName().c_str());
2932  }
2933  }
2934  }
2935  else
2936  {
2937  /* if no KineticLaw was defined for the reaction. */
2938  copasiReaction->setFunction(NULL);
2939  }
2940 
2941  //DebugFile << "Created reaction: " << copasiReaction->getObjectName() << std::endl;
2942  SBMLImporter::importMIRIAM(sbmlReaction, copasiReaction);
2943  SBMLImporter::importNotes(copasiReaction, sbmlReaction);
2944 
2945  return copasiReaction;
2946 }
2947 
2948 /**
2949  * Creates a map of each parameter of the function definition and its
2950  * corresponding parameter in the function call.
2951  */
2952 std::map<std::string, ASTNode*>
2953 SBMLImporter::createBVarMap(const ASTNode* uDefFunction, const ASTNode* function)
2954 {
2955  /* the first n-1 children, where n is the number of children, of a function definition ASTnode are the
2956  * arguments to the function. These correspond to the m=n-1 children of the
2957  * function call.
2958  */
2959  if (uDefFunction->getNumChildren() != function->getNumChildren() + 1)
2960  {
2961  std::string functionName = uDefFunction->getName();
2962  fatalError();
2963  }
2964 
2965  std::map<std::string, ASTNode*> varMap;
2966  unsigned int counter;
2967 
2968  for (counter = 0; counter < uDefFunction->getNumChildren() - 1; counter++)
2969  {
2970  varMap[uDefFunction->getChild(counter)->getName()] = function->getChild(counter);
2971  }
2972 
2973  return varMap;
2974 }
2975 
2976 /**
2977  * Returns the user defined SBML function definition that belongs to the given
2978  * name, or NULL if none can be found.
2979  */
2980 const FunctionDefinition*
2981 SBMLImporter::getFunctionDefinitionForName(const std::string name, const Model* sbmlModel)
2982 {
2983  const FunctionDefinition* fDef = NULL;
2984  unsigned int counter;
2985 
2986  for (counter = 0; counter < sbmlModel->getNumFunctionDefinitions(); counter++)
2987  {
2988  std::string functionName = sbmlModel->getFunctionDefinition(counter)->getName();
2989 
2990  if (sbmlModel->getFunctionDefinition(counter)->isSetId())
2991  {
2992  functionName = sbmlModel->getFunctionDefinition(counter)->getId();
2993  }
2994 
2995  if (functionName == name)
2996  {
2997  fDef = sbmlModel->getFunctionDefinition(counter);
2998  break;
2999  }
3000  }
3001 
3002  return fDef;
3003 }
3004 
3005 #ifdef XXXX
3006 /**
3007  * Replaces the variables in a function definition with the actual function
3008  * parameters that were used when the function was called. The function returns
3009  * a pointer to the ConverterAST node with the replaced variables.
3010  */
3012 SBMLImporter::replaceBvars(const ASTNode* node, std::map<std::string, ASTNode*> bvarMap)
3013 {
3014  ConverterASTNode* newNode = NULL;
3015 
3016  if (node->isName())
3017  {
3018  /* check if name matches any in bvarMap */
3019  /* if yes, replace node with node in bvarMap */
3020  /* node needs to be set to be a deep copy of the replacement */
3021  if (bvarMap.find(node->getName()) != bvarMap.end())
3022  {
3023  newNode = new ConverterASTNode(*bvarMap[node->getName()]);
3024  }
3025  }
3026  else
3027  {
3028  newNode = new ConverterASTNode(*node);
3029  newNode->setChildren(new List());
3030  unsigned int counter;
3031 
3032  for (counter = 0; counter < node->getNumChildren(); counter++)
3033  {
3034  newNode->addChild(this->replaceBvars(node->getChild(counter), bvarMap));
3035  }
3036  }
3037 
3038  return newNode;
3039 }
3040 #endif // XXXX
3041 
3042 /**
3043  * Constructor that initializes speciesMap and the FunctionDB object
3044  */
3046  mIgnoredSBMLMessages(),
3047  speciesMap(),
3048  functionDB(NULL),
3049  mIncompleteModel(false),
3050  mUnsupportedRuleFound(false),
3051  mUnsupportedRateRuleFound(false),
3052  mUnsupportedAssignmentRuleFound(false),
3053  mUnitOnNumberFound(false),
3054  mAssignmentToSpeciesReferenceFound(false),
3055  mLevel(0),
3056  mOriginalLevel(0),
3057  mVersion(0),
3058  sbmlIdMap(),
3059  mUsedFunctions(),
3060  mpDataModel(NULL),
3061  mpCopasiModel(NULL),
3062  mFunctionNameMapping(),
3063  mDivisionByCompartmentReactions(),
3064  mpImportHandler(NULL),
3065  mImportStep(0),
3066  mhImportStep(C_INVALID_INDEX),
3067  mTotalSteps(0),
3068  mSubstanceOnlySpecies(),
3069  mFastReactions(),
3070  mReactionsWithReplacedLocalParameters(),
3071  mExplicitelyTimeDependentFunctionDefinitions(),
3072  mIgnoredParameterUnits(),
3073  mStoichiometricExpressionMap(),
3074  mDelayFound(false),
3075  mPotentialAvogadroNumbers(),
3076  mAvogadroCreated(false),
3077  mImportCOPASIMIRIAM(true),
3078  mDelayNodeMap(),
3079  mUsedSBMLIds(),
3080  mUsedSBMLIdsPopulated(false),
3081  mAvogadroSet(false),
3082  mKnownCustomUserDefinedFunctions(),
3083 #if LIBSBML_VERSION >= 40100
3084  mpModelConversionFactor(NULL),
3085  mChemEqElementSpeciesIdMap(),
3086  mSpeciesConversionParameterMap(),
3087  mSBMLIdModelValueMap(),
3088  mSBMLSpeciesReferenceIds(),
3089  mRateRuleForSpeciesReferenceIgnored(false),
3090  mEventAssignmentForSpeciesReferenceIgnored(false),
3091  mConversionFactorFound(false),
3092 # if LIBSBML_VERSION >= 40200
3093  mEventPrioritiesIgnored(false),
3094  mInitialTriggerValues(false),
3095  mNonPersistentTriggerFound(false)
3096 # endif // LIBSBML_VERSION >= 40200
3097 #endif // LIBSBML_VERSION >= 40100
3098 {
3099  this->speciesMap = std::map<std::string, CMetab*>();
3100  this->functionDB = NULL;
3101  this->mIncompleteModel = false;
3102  this->mUnsupportedRuleFound = false;
3103  this->mUnsupportedRateRuleFound = false;
3104  this->mUnsupportedAssignmentRuleFound = false;
3105  this->mUnitOnNumberFound = false;
3106  this->mAssignmentToSpeciesReferenceFound = false;
3107  this->mpImportHandler = NULL;
3108  this->mDelayFound = false;
3109  this->mAvogadroCreated = false;
3110  // these data structures are used to handle the new conversion factores in
3111  // SBML L3 models and references to species references
3112  this->mpModelConversionFactor = NULL;
3113  this->mChemEqElementSpeciesIdMap.clear();
3114  this->mSpeciesConversionParameterMap.clear();
3115  this->mSBMLIdModelValueMap.clear();
3116  this->mRateRuleForSpeciesReferenceIgnored = false;
3117  this->mEventAssignmentForSpeciesReferenceIgnored = false;
3118  this->mConversionFactorFound = false;
3119  this->mEventPrioritiesIgnored = false;
3120  this->mInitialTriggerValues = false;
3121  this->mNonPersistentTriggerFound = false;
3122 
3123  this->mIgnoredSBMLMessages.insert(10501);
3124  this->mIgnoredSBMLMessages.insert(10512);
3125  this->mIgnoredSBMLMessages.insert(10513);
3126  this->mIgnoredSBMLMessages.insert(10522);
3127  this->mIgnoredSBMLMessages.insert(10533);
3128  this->mIgnoredSBMLMessages.insert(10541);
3129  this->mIgnoredSBMLMessages.insert(10551);
3130  this->mIgnoredSBMLMessages.insert(10562);
3131  this->mIgnoredSBMLMessages.insert(80701);
3132  this->mIgnoredSBMLMessages.insert(99505);
3133 }
3134 
3135 /**
3136  * Destructor that does nothing.
3137  */
3139 {}
3140 
3141 /**
3142  * This functions replaces all species nodes for species that are in the substanceOnlySpeciesVector.
3143  * With the node multiplied by the volume of the species compartment.
3144 void SBMLImporter::replaceSubstanceOnlySpeciesNodes(ConverterASTNode* node, const std::map<Species*, Compartment*>& substanceOnlySpecies)
3145 {
3146  if (node != NULL)
3147  {
3148  if (node->getType() == AST_NAME)
3149  {
3150  std::map<Species*, Compartment*>::const_iterator it = substanceOnlySpecies.begin();
3151  std::map<Species*, Compartment*>::const_iterator endIt = substanceOnlySpecies.end();
3152  while (it != endIt)
3153  {
3154  if (it->first->getId() == node->getName())
3155  {
3156  // replace node
3157  List* l = new List();
3158  ConverterASTNode* child1 = new ConverterASTNode(AST_NAME);
3159  child1->setName(node->getName());
3160  ConverterASTNode* child2 = new ConverterASTNode(AST_NAME);
3161  child2->setName(it->second->getId().c_str());
3162  l->add(child1);
3163  l->add(child2);
3164  node->setChildren(l);
3165  node->setType(AST_TIMES);
3166  break;
3167  }
3168  ++it;
3169  }
3170  }
3171  else
3172  {
3173  unsigned int counter;
3174  for (counter = 0;counter < node->getNumChildren();counter++)
3175  {
3176  this->replaceSubstanceOnlySpeciesNodes((ConverterASTNode*)node->getChild(counter), substanceOnlySpecies);
3177  }
3178  }
3179  }
3180 }
3181  */
3182 
3183 /**
3184  * Function reads an SBML file with libsbml and converts it to a Copasi CModel
3185  */
3186 CModel* SBMLImporter::readSBML(std::string filename,
3187  CFunctionDB* funDB,
3188  SBMLDocument*& pSBMLDocument,
3189  std::map<CCopasiObject*, SBase*>& copasi2sbmlmap,
3190  CListOfLayouts *& prLol,
3191  CCopasiDataModel* pDataModel)
3192 {
3193  // convert filename to the locale encoding
3194  std::ifstream file(CLocaleString::fromUtf8(filename).c_str());
3195 
3196  if (!file)
3197  {
3198  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 50, filename.c_str());
3199  }
3200 
3201  std::ostringstream stringStream;
3202  char c;
3203 
3204  while (file.get(c))
3205  {
3206  stringStream << c;
3207  }
3208 
3209  file.clear();
3210  file.close();
3211  return this->parseSBML(stringStream.str(), funDB,
3212  pSBMLDocument, copasi2sbmlmap, prLol, pDataModel);
3213 }
3214 
3215 /**
3216  * Function parses an SBML document with libsbml and converts it to a COPASI CModel
3217  * object which is returned. Deletion of the returned pointer is up to the
3218  * caller.
3219  */
3220 CModel*
3221 SBMLImporter::parseSBML(const std::string& sbmlDocumentText,
3222  CFunctionDB* funDB,
3223  SBMLDocument *& pSBMLDocument,
3224  std::map<CCopasiObject*, SBase*>& copasi2sbmlmap,
3225  CListOfLayouts *& prLol,
3226  CCopasiDataModel* pDataModel)
3227 {
3228  this->mUsedSBMLIdsPopulated = false;
3229  this->mAvogadroSet = false;
3230  mpDataModel = pDataModel;
3231  assert(mpDataModel != NULL);
3232 
3233  this->mpCopasiModel = NULL;
3234 
3235  if (funDB != NULL)
3236  {
3237  this->functionDB = funDB;
3238  SBMLReader* reader = new SBMLReader();
3239 
3240  mImportStep = 0;
3241 
3242  if (mpImportHandler)
3243  {
3244  mpImportHandler->setName("Importing SBML file...");
3245  mTotalSteps = 11;
3247  mImportStep,
3248  &mTotalSteps);
3249  }
3250 
3251  unsigned C_INT32 step = 0, totalSteps = 0;
3252  size_t hStep = C_INVALID_INDEX;
3253 
3254  if (this->mpImportHandler != 0)
3255  {
3256  step = 0;
3257  totalSteps = 1;
3258  hStep = mpImportHandler->addItem("Reading SBML file...",
3259  step,
3260  &totalSteps);
3261  }
3262 
3263  SBMLDocument* sbmlDoc = reader->readSBMLFromString(sbmlDocumentText);
3264 
3266 
3267  if (this->mpImportHandler != 0)
3268  {
3269  step = 0;
3270  totalSteps = 1;
3271  hStep = mpImportHandler->addItem("Checking consistency...",
3272  step,
3273  &totalSteps);
3274  }
3275 
3276  if (CCopasiRootContainer::getConfiguration()->validateUnits())
3277  sbmlDoc->setApplicableValidators(AllChecksON);
3278  else
3279  sbmlDoc->setApplicableValidators(AllChecksON & UnitsCheckOFF);
3280 
3281 #if LIBSBML_VERSION > 50800
3282  // libSBML is validating comp models after 5.8.0 this would throw an
3283  // error in case external references can't be resolved.
3284  sbmlDoc->setLocationURI(pDataModel->getReferenceDirectory());
3285 #endif
3286 
3287  bool checkResult = sbmlDoc->checkConsistency();
3288 
3289 #if LIBSBML_VERSION > 50800
3290 
3291  sbmlDoc->setLocationURI(pDataModel->getReferenceDirectory());
3292 
3293  // the new libsbml includes complete layout validation, and flags many
3294  // things as errors, that would cause COPASI to reject the model (even
3295  // though the layout code has been written to anticipate all these possible
3296  // error sources). Thus downgrade these error messages
3297  sbmlDoc->getErrorLog()->changeErrorSeverity(LIBSBML_SEV_ERROR, LIBSBML_SEV_WARNING, "layout");
3298 #endif
3299 
3301 
3302  if (checkResult != 0)
3303  {
3304  int fatal = -1;
3305  unsigned int i, iMax = sbmlDoc->getNumErrors();
3306 
3307  for (i = 0; (i < iMax) && (fatal == -1); ++i)
3308  {
3309  const XMLError* pSBMLError = sbmlDoc->getError(i);
3310 
3311  // we check if the model contained a required package
3312  if (sbmlDoc->getLevel() > 2 && pSBMLError->getErrorId() == 99107)
3313  {
3315  }
3316 
3318 
3319  switch (pSBMLError->getSeverity())
3320  {
3321  case LIBSBML_SEV_INFO:
3322 
3323  if (mIgnoredSBMLMessages.find(pSBMLError->getErrorId()) != mIgnoredSBMLMessages.end())
3324  {
3325  messageType = CCopasiMessage::WARNING_FILTERED;
3326  }
3327  else
3328  {
3329  messageType = CCopasiMessage::WARNING;
3330  }
3331 
3332  CCopasiMessage(messageType, MCSBML + 40, "INFO", pSBMLError->getErrorId(), pSBMLError->getLine(), pSBMLError->getColumn(), pSBMLError->getMessage().c_str());
3333  break;
3334 
3335  case LIBSBML_SEV_WARNING:
3336 
3337 #if LIBSBML_VERSION > 50800
3338 
3339  // filter layout warnings always
3340  if (pSBMLError->getPackage() == "layout")
3341  messageType = CCopasiMessage::WARNING_FILTERED;
3342  else
3343 #endif
3344  if (mIgnoredSBMLMessages.find(pSBMLError->getErrorId()) != mIgnoredSBMLMessages.end())
3345  {
3346  messageType = CCopasiMessage::WARNING_FILTERED;
3347  }
3348  else
3349  {
3350  messageType = CCopasiMessage::WARNING;
3351  }
3352 
3353  CCopasiMessage(messageType, MCSBML + 40, "WARNING", pSBMLError->getErrorId(), pSBMLError->getLine(), pSBMLError->getColumn(), pSBMLError->getMessage().c_str());
3354  break;
3355 
3356  case LIBSBML_SEV_ERROR:
3357 
3358  if (mIgnoredSBMLMessages.find(pSBMLError->getErrorId()) != mIgnoredSBMLMessages.end())
3359  {
3360  messageType = CCopasiMessage::ERROR_FILTERED;
3361  }
3362 
3363  CCopasiMessage(messageType, MCSBML + 40, "ERROR", pSBMLError->getErrorId(), pSBMLError->getLine(), pSBMLError->getColumn(), pSBMLError->getMessage().c_str());
3364  break;
3365 
3366  case LIBSBML_SEV_FATAL:
3367 
3368  // treat unknown as fatal
3369  default:
3370 
3371  //CCopasiMessage(CCopasiMessage::TRACE, MCSBML + 40,"FATAL",pSBMLError->getLine(),pSBMLError->getColumn(),pSBMLError->getMessage().c_str());
3372  if (pSBMLError->getErrorId() == 10804)
3373  {
3374  // this error indicates a problem with a notes element
3375  // although libsbml flags this as fatal, we would still
3376  // like to read the model
3377  CCopasiMessage(messageType, MCSBML + 40, "ERROR", pSBMLError->getErrorId(), pSBMLError->getLine(), pSBMLError->getColumn(), pSBMLError->getMessage().c_str());
3378  }
3379  else
3380  {
3381  fatal = i;
3382  }
3383 
3384  break;
3385  }
3386 
3387  //std::cerr << pSBMLError->getMessage() << std::endl;
3388  }
3389 
3390  if (fatal != -1)
3391  {
3392  const XMLError* pSBMLError = sbmlDoc->getError(fatal);
3394  pSBMLError->getLine(),
3395  pSBMLError->getColumn(),
3396  pSBMLError->getMessage().c_str());
3397 
3399 
3400  return NULL;
3401  }
3402  }
3403  else
3404  {
3405  if (sbmlDoc->getLevel() > 2)
3406  {
3407  // we check if the model contained a required package
3408  unsigned int i, iMax = sbmlDoc->getNumErrors();
3409 
3410  for (i = 0; i < iMax ; ++i)
3411  {
3412  const XMLError* pSBMLError = sbmlDoc->getError(i);
3413 
3414  if (pSBMLError->getErrorId() == 99107)
3415  {
3416  if (pSBMLError->getMessage().find("'spatial'") != std::string::npos)
3417  {
3419  "The model you tried to open requires the SBML spatial package. "
3420  "This version of COPASI does not support spatial models.");
3421  }
3422  else
3423  {
3425  }
3426  }
3427  }
3428  }
3429  }
3430 
3431 #if LIBSBML_VERSION >= 50700
3432 
3433  if (sbmlDoc->getPlugin("comp") != NULL && sbmlDoc->isSetPackageRequired("comp"))
3434  {
3435 
3436  if (sbmlDoc->getNumErrors(LIBSBML_SEV_ERROR) > 0)
3437  {
3438  std::string message =
3439  "The SBML model you are trying to import uses the Hierarchical Modeling extension. "
3440  "In order to import this model in COPASI, it has to be flattened, however flattening is "
3441  "not possible because of the following errors:\n" + sbmlDoc->getErrorLog()->toString();
3442 
3443  // escape potential percent signs that would lead to an error when passed to CCopasiMessage
3444  CPrefixNameTransformer::replaceStringInPlace(message, "%", "%%");
3445 
3446  CCopasiMessage(CCopasiMessage::EXCEPTION, message.c_str());
3447  }
3448 
3449 #if LIBSBML_VERSION >= 50903 && LIBSBML_HAS_PACKAGE_COMP
3450  // apply the name transformer
3451  CompModelPlugin* mPlug = dynamic_cast<CompModelPlugin*>(sbmlDoc->getModel()->getPlugin("comp"));
3452  CPrefixNameTransformer trans;
3453 
3454  if (mPlug != NULL)
3455  {
3456  mPlug->setTransformer(&trans);
3457  }
3458 
3459 #endif //LIBSBML_VERSION >= 50903 && LIBSBML_HAS_PACKAGE_COMP
3460 
3461  // the sbml comp package is used, and the required flag is set, so it stands to reason
3462  // that we need to flatten the document
3463  sbmlDoc->getErrorLog()->clearLog();
3464 
3465  ConversionProperties props;
3466  props.addOption("flatten comp");
3467  props.addOption("leavePorts", false);
3468  props.addOption("basePath", pDataModel->getReferenceDirectory());
3469 
3470  if (sbmlDoc->convert(props) != LIBSBML_OPERATION_SUCCESS)
3471  {
3472  std::string message =
3473  "The SBML model you are trying to import uses the Hierarchical Modeling extension. "
3474  "In order to import this model in COPASI, it has to be flattened, however flattening failed";
3475 
3476  if (sbmlDoc->getNumErrors() == 0)
3477  message += ".";
3478  else
3479  message += " with the following errors:\n" + sbmlDoc->getErrorLog()->toString();
3480 
3481  CCopasiMessage(CCopasiMessage::EXCEPTION, message.c_str());
3482  }
3483  }
3484 
3485 #endif
3486 
3487  if (sbmlDoc->getModel() == NULL)
3488  {
3490 
3492 
3493  return NULL;
3494  }
3495 
3496  delete reader;
3497  pSBMLDocument = sbmlDoc;
3498  this->mLevel = pSBMLDocument->getLevel();
3499  // remember the original level of the document because we convert Level 1 documents to Level 2
3500  // For the import of rules, we need to remember that is was actually a level 1 document
3501  // because otherwise we throw error messages on rules on parameters since the parameters
3502  // have been set to constant by the conversation to Level 2
3503  this->mOriginalLevel = this->mLevel;
3504  this->mVersion = pSBMLDocument->getVersion();
3505 
3506  if (mLevel == 1)
3507  {
3508  unsigned int i, iMax = pSBMLDocument->getModel()->getNumCompartments();
3509 
3510  for (i = 0; i < iMax; ++i)
3511  {
3512  Compartment* pCompartment = pSBMLDocument->getModel()->getCompartment(i);
3513  pCompartment->setSize(pCompartment->getVolume());
3514  }
3515 
3516  pSBMLDocument->setLevelAndVersion(2, 1);
3517  mLevel = pSBMLDocument->getLevel();
3518  }
3519 
3520  this->mpCopasiModel = this->createCModelFromSBMLDocument(sbmlDoc, copasi2sbmlmap);
3521 
3522  prLol = new CListOfLayouts("ListOfLayouts", mpDataModel);
3523  Model* sbmlmodel = pSBMLDocument->getModel();
3524 
3525  if (sbmlmodel && prLol)
3526  {
3527  LayoutModelPlugin *lmPlugin = (LayoutModelPlugin*)sbmlmodel->getPlugin("layout");
3528 
3529  if (lmPlugin != NULL)
3531  *lmPlugin->getListOfLayouts(),
3532  copasi2sbmlmap);
3533  }
3534  }
3535  else
3536  {
3538 
3539  fatalError();
3540  }
3541 
3543 
3544  return this->mpCopasiModel;
3545 }
3546 
3547 /**
3548  * Returns the copasi QuantityUnit corresponding to the given SBML
3549  * Substance UnitDefinition.
3550  */
3551 std::pair<CModel::QuantityUnit, bool>
3552 SBMLImporter::handleSubstanceUnit(const UnitDefinition* uDef)
3553 {
3554  bool result = false;
3556 
3557  if (uDef == NULL)
3558  {
3559  //DebugFile << "Argument to handleSubstanceUnit is NULL pointer." << std::endl;
3560  fatalError();
3561  }
3562 
3563  if (uDef->getNumUnits() == 1)
3564  {
3565  const Unit* u = uDef->getUnit(0);
3566 
3567  if (u == NULL)
3568  {
3569  //DebugFile << "Expected Unit, got NULL pointer." << std::endl;
3570  fatalError();
3571  }
3572 
3573  if ((u->getKind() == UNIT_KIND_MOLE)
3574  || u->getKind() == UNIT_KIND_AVOGADRO
3575  )
3576  {
3577  double multiplier = u->getMultiplier();
3578  int scale = u->getScale();
3579 
3580  if (multiplier != 1)
3581  {
3582  // check if the multiplier is a multiple of 10
3583  // so that we might be able to convert it to a scale that makes
3584  // sense
3585  double tmp = log10(multiplier);
3586 
3587  if (areApproximatelyEqual(tmp, round(tmp)))
3588  {
3589  scale += (int) round(tmp);
3590  multiplier = 1;
3591  }
3592  }
3593 
3594  if ((u->getExponent() == 1) &&
3595  areApproximatelyEqual(multiplier, 1.0) &&
3596  ((scale % 3) == 0) &&
3597  (scale < 1) &&
3598  (scale > -16))
3599  {
3600  switch (scale)
3601  {
3602  case 0:
3603  qUnit = CModel::Mol;
3604  result = true;
3605  break;
3606 
3607  case - 3:
3608  qUnit = CModel::mMol;
3609  result = true;
3610  break;
3611 
3612  case - 6:
3613  qUnit = CModel::microMol;
3614  result = true;
3615  break;
3616 
3617  case - 9:
3618  qUnit = CModel::nMol;
3619  result = true;
3620  break;
3621 
3622  case - 12:
3623  qUnit = CModel::pMol;
3624  result = true;
3625  break;
3626 
3627  case - 15:
3628  qUnit = CModel::fMol;
3629  result = true;
3630  break;
3631 
3632  default:
3633  //DebugFile << "Error. This value should never have been reached for the scale of the liter unit." << std::endl;
3634  result = false;
3635  break;
3636  }
3637  }
3638  else
3639  {
3640  result = false;
3641  }
3642  }
3643  else if ((u->getKind() == UNIT_KIND_ITEM))
3644  {
3645  double multiplier = u->getMultiplier();
3646  int scale = u->getScale();
3647 
3648  if (multiplier != 1)
3649  {
3650  // check if the multiplier is a multiple of 10
3651  // so that we might be able to convert it to a scale that makes
3652  // sense
3653  double tmp = log10(multiplier);
3654 
3655  if (areApproximatelyEqual(tmp, round(tmp)))
3656  {
3657  scale += (int)round(tmp);
3658  multiplier = 1;
3659  }
3660  }
3661 
3662  if ((u->getExponent() == 1) &&
3663  areApproximatelyEqual(multiplier, 1.0) &&
3664  (scale == 0 || scale == 1))
3665  {
3666  if (u->getScale() == 1)
3667  {
3669  }
3670  else
3671  {
3672  result = true;
3673  qUnit = CModel::number;
3674  }
3675  }
3676  else
3677  {
3678  result = false;
3679  }
3680  }
3681  else if ((u->getKind() == UNIT_KIND_DIMENSIONLESS))
3682  {
3683  double multiplier = u->getMultiplier();
3684  int scale = u->getScale();
3685 
3686  if (multiplier != 1)
3687  {
3688  // check if the multiplier is a multiple of 10
3689  // so that we might be able to convert it to a scale that makes
3690  // sense
3691  double tmp = log10(multiplier);
3692 
3693  if (areApproximatelyEqual(tmp, round(tmp)))
3694  {
3695  scale += (int)round(tmp);
3696  multiplier = 1;
3697  }
3698  }
3699 
3700  if ((u->getExponent() == 1) &&
3701  areApproximatelyEqual(multiplier, 1.0) &&
3702  scale == 0)
3703  {
3704  result = true;
3706  }
3707  else
3708  {
3709  result = false;
3710  }
3711  }
3712  else
3713  {
3714  result = false;
3715  }
3716  }
3717  else
3718  {
3719  result = false;
3720  }
3721 
3722  return std::make_pair(qUnit, result);
3723 }
3724 
3725 /**
3726  * Returns the copasi TimeUnit corresponding to the given SBML Time
3727  * UnitDefinition.
3728  */
3729 std::pair<CModel::TimeUnit, bool>
3730 SBMLImporter::handleTimeUnit(const UnitDefinition* uDef)
3731 {
3732  bool result = false;
3733  CModel::TimeUnit tUnit = CModel::s;
3734 
3735  if (uDef == NULL)
3736  {
3737  //DebugFile << "Argument to handleTimeUnit is NULL pointer." << std::endl;
3738  fatalError();
3739  }
3740 
3741  if (uDef->getNumUnits() == 1)
3742  {
3743  const Unit* u = uDef->getUnit(0);
3744 
3745  if (u == NULL)
3746  {
3747  //DebugFile << "Expected Unit, got NULL pointer." << std::endl;
3748  fatalError();
3749  }
3750 
3751  if ((u->getKind() == UNIT_KIND_SECOND))
3752  {
3753  double multiplier = u->getMultiplier();
3754  int scale = u->getScale();
3755 
3756  // this is more difficult, we have to check several possible
3757  // combinations of scales and multipliers.
3758  // Valid multipliers are 60, 3600 and 86400 which correspond to a
3759  // minute an hour and a day each of those has to have a scale of 0
3760  if (scale == 0)
3761  {
3762  // check if the multiplier is 1, 60, 3600 or 86400
3763  // if not, try to make the multiplier 1
3764  if (!areApproximatelyEqual(multiplier, 1.0) &&
3765  !areApproximatelyEqual(multiplier, 60.0) &&
3766  !areApproximatelyEqual(multiplier, 3600.0) &&
3767  !areApproximatelyEqual(multiplier, 86400.0))
3768  {
3769  double tmp = log10(multiplier);
3770 
3771  if (!areApproximatelyEqual(tmp, round(tmp)))
3772  {
3773  result = false;
3774  }
3775  else
3776  {
3777  multiplier = 1;
3778  scale += (int)round(tmp);
3779  }
3780  }
3781  }
3782  else
3783  {
3784  // the multiplier must be 1
3785  if (!areApproximatelyEqual(multiplier, 1.0))
3786  {
3787  // make the multiplier 1 and check if the scale is an integer,
3788  // if not, try to make the 0 and check if the multiplier becomes one
3789  // of the valid multipliers 1,60, 3600 or 86400
3790  double tmp = log10(multiplier);
3791 
3792  if (!areApproximatelyEqual(tmp, round(tmp)))
3793  {
3794  // try to make the scale 0
3795  multiplier *= pow(10.0, (double)scale);
3796  scale = 0;
3797  }
3798  else
3799  {
3800  // make the multiplier 1
3801  multiplier = 1;
3802  scale += (int) round(tmp);
3803  }
3804  }
3805  }
3806 
3807  if ((u->getExponent() == 1) && ((scale % 3) == 0) && (scale < 1) && (scale > -16))
3808  {
3809  if (areApproximatelyEqual(multiplier, 1.0))
3810  {
3811  switch (scale)
3812  {
3813  case 0:
3814  tUnit = CModel::s;
3815  result = true;
3816  break;
3817 
3818  case - 3:
3819  tUnit = CModel::ms;
3820  result = true;
3821  break;
3822 
3823  case - 6:
3824  tUnit = CModel::micros;
3825  result = true;
3826  break;
3827 
3828  case - 9:
3829  tUnit = CModel::ns;
3830  result = true;
3831  break;
3832 
3833  case - 12:
3834  tUnit = CModel::ps;
3835  result = true;
3836  break;
3837 
3838  case - 15:
3839  tUnit = CModel::fs;
3840  result = true;
3841  break;
3842 
3843  default:
3844  //DebugFile << "Error. This value should never have been reached for the scale of the time unit." << std::endl;
3845  result = false;
3846  break;
3847  }
3848  }
3849  else if ((scale == 0) &&
3850  areApproximatelyEqual(multiplier, 60.0))
3851  {
3852  tUnit = CModel::min;
3853  result = true;
3854  }
3855  else if ((scale == 0) &&
3856  areApproximatelyEqual(multiplier, 3600.0))
3857  {
3858  tUnit = CModel::h;
3859  result = true;
3860  }
3861  else if ((scale == 0) &&
3862  areApproximatelyEqual(multiplier, 86400.0))
3863  {
3864  tUnit = CModel::d;
3865  result = true;
3866  }
3867  else
3868  {
3869  result = false;
3870  }
3871  }
3872  else
3873  {
3874  result = false;
3875  }
3876  }
3877  else if ((u->getKind() == UNIT_KIND_DIMENSIONLESS))
3878  {
3879  double multiplier = u->getMultiplier();
3880  int scale = u->getScale();
3881 
3882  if (multiplier != 1)
3883  {
3884  // check if the multiplier is a multiple of 10
3885  // so that we might be able to convert it to a scale that makes
3886  // sense
3887  double tmp = log10(multiplier);
3888 
3889  if (areApproximatelyEqual(tmp, round(tmp)))
3890  {
3891  scale += (int)round(tmp);
3892  multiplier = 1;
3893  }
3894  }
3895 
3896  if ((u->getExponent() == 1) &&
3897  areApproximatelyEqual(multiplier, 1.0) &&
3898  scale == 0)
3899  {
3900  result = true;
3901  tUnit = CModel::dimensionlessTime;
3902  }
3903  else
3904  {
3905  result = false;
3906  }
3907  }
3908  else
3909  {
3910  result = false;
3911  }
3912  }
3913  else
3914  {
3915  result = false;
3916  }
3917 
3918  return std::make_pair(tUnit, result);
3919 }
3920 
3921 /**
3922  * Returns the copasi LengthUnit corresponding to the given SBML length
3923  * UnitDefinition.
3924  */
3925 std::pair<CModel::LengthUnit, bool>
3926 SBMLImporter::handleLengthUnit(const UnitDefinition* uDef)
3927 {
3928  bool result = false;
3929  CModel::LengthUnit lUnit = CModel::m;
3930 
3931  if (uDef == NULL)
3932  {
3933  //DebugFile << "Argument to handleLengthUnit is NULL pointer." << std::endl;
3934  fatalError();
3935  }
3936 
3937  if (uDef->getNumUnits() == 1)
3938  {
3939  const Unit* u = uDef->getUnit(0);
3940 
3941  if (u == NULL)
3942  {
3943  //DebugFile << "Expected Unit, got NULL pointer." << std::endl;
3944  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 54, "length", uDef->getId().c_str());
3945  }
3946 
3947  if ((u->getKind() == UNIT_KIND_METRE || u->getKind() == UNIT_KIND_METER) && u->getExponent() == 1)
3948  {
3949  double multiplier = u->getMultiplier();
3950  int scale = u->getScale();
3951 
3952  if (multiplier != 1)
3953  {
3954  // check if the multiplier is a multiple of 10
3955  // so that we might be able to convert it to a scale that makes
3956  // sense
3957  double tmp = log10(multiplier);
3958 
3959  if (areApproximatelyEqual(tmp, round(tmp)))
3960  {
3961  scale += (int) round(tmp);
3962  multiplier = 1;
3963  }
3964  }
3965 
3966  if (areApproximatelyEqual(multiplier, 1.0) &&
3967  ((scale % 3) == 0 || scale == -1 || scale == -2) &&
3968  (scale < 1) &&
3969  (scale > -16))
3970  {
3971  switch (scale)
3972  {
3973  case 0:
3974  lUnit = CModel::m;
3975  result = true;
3976  break;
3977 
3978  case - 1:
3979  lUnit = CModel::dm;
3980  result = true;
3981  break;
3982 
3983  case - 2:
3984  lUnit = CModel::cm;
3985  result = true;
3986  break;
3987 
3988  case - 3:
3989  lUnit = CModel::mm;
3990  result = true;
3991  break;
3992 
3993  case - 6:
3994  lUnit = CModel::microm;
3995  result = true;
3996  break;
3997 
3998  case - 9:
3999  lUnit = CModel::nm;
4000  result = true;
4001  break;
4002 
4003  case - 12:
4004  lUnit = CModel::pm;
4005  result = true;
4006  break;
4007 
4008  case - 15:
4009  lUnit = CModel::fm;
4010  result = true;
4011  break;
4012 
4013  default:
4014  //DebugFile << "Error. This value should never have been reached for the scale of the liter unit." << std::endl;
4015  result = false;
4016  break;
4017  }
4018  }
4019  else
4020  {
4021  result = false;
4022  }
4023  }
4024  else if ((u->getKind() == UNIT_KIND_DIMENSIONLESS))
4025  {
4026  double multiplier = u->getMultiplier();
4027  int scale = u->getScale();
4028 
4029  if (multiplier != 1)
4030  {
4031  // check if the multiplier is a multiple of 10
4032  // so that we might be able to convert it to a scale that makes
4033  // sense
4034  double tmp = log10(multiplier);
4035 
4036  if (areApproximatelyEqual(tmp, round(tmp)))
4037  {
4038  scale += (int)round(tmp);
4039  multiplier = 1;
4040  }
4041  }
4042 
4043  if ((u->getExponent() == 1) &&
4044  areApproximatelyEqual(multiplier, 1.0) &&
4045  scale == 0)
4046  {
4047  result = true;
4049  }
4050  else
4051  {
4052  result = false;
4053  }
4054  }
4055  else
4056  {
4057  result = false;
4058  }
4059  }
4060  else
4061  {
4062  result = false;
4063  }
4064 
4065  return std::make_pair(lUnit, result);
4066 }
4067 
4068 /**
4069  * Returns the copasi AreaUnit corresponding to the given SBML area
4070  * UnitDefinition.
4071  */
4072 std::pair<CModel::AreaUnit, bool>
4073 SBMLImporter::handleAreaUnit(const UnitDefinition* uDef)
4074 {
4075  bool result = false;
4076  CModel::AreaUnit aUnit = CModel::m2;
4077 
4078  if (uDef == NULL)
4079  {
4080  //DebugFile << "Argument to handleLengthUnit is NULL pointer." << std::endl;
4081  fatalError();
4082  }
4083 
4084  if (uDef->getNumUnits() == 1)
4085  {
4086  const Unit* u = uDef->getUnit(0);
4087 
4088  if (u == NULL)
4089  {
4090  //DebugFile << "Expected Unit, got NULL pointer." << std::endl;
4091  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 54, "area", uDef->getId().c_str());
4092  }
4093 
4094  if ((u->getKind() == UNIT_KIND_METRE || u->getKind() == UNIT_KIND_METER) && u->getExponent() == 2)
4095  {
4096  double multiplier = u->getMultiplier();
4097  int scale = u->getScale();
4098 
4099  if (multiplier != 1)
4100  {
4101  // check if the multiplier is a multiple of 10
4102  // so that we might be able to convert it to a scale that makes
4103  // sense
4104  double tmp = log10(multiplier);
4105 
4106  if (areApproximatelyEqual(tmp, round(tmp)))
4107  {
4108  scale += (int) round(tmp);
4109  multiplier = 1;
4110  }
4111  }
4112 
4113  if (areApproximatelyEqual(multiplier, 1.0) &&
4114  ((scale % 3) == 0 || scale == -1 || scale == -2) &&
4115  (scale < 1) &&
4116  (scale > -16))
4117  {
4118  switch (scale)
4119  {
4120  case 0:
4121  aUnit = CModel::m2;
4122  result = true;
4123  break;
4124 
4125  case - 1:
4126  aUnit = CModel::dm2;
4127  result = true;
4128  break;
4129 
4130  case - 2:
4131  aUnit = CModel::cm2;
4132  result = true;
4133  break;
4134 
4135  case - 3:
4136  aUnit = CModel::mm2;
4137  result = true;
4138  break;
4139 
4140  case - 6:
4141  aUnit = CModel::microm2;
4142  result = true;
4143  break;
4144 
4145  case - 9:
4146  aUnit = CModel::nm2;
4147  result = true;
4148  break;
4149 
4150  case - 12:
4151  aUnit = CModel::pm2;
4152  result = true;
4153  break;
4154 
4155  case - 15:
4156  aUnit = CModel::fm2;
4157  result = true;
4158  break;
4159 
4160  default:
4161  //DebugFile << "Error. This value should never have been reached for the scale of the liter unit." << std::endl;
4162  result = false;
4163  break;
4164  }
4165  }
4166  else
4167  {
4168  result = false;
4169  }
4170  }
4171  else if ((u->getKind() == UNIT_KIND_DIMENSIONLESS))
4172  {
4173  double multiplier = u->getMultiplier();
4174  int scale = u->getScale();
4175 
4176  if (multiplier != 1)
4177  {
4178  // check if the multiplier is a multiple of 10
4179  // so that we might be able to convert it to a scale that makes
4180  // sense
4181  double tmp = log10(multiplier);
4182 
4183  if (areApproximatelyEqual(tmp, round(tmp)))
4184  {
4185  scale += (int)round(tmp);
4186  multiplier = 1;
4187  }
4188  }
4189 
4190  if ((u->getExponent() == 1) &&
4191  areApproximatelyEqual(multiplier, 1.0) &&
4192  scale == 0)
4193  {
4194  result = true;
4195  aUnit = CModel::dimensionlessArea;
4196  }
4197  else
4198  {
4199  result = false;
4200  }
4201  }
4202  else
4203  {
4204  result = false;
4205  }
4206  }
4207  else
4208  {
4209  result = false;
4210  }
4211 
4212  return std::make_pair(aUnit, result);
4213 }
4214 
4215 /**
4216  * Returns the copasi VolumeUnit corresponding to the given SBML Volume
4217  * UnitDefinition.
4218  */
4219 std::pair<CModel::VolumeUnit, bool>
4220 SBMLImporter::handleVolumeUnit(const UnitDefinition* uDef)
4221 {
4222  // simplify the Unitdefiniton first if this normalizes
4223  // the scale and the multiplier
4224  bool result = false;
4225  CModel::VolumeUnit vUnit = CModel::l;
4226 
4227  if (uDef == NULL)
4228  {
4229  //DebugFile << "Argument to handleVolumeUnit is NULL pointer." << std::endl;
4230  fatalError();
4231  }
4232 
4233  if (uDef->getNumUnits() == 1)
4234  {
4235  const Unit* u = uDef->getUnit(0);
4236 
4237  if (u == NULL)
4238  {
4239  //DebugFile << "Expected Unit, got NULL pointer." << std::endl;
4240  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 54, "volume", uDef->getId().c_str());
4241  }
4242 
4243  if (((u->getKind() == UNIT_KIND_LITER) || (u->getKind() == UNIT_KIND_LITRE)) && u->getExponent() == 1)
4244  {
4245  double multiplier = u->getMultiplier();
4246  int scale = u->getScale();
4247 
4248  if (multiplier != 1)
4249  {
4250  // check if the multiplier is a multiple of 10
4251  // so that we might be able to convert it to a scale that makes
4252  // sense
4253  double tmp = log10(multiplier);
4254 
4255  if (areApproximatelyEqual(tmp, round(tmp)))
4256  {
4257  scale += (int)round(tmp);
4258  multiplier = 1;
4259  }
4260  }
4261 
4262  if (areApproximatelyEqual(multiplier, 1.0) &&
4263  ((scale % 3) == 0) &&
4264  (scale < 1) &&
4265  (scale > -16))
4266  {
4267  switch (scale)
4268  {
4269  case 0:
4270  vUnit = CModel::l;
4271  result = true;
4272  break;
4273 
4274  case - 3:
4275  vUnit = CModel::ml;
4276  result = true;
4277  break;
4278 
4279  case - 6:
4280  vUnit = CModel::microl;
4281  result = true;
4282  break;
4283 
4284  case - 9:
4285  vUnit = CModel::nl;
4286  result = true;
4287  break;
4288 
4289  case - 12:
4290  vUnit = CModel::pl;
4291  result = true;
4292  break;
4293 
4294  case - 15:
4295  vUnit = CModel::fl;
4296  result = true;
4297  break;
4298 
4299  default:
4300  //DebugFile << "Error. This value should never have been reached for the scale of the liter unit." << std::endl;
4301  result = false;
4302  break;
4303  }
4304  }
4305  else
4306  {
4307  result = false;
4308  }
4309  }
4310  else if (((u->getKind() == UNIT_KIND_METER) || (u->getKind() == UNIT_KIND_METRE)) && u->getExponent() == 3)
4311  {
4312  double multiplier = u->getMultiplier();
4313  int scale = u->getScale();
4314 
4315  if (multiplier != 1)
4316  {
4317  // check if the multiplier is a multiple of 10
4318  // so that we might be able to convert it to a scale that makes
4319  // sense
4320  double tmp = log10(multiplier);
4321 
4322  if (areApproximatelyEqual(tmp, round(tmp)))
4323  {
4324  scale += (int)round(tmp);
4325  multiplier = 1;
4326  }
4327  }
4328 
4329  if (areApproximatelyEqual(multiplier, 1.0) &&
4330  (scale == 0))
4331  {
4332  vUnit = CModel::m3;
4333  result = true;
4334  }
4335  else
4336  {
4337  // try to convert to liter
4338  Unit* pLitreUnit = convertSBMLCubicmetresToLitres(u);
4339 
4340  if (pLitreUnit != NULL &&
4341  pLitreUnit->getExponent() == 1 &&
4342  (pLitreUnit->getScale() % 3 == 0) &&
4343  (pLitreUnit->getScale() < 1) &&
4344  (pLitreUnit->getScale() > -16) &&
4345  areApproximatelyEqual(pLitreUnit->getMultiplier(), 1.0))
4346  {
4347  switch (pLitreUnit->getScale())
4348  {
4349  case 0:
4350  vUnit = CModel::l;
4351  result = true;
4352  break;
4353 
4354  case - 3:
4355  vUnit = CModel::ml;
4356  result = true;
4357  break;
4358 
4359  case - 6:
4360  vUnit = CModel::microl;
4361  result = true;
4362  break;
4363 
4364  case - 9:
4365  vUnit = CModel::nl;
4366  result = true;
4367  break;
4368 
4369  case - 12:
4370  vUnit = CModel::pl;
4371  result = true;
4372  break;
4373 
4374  case - 15:
4375  vUnit = CModel::fl;
4376  result = true;
4377  break;
4378 
4379  default:
4380  result = false;
4381  break;
4382  }
4383  }
4384  else
4385  {
4386  result = false;
4387  }
4388 
4389  delete pLitreUnit;
4390  }
4391  }
4392  else if ((u->getKind() == UNIT_KIND_DIMENSIONLESS))
4393  {
4394  double multiplier = u->getMultiplier();
4395  int scale = u->getScale();
4396 
4397  if (multiplier != 1)
4398  {
4399  // check if the multiplier is a multiple of 10
4400  // so that we might be able to convert it to a scale that makes
4401  // sense
4402  double tmp = log10(multiplier);
4403 
4404  if (areApproximatelyEqual(tmp, round(tmp)))
4405  {
4406  scale += (int)round(tmp);
4407  multiplier = 1;
4408  }
4409  }
4410 
4411  if ((u->getExponent() == 1) &&
4412  areApproximatelyEqual(multiplier, 1.0) &&
4413  scale == 0)
4414  {
4415  result = true;
4417  }
4418  else
4419  {
4420  result = false;
4421  }
4422  }
4423  else
4424  {
4425  result = false;
4426  }
4427  }
4428  else
4429  {
4430  result = false;
4431  }
4432 
4433  return std::make_pair(vUnit, result);
4434 }
4435 
4436 CModelValue* SBMLImporter::createCModelValueFromParameter(const Parameter* sbmlParameter, CModel* copasiModel, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
4437 {
4438  if (sbmlParameter->isSetUnits())
4439  {
4440  /* !!!! create a warning that the units will be ignored. */
4441  //CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 26, sbmlParameter->getId().c_str());
4442  mIgnoredParameterUnits.push_back(sbmlParameter->getId());
4443  const_cast<Parameter*>(sbmlParameter)->unsetUnits();
4444  }
4445 
4446  std::string name = sbmlParameter->getName();
4447 
4448  if (!sbmlParameter->isSetName())
4449  {
4450  name = sbmlParameter->getId();
4451  }
4452 
4453  std::string appendix = "";
4454  unsigned int counter = 2;
4455  std::ostringstream numberStream;
4456 
4457  while (copasiModel->getModelValues().getIndex(name + appendix) != C_INVALID_INDEX)
4458  {
4459  numberStream.str("");
4460  numberStream << "_" << counter;
4461  counter++;
4462  appendix = numberStream.str();
4463  }
4464 
4465  std::string sbmlId;
4466 
4467  if (this->mLevel == 1)
4468  {
4469  sbmlId = sbmlParameter->getName();
4470  }
4471  else
4472  {
4473  sbmlId = sbmlParameter->getId();
4474  }
4475 
4476  CModelValue* pMV = copasiModel->createModelValue(name + appendix, 0.0);
4477  copasi2sbmlmap[pMV] = const_cast<Parameter*>(sbmlParameter);
4478  pMV->setSBMLId(sbmlId);
4479  SBMLImporter::importMIRIAM(sbmlParameter, pMV);
4480  SBMLImporter::importNotes(pMV, sbmlParameter);
4481 
4482  if (this->mLevel > 2)
4483  {
4484  this->mSBMLIdModelValueMap[sbmlId] = pMV;
4485  }
4486 
4487  return pMV;
4488 }
4489 
4490 bool SBMLImporter::sbmlId2CopasiCN(ASTNode* pNode, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap, CCopasiParameterGroup& pParamGroup)
4491 {
4492  // TODO CRITICAL We need to use a node iterator
4493 
4494  bool success = true;
4495  unsigned int i, iMax = pNode->getNumChildren();
4496 
4497  if (pNode->getType() == AST_NAME)
4498  {
4499  Compartment* pSBMLCompartment = NULL;
4500  Species* pSBMLSpecies = NULL;
4501  Reaction* pSBMLReaction = NULL;
4502  Parameter* pSBMLParameter = NULL;
4503  std::string sbmlId;
4504  std::string name = pNode->getName();
4505  CCopasiParameter* pParam = pParamGroup.getParameter(name);
4506 
4507  std::map<std::string, double>::const_iterator speciesReference = mSBMLSpeciesReferenceIds.find(name);
4508 
4509  if (speciesReference != mSBMLSpeciesReferenceIds.end())
4510  {
4511  // replace the name with the value
4512  pNode->setType(AST_REAL);
4513  pNode->setValue(mSBMLSpeciesReferenceIds[name]);
4514  }
4515  else if (pParam)
4516  {
4517  pNode->setName(pParam->getCN().c_str());
4518  }
4519  else
4520  {
4521  std::map<CCopasiObject*, SBase*>::iterator it = copasi2sbmlmap.begin();
4522  std::map<CCopasiObject*, SBase*>::iterator endIt = copasi2sbmlmap.end();
4523  bool found = false;
4524 
4525  while (it != endIt)
4526  {
4527  int type = it->second->getTypeCode();
4528 
4529  switch (type)
4530  {
4531  case SBML_COMPARTMENT:
4532  pSBMLCompartment = dynamic_cast<Compartment*>(it->second);
4533 
4534  if (this->mLevel == 1)
4535  {
4536  sbmlId = pSBMLCompartment->getName();
4537  }
4538  else
4539  {
4540  sbmlId = pSBMLCompartment->getId();
4541  }
4542 
4543  if (sbmlId == pNode->getName())
4544  {
4545  pNode->setName(dynamic_cast<CCompartment*>(it->first)->getObject(CCopasiObjectName("Reference=InitialVolume"))->getCN().c_str());
4546  found = true;
4547  }
4548 
4549  break;
4550 
4551  case SBML_SPECIES:
4552  pSBMLSpecies = dynamic_cast<Species*>(it->second);
4553 
4554  if (this->mLevel == 1)
4555  {
4556  sbmlId = pSBMLSpecies->getName();
4557  }
4558  else
4559  {
4560  sbmlId = pSBMLSpecies->getId();
4561  }
4562 
4563  if (sbmlId == pNode->getName())
4564  {
4565  pNode->setName(dynamic_cast<CMetab*>(it->first)->getObject(CCopasiObjectName("Reference=InitialConcentration"))->getCN().c_str());
4566  found = true;
4567  }
4568 
4569  break;
4570 
4571  case SBML_REACTION:
4572  pSBMLReaction = dynamic_cast<Reaction*>(it->second);
4573 
4574  if (this->mLevel == 1)
4575  {
4576  sbmlId = pSBMLReaction->getName();
4577  }
4578  else
4579  {
4580  sbmlId = pSBMLReaction->getId();
4581  }
4582 
4583  if (sbmlId == pNode->getName())
4584  {
4585  pNode->setName(dynamic_cast<CReaction*>(it->first)->getObject(CCopasiObjectName("Reference=ParticleFlux"))->getCN().c_str());
4586  found = true;
4587  }
4588 
4589  break;
4590 
4591  case SBML_PARAMETER:
4592  pSBMLParameter = dynamic_cast<Parameter*>(it->second);
4593 
4594  if (this->mLevel == 1)
4595  {
4596  sbmlId = pSBMLParameter->getName();
4597  }
4598  else
4599  {
4600  sbmlId = pSBMLParameter->getId();
4601  }
4602 
4603  if (sbmlId == pNode->getName())
4604  {
4605  pNode->setName(dynamic_cast<CModelValue*>(it->first)->getValueReference()->getCN().c_str());
4606  found = true;
4607  }
4608 
4609  break;
4610 
4611  default:
4612  break;
4613  }
4614 
4615  ++it;
4616  }
4617 
4618  if (!found) success = false;
4619  }
4620  }
4621 
4622  for (i = 0; i < iMax; ++i)
4623  {
4624  if (!this->sbmlId2CopasiCN(pNode->getChild(i), copasi2sbmlmap, pParamGroup))
4625  {
4626  success = false;
4627  break;
4628  }
4629  }
4630 
4631  return success;
4632 }
4633 
4634 #ifdef COPASI_DEBUG
4635 void SBMLImporter::printMap(const std::map<CCopasiObject*, SBase*> & copasi2sbml)
4636 {
4637  std::map<CCopasiObject*, SBase*>::const_iterator it = copasi2sbml.begin();
4638  std::map<CCopasiObject*, SBase*>::const_iterator end = copasi2sbml.end();
4639  std::cout << "Number of elements: " << copasi2sbml.size() << std::endl;
4640 
4641  while (it != end)
4642  {
4643  std::cout << "(@" << it->first << ")" << it->first->getObjectName() << " : " << "(@" << it->second << ")" << it->second->getTypeCode() << std::endl;
4644  ++it;
4645  }
4646 
4647  std::cout << std::endl;
4648 }
4649 #endif // COPASI_DEBUG
4650 
4652 {
4653  // set all the old sbml ids
4654  std::map<CFunction*, std::string>::iterator it = this->sbmlIdMap.begin();
4655  std::map<CFunction*, std::string>::iterator endIt = this->sbmlIdMap.end();
4656 
4657  while (it != endIt)
4658  {
4659  it->first->setSBMLId(it->second);
4660  ++it;
4661  }
4662 
4663  // remove all the functions that were added during import
4664  std::set<std::string>::iterator it2 = this->mUsedFunctions.begin();
4665  std::set<std::string>::iterator endIt2 = this->mUsedFunctions.end();
4666 
4667  while (it2 != endIt2)
4668  {
4669  CEvaluationTree* pTree = this->functionDB->findFunction(*it2);
4670  assert(pTree);
4671 
4672  if (pTree->getType() == CEvaluationTree::UserDefined)
4673  {
4674  this->functionDB->removeFunction(pTree->getKey());
4675  }
4676 
4677  ++it2;
4678  }
4679 }
4680 
4681 void SBMLImporter::preprocessNode(ConverterASTNode* pNode, Model* pSBMLModel, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap, Reaction* pSBMLReaction)
4682 {
4683  // this function goes through the tree several times.
4684  // this can probably be handled more intelligently
4685  // maybe we should use some kind of visitor schema
4686  //
4687  // check if there are units on pure numbers
4688  // so that we can create a warning that they have been ignored
4689 
4690  if (this->mLevel > 2 && !this->mUnitOnNumberFound)
4691  {
4692  this->mUnitOnNumberFound = SBMLImporter::checkForUnitsOnNumbers(pNode);
4693  }
4694 
4695  // first replace the calls to explicitly time dependent functions
4696  this->replaceTimeDependentFunctionCalls(pNode);
4697 
4698  if (!this->mDelayFound || pSBMLReaction != NULL)
4699  {
4700  bool result = isDelayFunctionUsed(pNode);
4701 
4702  if (pSBMLReaction != NULL && result)
4703  {
4704  // if this is the first delay we find, we have to populate the id set
4705  // in order to be able to create new global parameters with unique ids
4706  if (!this->mUsedSBMLIdsPopulated)
4707  {
4708  std::map<std::string, const SBase*> idMap;
4709  std::map<std::string, const SBase*> metaIdMap;
4710  // this is overkill, but better than writing yet another function to
4711  // collect ids
4712  SBMLUtils::collectIds(pSBMLModel, idMap, metaIdMap);
4713  std::map<std::string, const SBase*>::iterator it = idMap.begin(), endit = idMap.end();
4714 
4715  while (it != endit)
4716  {
4717  this->mUsedSBMLIds.insert(it->first);
4718  ++it;
4719  }
4720 
4721  this->mUsedSBMLIdsPopulated = true;
4723  }
4724 
4725  // we need a map to store the replacements for local parameters
4726  // which occur within delay calls
4727  std::map<std::string, std::string> replacementMap;
4728  this->replace_delay_nodes(pNode, pSBMLModel, copasi2sbmlmap, pSBMLReaction, replacementMap);
4729 
4730  if (!replacementMap.empty())
4731  {
4732  // replace all local parameter nodes which have been converted to global parameters
4733  // and that were not in delay calls
4734  this->replace_name_nodes(pNode, replacementMap);
4735  // delete the local parameters that have been replaced from the
4736  // reaction
4737  std::map<std::string, std::string>::const_iterator it = replacementMap.begin(), endit = replacementMap.end();
4738  // starting with SBML Level 3, the parameters of a kinetic law are expressed in
4739  // terms of a new class called LocalParameter instead of Parameter
4740  // Unfortunatelly libsbml 4.1 uses separate data structures for
4741  // the Parameters and the LocalParameters which mandates a small
4742  // code change to be on the safe side
4743  ListOfParameters* pList = NULL;
4744 
4745  if (this->mLevel > 2)
4746  {
4747  pList = pSBMLReaction->getKineticLaw()->getListOfLocalParameters();
4748  }
4749  else
4750  {
4751  pList = pSBMLReaction->getKineticLaw()->getListOfParameters();
4752  }
4753 
4754  Parameter* pParam = NULL;
4755 
4756  while (it != endit)
4757  {
4758  pParam = pList->remove(it->first);
4759  assert(pParam != NULL);
4760  pdelete(pParam);
4761  ++it;
4762  }
4763 
4764  this->mReactionsWithReplacedLocalParameters.insert(pSBMLReaction->getId());
4765  }
4766  }
4767 
4768  this->mDelayFound = result;
4769  }
4770 
4771  this->replaceCallNodeNames(pNode);
4772  this->replaceTimeAndAvogadroNodeNames(pNode);
4773 
4774  if (pSBMLReaction != NULL && !this->mSubstanceOnlySpecies.empty())
4775  {
4777  }
4778 
4779  if (!this->mSubstanceOnlySpecies.empty() && this->mpCopasiModel->getQuantityUnitEnum() != CModel::number && pSBMLReaction == NULL)
4780  {
4781  this->replaceAmountReferences(pNode, pSBMLModel, this->mpCopasiModel->getQuantity2NumberFactor(), copasi2sbmlmap);
4782  }
4783 }
4784 
4785 /**
4786  * This method replaces references to the id of species which have the
4787  * hasOnlySubstanceUnits flag set with the reference divided by avogadros
4788  * number.
4789  */
4790 void SBMLImporter::replaceAmountReferences(ConverterASTNode* pASTNode, Model* pSBMLModel, double factor, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
4791 {
4792  CNodeIterator< ConverterASTNode > itNode(pASTNode);
4794 
4795  while (itNode.next() != itNode.end())
4796  {
4797  if (*itNode == NULL)
4798  {
4799  continue;
4800  }
4801 
4802  // check if there is already a multiplication with the avogadro number
4803  // if so, we just replace that multiplication instead of adding an additional division
4804  //
4805  // we only do this, if there are two children to the multiplication because otherwise we would have to do sophisticated checks
4806  // whether the multiplication is really only between a certain species and the avogadro number
4807  if (itNode->getType() == AST_TIMES && itNode->getNumChildren() == 2)
4808  {
4809  std::set<const Parameter*>::const_iterator avoIt = this->mPotentialAvogadroNumbers.begin();
4810  std::set<const Parameter*>::const_iterator avoEndit = this->mPotentialAvogadroNumbers.end();
4811  // check if one of the potantial avogrador numbers is a child to this multiplication
4812  ASTNode* pChild1 = itNode->getChild(0);
4813  ASTNode* pChild2 = itNode->getChild(1);
4814  assert(pChild1 != NULL);
4815  assert(pChild2 != NULL);
4816 
4817  if (pChild1->getType() == AST_NAME && pChild2->getType() == AST_NAME)
4818  {
4819  std::string id1 = pChild1->getName();
4820  std::string id2 = pChild2->getName();
4821  ASTNode *pAvogadro = NULL, *pSpecies = NULL;
4822 
4823  if (pChild1 != NULL && pChild2 != NULL)
4824  {
4825  while (avoIt != avoEndit)
4826  {
4827  if (id1 == (*avoIt)->getId())
4828  {
4829  pAvogadro = pChild1;
4830  pSpecies = pChild2;
4831  // store the id of the species in id1 for use below
4832  id1 = id2;
4833  break;
4834  }
4835  else if (id2 == (*avoIt)->getId())
4836  {
4837  pAvogadro = pChild2;
4838  pSpecies = pChild1;
4839  break;
4840  }
4841 
4842  ++avoIt;
4843  }
4844 
4845  if (pAvogadro != NULL)
4846  {
4847  // check if the potential speices node is really a substance only species
4848  //
4849  std::map<Species*, Compartment*>::const_iterator it = this->mSubstanceOnlySpecies.begin();
4850  std::map<Species*, Compartment*>::const_iterator endit = this->mSubstanceOnlySpecies.end();
4851 
4852  while (it != endit)
4853  {
4854  if (it->first->getId() == id1)
4855  {
4856  // delete the two children
4857  itNode->removeChild(1);
4858  itNode->removeChild(0);
4859  pdelete(pChild1);
4860  pdelete(pChild2);
4861  // now we know that we can change the current node
4862  // to represent the species
4863  itNode->setType(AST_NAME);
4864  itNode->setName(id1.c_str());
4865  itNode.skipChildren();
4866  break;
4867  }
4868 
4869  ++it;
4870  }
4871  }
4872  }
4873  }
4874  else if ((pChild1->getType() == AST_NAME && (pChild2->getType() == AST_REAL || pChild2->getType() == AST_REAL_E)) ||
4875  (pChild2->getType() == AST_NAME && (pChild1->getType() == AST_REAL || pChild1->getType() == AST_REAL_E)))
4876  {
4877  double value = 0.0;
4878  std::string id;
4879 
4880  if (pChild1->getType() == AST_NAME)
4881  {
4882  value = pChild2->getMantissa() * pow(10.0, (double)pChild2->getExponent());
4883  id = pChild1->getName();
4884  }
4885  else
4886  {
4887  value = pChild1->getMantissa() * pow(10.0, (double)pChild1->getExponent());
4888  id = pChild2->getName();
4889  }
4890 
4891  if (areApproximatelyEqual(factor, value, 1e-3))
4892  {
4893  std::map<Species*, Compartment*>::const_iterator it = this->mSubstanceOnlySpecies.begin();
4894  std::map<Species*, Compartment*>::const_iterator endit = this->mSubstanceOnlySpecies.end();
4895 
4896  while (it != endit)
4897  {
4898  if (it->first->getId() == id)
4899  {
4900  // delete the two children
4901  itNode->removeChild(1);
4902  itNode->removeChild(0);
4903  pdelete(pChild1);
4904  pdelete(pChild2);
4905  // now we know that we can change the current node
4906  // to represent the species
4907  itNode->setType(AST_NAME);
4908  itNode->setName(id.c_str());
4909  itNode.skipChildren();
4910  break;
4911  }
4912 
4913  ++it;
4914  }
4915  }
4916  }
4917  }
4918  else if (itNode->getType() == AST_NAME)
4919  {
4920  // check if pNode is a reference to a hasOnlySubstance species
4921  std::string id = itNode->getName();
4922  std::map<Species*, Compartment*>::const_iterator it = this->mSubstanceOnlySpecies.begin();
4923  std::map<Species*, Compartment*>::const_iterator endit = this->mSubstanceOnlySpecies.end();
4924 
4925  while (it != endit)
4926  {
4927  if (it->first->getId() == id)
4928  {
4929  break;
4930  }
4931 
4932  ++it;
4933  }
4934 
4935  if (it != endit)
4936  {
4937  // replace pNode by a division by the quantity to number
4938  // factor
4939  if (this->mPotentialAvogadroNumbers.empty())
4940  {
4941  this->createHasOnlySubstanceUnitFactor(pSBMLModel, factor, copasi2sbmlmap);
4942  }
4943 
4944  itNode->setType(AST_DIVIDE);
4945  itNode->setCharacter('/');
4946  ConverterASTNode* pChild = new ConverterASTNode(AST_NAME);
4947  pChild->setName(id.c_str());
4948  itNode->addChild(pChild);
4949  id = (*this->mPotentialAvogadroNumbers.begin())->getId();
4950  pChild = new ConverterASTNode(AST_NAME);
4951  pChild->setName(id.c_str());
4952  itNode->addChild(pChild);
4953  }
4954  }
4955  }
4956 }
4957 
4959 {
4960  bool result = false;
4961  CNodeIterator< ConverterASTNode > itNode(pASTNode);
4962 
4963  while (itNode.next() != itNode.end())
4964  {
4965  if (*itNode == NULL)
4966  {
4967  continue;
4968  }
4969 
4970  if (itNode->getType() == AST_FUNCTION_DELAY)
4971  {
4972  result = true;
4973  break;
4974  }
4975  }
4976 
4977  return result;
4978 }
4979 
4981 {
4982  CNodeIterator< ASTNode > itNode(pASTNode);
4983 
4984  while (itNode.next() != itNode.end())
4985  {
4986  if (*itNode == NULL)
4987  {
4988  continue;
4989  }
4990 
4991  if (itNode->getType() == AST_NAME_TIME)
4992  {
4993  itNode->setName(this->mpCopasiModel->getObject(CCopasiObjectName("Reference=Time"))->getCN().c_str());
4994  }
4995  else if (itNode->getType() == AST_NAME_AVOGADRO)
4996  {
4997  itNode->setName(this->mpCopasiModel->getObject(CCopasiObjectName("Reference=Avogadro Constant"))->getCN().c_str());
4998  }
4999  }
5000 }
5001 
5002 void SBMLImporter::replaceCallNodeNames(ASTNode* pASTNode)
5003 {
5004  CNodeIterator< ASTNode > itNode(pASTNode);
5005 
5006  while (itNode.next() != itNode.end())
5007  {
5008  if (*itNode == NULL)
5009  {
5010  continue;
5011  }
5012 
5013  if (itNode->getType() == AST_FUNCTION)
5014  {
5015  std::map<std::string, std::string>::const_iterator pos = this->mFunctionNameMapping.find(itNode->getName());
5016 
5017  std::map<std::string, std::string>::const_iterator knownPos = mKnownCustomUserDefinedFunctions.find(itNode->getName());
5018 
5019  if (pos == this->mFunctionNameMapping.end())
5020  {
5021  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 47, itNode->getName());
5022  }
5023 
5024  if (knownPos != mKnownCustomUserDefinedFunctions.end() &&
5025  (knownPos->second == "RATE" && itNode->getNumChildren() == 1))
5026  {
5027  // replace the known function with the correct calls.
5028  // replace rateOf with the rate of whatever is inside
5029  std::string symbol = itNode->getChild(0)->getName();
5030  itNode->removeChild(0);
5031  itNode->setType(AST_NAME);
5032  itNode->setName(symbol.c_str());
5033  itNode->setUserData(strdup("RATE"));
5034  }
5035  else
5036  {
5037  std::string newName = pos->second;
5038  itNode->setName(newName.c_str());
5039  this->mUsedFunctions.insert(newName);
5040  }
5041  }
5042  }
5043 }
5044 
5045 /**
5046  * The methods gets a function where all the parameters have a usage of "PARAMETER".
5047  * In addition it get the root node of a call to that function which is an expression
5048  * and contains the actual objects with which the function is called in a certain reaction.
5049  * From this expression we can determine if there already is a function in the database
5050  * that does the same. Or we can find out if this function is a Mass Action kinetic.
5051  */
5052 CFunction* SBMLImporter::findCorrespondingFunction(const CExpression * pExpression, const CReaction* pCopasiReaction)
5053 {
5054  // We are certain at this point that pExpression is a simple function call.
5055 
5056  // Check whether a function with the name exists in COPASI and whether it is suitable
5057  CFunction* pCorrespondingFunction = NULL;
5058 
5059  std::string Name = pExpression->getRoot()->getData();
5060  pCorrespondingFunction = functionDB->findFunction(Name);
5061 
5062  if (pCorrespondingFunction != NULL)
5063  {
5064  const CFunctionParameters & Variables = pCorrespondingFunction->getVariables();
5065 
5066  for (size_t i = 0; i < Variables.size(); ++i)
5067  {
5068  // This fails for global parameters but those are handled in a second attempt
5069  if (pCopasiReaction->getParameterIndex(Variables[i]->getObjectName()) == C_INVALID_INDEX)
5070  {
5071  pCorrespondingFunction = NULL;
5072  break;
5073  }
5074  }
5075  }
5076 
5077  if (pCorrespondingFunction != NULL)
5078  {
5079  return pCorrespondingFunction;
5080  }
5081 
5082  // TODO This has never worked
5083 #ifdef XXXX
5084  // Check for suitable functions which have a different name but have the same call tree as
5085  std::vector< CFunction * > functions = functionDB->suitableFunctions(pCopasiReaction->getChemEq().getSubstrates().size(),
5086  pCopasiReaction->getChemEq().getProducts().size(),
5087  pCopasiReaction->isReversible() ? TriTrue : TriFalse);
5088  size_t i, iMax = functions.size();
5089 
5090  for (i = 0; i < iMax; ++i)
5091  {
5092  pCorrespondingFunction = (functions[i]);
5093 
5094  // make sure the function is not compared to itself since it can already
5095  // be in the database if it has been used a call in another function
5096  // don't compare the mass action kinetics
5097  if ((pCorrespondingFunction != tree) && (!dynamic_cast<CMassAction*>(pCorrespondingFunction)) && areEqualFunctions(pCorrespondingFunction, tree))
5098  {
5099  const CFunctionParameters & Variables = pFun->getVariables();
5100 
5101  for (size_t i = 0; i < Variables.size(); ++i)
5102  {
5103  if (pCopasiReaction->getParameterIndex(Variables[i]->getObjectName()) == C_INVALID_INDEX)
5104  {
5105  pCorrespondingFunction = NULL;
5106  break;
5107  }
5108  }
5109 
5110  pCorrespondingFunction = pFun;
5111  break;
5112  }
5113  }
5114 
5115 #endif // XXXX
5116 
5117  return pCorrespondingFunction;
5118 }
5119 
5120 // static
5121 bool SBMLImporter::areEqualFunctions(const CFunction* pFun, const CFunction* pFun2)
5122 {
5123  bool result = true;
5124  const CFunctionParameters& funParams1 = pFun->getVariables();
5125  const CFunctionParameters& funParams2 = pFun2->getVariables();
5126 
5127  if (funParams1.size() == funParams2.size())
5128  {
5129  size_t i, iMax = funParams1.size();
5130 
5131  for (i = 0; i < iMax; ++i)
5132  {
5133  const CFunctionParameter* pFunParam1 = funParams1[i];
5134  const CFunctionParameter* pFunParam2 = funParams2[i];
5135 
5136  if (pFunParam1->getObjectName() != pFunParam2->getObjectName())
5137  {
5138  result = false;
5139  break;
5140  }
5141  }
5142 
5143  if (result == true)
5144  {
5145  const CEvaluationNode* pNodeFun1 = static_cast<const CEvaluationNode*>(pFun->getRoot());
5146  const CEvaluationNode* pNodeFun2 = static_cast<const CEvaluationNode*>(pFun2->getRoot());
5147  result = areEqualSubtrees(pNodeFun1, pNodeFun2);
5148  }
5149  }
5150  else
5151  {
5152  result = false;
5153  }
5154 
5155  return result;
5156 }
5157 
5158 // static
5160 {
5161  // TODO CRITICAL We need to use a node iterator
5162  bool result = ((pNode1->getType() == pNode2->getType()) && (pNode1->getData() == pNode2->getData()));
5163  const CEvaluationNode* pChild1 = static_cast<const CEvaluationNode*>(pNode1->getChild());
5164  const CEvaluationNode* pChild2 = static_cast<const CEvaluationNode*>(pNode2->getChild());
5165 
5166  while (result && pChild1 && pChild2)
5167  {
5168  result = areEqualSubtrees(pChild1, pChild2);
5169  pChild1 = static_cast<const CEvaluationNode*>(pChild1->getSibling());
5170  pChild2 = static_cast<const CEvaluationNode*>(pChild2->getSibling());
5171  }
5172 
5173  result = (result && !pChild1 && !pChild2);
5174  return result;
5175 }
5176 
5177 std::vector<CEvaluationNodeObject*>* SBMLImporter::isMassAction(const CEvaluationTree* pTree, const CChemEq& chemicalEquation, const CEvaluationNodeCall* pCallNode)
5178 {
5179  // TODO CRITICAL We need to use a node iterator
5180  CEvaluationTree::Type type = pTree->getType();
5181  std::vector< std::vector< std::string > > functionArgumentCNs;
5182  const CEvaluationNode* pChildNode = NULL;
5183  std::string str;
5184  std::vector<CEvaluationNodeObject*>* result = NULL;
5185 
5186  switch (type)
5187  {
5191  pChildNode = static_cast<const CEvaluationNode*>(pCallNode->getChild());
5192 
5193  while (pChildNode)
5194  {
5195  if (pChildNode->getType() == CEvaluationNode::OBJECT)
5196  {
5197  str = pChildNode->getData().substr(1, pChildNode->getData().length() - 2);
5198  functionArgumentCNs.push_back(std::vector<std::string>());
5199  functionArgumentCNs[functionArgumentCNs.size() - 1].push_back(str);
5200  pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
5201  }
5202  else
5203  {
5204  fatalError();
5205  }
5206  }
5207 
5208  result = this->isMassActionFunction(dynamic_cast<const CFunction*>(pTree), chemicalEquation, functionArgumentCNs);
5209  break;
5210 
5212  result = this->isMassActionExpression(pTree->getRoot(), chemicalEquation);
5213  break;
5214 
5215  default:
5216  fatalError();
5217  break;
5218  }
5219 
5220  return result;
5221 }
5222 
5223 std::vector<CEvaluationNodeObject*>* SBMLImporter::isMassActionExpression(const CEvaluationNode* pRootNode, const CChemEq& chemicalEquation)
5224 {
5225  bool result = true;
5226  std::vector<CEvaluationNodeObject*>* v = NULL;
5227 
5228  if (chemicalEquation.getReversibility())
5229  {
5230 #ifdef COPASI_DEBUG
5231  CEvaluationNode* pTmpNode = pRootNode->copyBranch();
5232 // CEvaluationNode* pTmpNode = CEvaluationNodeNormalizer::normalize(pRootNode);
5233 #else
5234  CEvaluationNode* pTmpNode = pRootNode->copyBranch();
5235 #endif /* COPASI_DEBUG */
5236  assert(pTmpNode != NULL);
5237  // the root node must be a minus operator
5238  // the two children must be irreversible mass action terms
5240 
5241  if (result)
5242  {
5243  const CEvaluationNode* pChildNode = static_cast<const CEvaluationNode*>(pTmpNode->getChild());
5244  result = (pChildNode != NULL);
5245 
5246  if (result)
5247  {
5248  CChemEq tmpEq;
5249  const CCopasiVector<CChemEqElement>* metabolites = &chemicalEquation.getSubstrates();
5250  size_t i, iMax = metabolites->size();
5251  result = (iMax > 0);
5252 
5253  if (result)
5254  {
5255  for (i = 0; i < iMax; ++i)
5256  {
5257  const CChemEqElement* element = (*metabolites)[i];
5258  tmpEq.addMetabolite(element->getMetaboliteKey(), element->getMultiplicity(), CChemEq::SUBSTRATE);
5259  }
5260 
5261  v = this->isMassActionExpression(pChildNode, tmpEq);
5262 
5263  if (!v)
5264  {
5265  fatalError();
5266  }
5267 
5268  result = !v->empty();
5269 
5270  if (result)
5271  {
5272  pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
5273  result = (pChildNode != NULL);
5274 
5275  if (result)
5276  {
5277  CChemEq tmpEq2;
5278  metabolites = &chemicalEquation.getProducts();
5279  iMax = metabolites->size();
5280  result = (iMax > 0);
5281 
5282  for (i = 0; i < iMax; ++i)
5283  {
5284  const CChemEqElement* element = (*metabolites)[i];
5285  tmpEq2.addMetabolite(element->getMetaboliteKey(), element->getMultiplicity(), CChemEq::SUBSTRATE);
5286  }
5287 
5288  std::vector<CEvaluationNodeObject*>* v2 = this->isMassActionExpression(pChildNode, tmpEq2);
5289 
5290  if (!v2)
5291  {
5292  fatalError();
5293  }
5294 
5295  result = !v2->empty();
5296 
5297  if (result)
5298  {
5299  v->push_back((*v2)[0]);
5300  }
5301  else
5302  {
5303  v->clear();
5304  }
5305 
5306  pdelete(v2);
5307  }
5308  }
5309  }
5310  else
5311  {
5312  v = new std::vector<CEvaluationNodeObject*>;
5313  }
5314  }
5315  }
5316  else
5317  {
5318  v = new std::vector<CEvaluationNodeObject*>;
5319  }
5320 
5321  pdelete(pTmpNode);
5322  }
5323  else
5324  {
5325  // the expression must contain exactly one global or local parameter
5326  // the expression must contain each substrate in the CChemicalReaction
5327  std::vector<const CEvaluationNode*> arguments;
5328  std::map<const CMetab*, C_FLOAT64> multiplicityMap;
5329  this->separateProductArguments(pRootNode, arguments);
5330  size_t numParameters = 0, i, iMax = arguments.size();
5331  v = new std::vector<CEvaluationNodeObject*>;
5332 
5333  if (iMax != 0)
5334  {
5335  std::vector<CCopasiContainer*> listOfContainers;
5336  listOfContainers.push_back(this->mpCopasiModel);
5337 
5338  for (i = 0; (i < iMax) && (numParameters < 2); ++i)
5339  {
5340  const CEvaluationNode* pNode = arguments[i];
5341 
5342  // the node can either be an object node
5343  // or it can be a power function node
5344  if (pNode->getType() == CEvaluationNode::OBJECT)
5345  {
5346  // it can be a global or a local parameter or an metabolite
5347  std::string objectCN = pNode->getData().substr(1, pNode->getData().length() - 2);
5348  const CCopasiObject* pObject = mpDataModel->ObjectFromName(listOfContainers, objectCN);
5349 
5350  if (!pObject)
5351  {
5352  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 39, objectCN.c_str());
5353  }
5354 
5355  if (pObject->isReference())
5356  {
5357  pObject = pObject->getObjectParent();
5358 
5359  if (dynamic_cast<const CMetab*>(pObject))
5360  {
5361  const CMetab* pMetab = static_cast<const CMetab*>(pObject);
5362 
5363  if (multiplicityMap.find(pMetab) != multiplicityMap.end())
5364  {
5365  multiplicityMap[pMetab] = multiplicityMap[pMetab] + 1.0;
5366  }
5367  else
5368  {
5369  multiplicityMap[pMetab] = 1.0;
5370  }
5371  }
5372  else if (dynamic_cast<const CModelValue*>(pObject))
5373  {
5374  ++numParameters;
5376  }
5377  else
5378  {
5379  result = false;
5380  break;
5381  }
5382  }
5383  else if (dynamic_cast<const CCopasiParameter*>(pObject))
5384  {
5385  ++numParameters;
5387  }
5388  else
5389  {
5390  result = false;
5391  break;
5392  }
5393  }
5395  {
5396  // the two children must be a metabolite node and a number node in this order
5397  const CEvaluationNode* pChildNode = static_cast<const CEvaluationNode*>(pNode->getChild());
5398 
5399  if (pChildNode->getType() == CEvaluationNode::OBJECT)
5400  {
5401  std::string objectCN = pChildNode->getData().substr(1, pChildNode->getData().length() - 2);
5402  const CCopasiObject* pObject = mpDataModel->ObjectFromName(listOfContainers, objectCN);
5403  assert(pObject);
5404 
5405  if (pObject->isReference())
5406  {
5407  pObject = pObject->getObjectParent();
5408 
5409  if (dynamic_cast<const CMetab*>(pObject))
5410  {
5411  pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
5412  assert(pChildNode);
5413 
5414  CEvaluationNode::Type type = CEvaluationNode::type(pChildNode->getType());
5415 
5416  if (type == CEvaluationNode::NUMBER)
5417  {
5418  const CMetab* pMetab = static_cast<const CMetab*>(pObject);
5419 
5420  if (multiplicityMap.find(pMetab) != multiplicityMap.end())
5421  {
5422  multiplicityMap[pMetab] = multiplicityMap[pMetab] + pChildNode->getValue();
5423  }
5424  else
5425  {
5426  multiplicityMap[pMetab] = pChildNode->getValue();
5427  }
5428  }
5429  else if (type == CEvaluationNode::FUNCTION &&
5431  CEvaluationNode::type(((CEvaluationNode*)pChildNode->getChild())->getType()) == CEvaluationNode::NUMBER)
5432  {
5433  const CMetab* pMetab = static_cast<const CMetab*>(pObject);
5434  multiplicityMap[pMetab] = -1 * ((CEvaluationNodeNumber*)(pChildNode->getChild()))->getValue();
5435  }
5436  else
5437  {
5438  // not a math action
5439  }
5440  }
5441  }
5442  else
5443  {
5444  result = false;
5445  break;
5446  }
5447  }
5448  else
5449  {
5450  result = false;
5451  break;
5452  }
5453  }
5454  else
5455  {
5456  result = false;
5457  break;
5458  }
5459  }
5460  }
5461 
5462  if (numParameters != 1)
5463  {
5464  result = false;
5465  }
5466 
5467  if (result)
5468  {
5469  const CCopasiVector<CChemEqElement>& metabolites = chemicalEquation.getSubstrates();
5470  size_t i, iMax = metabolites.size();
5471 
5472  // all metabolites must occur in the muliplicityMap so they have to have the same size
5473  // and a mass action must have at least one metabolite
5474  if (iMax == 0 || iMax != multiplicityMap.size()) result = false;
5475 
5476  for (i = 0; i < iMax && result; ++i)
5477  {
5478  // the metabolite has to be present in the multiplicityMap, otherwise it is not a mass action
5479  // the stoichiometry also has to fit
5480  std::map<const CMetab*, C_FLOAT64>::iterator pos = multiplicityMap.find(metabolites[i]->getMetabolite());
5481 
5482  if (pos == multiplicityMap.end() ||
5483  !areApproximatelyEqual(pos->second, metabolites[i]->getMultiplicity(), 0.01))
5484  {
5485  result = false;
5486  break;
5487  }
5488  }
5489  }
5490 
5491  if (!result)
5492  {
5493  v->clear();
5494  }
5495  }
5496 
5497  return v;
5498 }
5499 
5500 std::vector<CEvaluationNodeObject*>* SBMLImporter::isMassActionFunction(const CFunction* pFun, const CChemEq& chemicalEquation, const std::vector<std::vector< std::string > >& functionArgumentCNs)
5501 {
5502  // create an expression from the function and call isMassActionExpression
5503  CEvaluationTree* pExpressionTree = this->createExpressionFromFunction(pFun, functionArgumentCNs);
5504 
5505  if (!pExpressionTree)
5506  {
5507  return NULL;
5508  }
5509 
5510  std::vector<CEvaluationNodeObject*>* v = this->isMassActionExpression(pExpressionTree->getRoot(), chemicalEquation);
5511 
5512  delete pExpressionTree;
5513  return v;
5514 }
5515 
5516 CEvaluationTree* SBMLImporter::createExpressionFromFunction(const CFunction* pFun, const std::vector<std::vector<std::string > >& functionArgumentCNs)
5517 {
5518  CEvaluationTree* pTree = NULL;
5519  const CFunctionParameters& pFunParams = pFun->getVariables();
5520  std::string str;
5521 
5522  if (pFunParams.size() == functionArgumentCNs.size())
5523  {
5524  std::map<std::string , std::string> variable2CNMap;
5525  size_t i, iMax = pFunParams.size();
5526 
5527  for (i = 0; i < iMax; ++i)
5528  {
5529  // vectors should not occur here
5530  assert(functionArgumentCNs[i].size() == 1);
5531  variable2CNMap[pFunParams[i]->getObjectName()] = functionArgumentCNs[i][0];
5532  }
5533 
5534  CEvaluationNode* pTmpNode = this->variables2objects(pFun->getRoot(), variable2CNMap);
5535  assert(pTmpNode);
5537  pTree->setRoot(pTmpNode);
5538  }
5539 
5540  return pTree;
5541 }
5542 
5543 void SBMLImporter::separateProductArguments(const CEvaluationNode* pRootNode, std::vector<const CEvaluationNode*>& arguments)
5544 {
5545  // TODO CRITICAL We need to use a node iterator
5546  const CEvaluationNodeOperator* pMultiplyNode = dynamic_cast<const CEvaluationNodeOperator*>(pRootNode);
5547 
5549  {
5550  // check if one if the children is an object node or a power operator, if so,
5551  // add the node to the vector
5552  // the nodes not one of those two are passed to this function recursively.
5553  const CEvaluationNode* pChildNode = static_cast<const CEvaluationNode*>(pMultiplyNode->getChild());
5554 
5555  while (pChildNode)
5556  {
5557  const CEvaluationNodeObject* pObjectNode = dynamic_cast<const CEvaluationNodeObject*>(pChildNode);
5558  const CEvaluationNodeOperator* pOperatorNode = dynamic_cast<const CEvaluationNodeOperator*>(pChildNode);
5559 
5560  if (pObjectNode)
5561  {
5562  arguments.push_back(pObjectNode);
5563  }
5564  else if (pOperatorNode && (((CEvaluationNodeOperator::SubType)CEvaluationNode::subType(pOperatorNode->getType())) == CEvaluationNodeOperator::POWER))
5565  {
5566  arguments.push_back(pOperatorNode);
5567  }
5568  else
5569  {
5570  this->separateProductArguments(pChildNode, arguments);
5571 
5572  if (arguments.empty())
5573  {
5574  // it is not a mass action kinetic, so we can stop here
5575  break;
5576  }
5577  }
5578 
5579  pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
5580  }
5581  }
5582  else
5583  {
5584  arguments.clear();
5585  }
5586 }
5587 
5588 void SBMLImporter::setCorrectUsage(CReaction* pCopasiReaction, const CEvaluationNodeCall* pCallNode)
5589 {
5590  // find out what type each argument in the call node has.
5591  // it can be a local parameter, a global parameter, a compartment or a metabolite
5592  // if it is a metabolite, try to find out if it is a substrate, product, or modifier
5593  if (!pCallNode)
5594  {
5595  fatalError();
5596  }
5597 
5598  const CEvaluationNode* pChildNode = static_cast<const CEvaluationNode*>(pCallNode->getChild());
5599 
5600  std::vector<CCopasiContainer*> listOfContainers = std::vector<CCopasiContainer*>();
5601 
5602  listOfContainers.push_back(this->mpCopasiModel);
5603 
5604  CChemEq& pChemEq = pCopasiReaction->getChemEq();
5605 
5606  std::map<const CChemEqElement*, CChemEq::MetaboliteRole> v;
5607 
5608  const CCopasiVector<CChemEqElement>* pV = &pChemEq.getSubstrates();
5609 
5610  size_t i, iMax = pV->size();
5611 
5612  for (i = 0; i < iMax; ++i)
5613  {
5614  v[(*pV)[i]] = CChemEq::SUBSTRATE;
5615  }
5616 
5617  pV = &pChemEq.getProducts();
5618  iMax = pV->size();
5619 
5620  for (i = 0; i < iMax; ++i)
5621  {
5622  v[(*pV)[i]] = CChemEq::PRODUCT;
5623  }
5624 
5625  pV = &pChemEq.getModifiers();
5626  iMax = pV->size();
5627 
5628  for (i = 0; i < iMax; ++i)
5629  {
5630  v[(*pV)[i]] = CChemEq::MODIFIER;
5631  }
5632 
5633  size_t parameterIndex = 0;
5634 
5635  while (pChildNode)
5636  {
5637  const CEvaluationNodeObject* pObjectNode = dynamic_cast<const CEvaluationNodeObject*>(pChildNode);
5638 
5639  if (!pObjectNode)
5640  {
5641  fatalError();
5642  }
5643 
5644  const CCopasiObject* object = mpDataModel->ObjectFromName(listOfContainers, pObjectNode->getData().substr(1, pObjectNode->getData().length() - 2));
5645 
5646  if (!object)
5647  {
5648  fatalError();
5649  }
5650 
5651  CFunctionParameter* pFunParam = const_cast<CFunctionParameter*>(pCopasiReaction->getFunction()->getVariables()[parameterIndex]);
5652 
5653  if (dynamic_cast<const CCopasiObjectReference<C_FLOAT64>*>(object))
5654  {
5655  object = object->getObjectParent();
5656 
5657  if (dynamic_cast<const CMetab*>(object))
5658  {
5659  std::map<const CChemEqElement*, CChemEq::MetaboliteRole>::iterator it = v.begin();
5660  std::map<const CChemEqElement*, CChemEq::MetaboliteRole>::iterator endIt = v.end();
5661 
5662  while (it != endIt)
5663  {
5664  if (it->first == object)
5665  {
5666  // get the role of the metabolite
5667  switch (it->second)
5668  {
5669  case CChemEq::SUBSTRATE:
5670  // it is a substrate
5672  break;
5673 
5674  case CChemEq::PRODUCT:
5675  // it is a product
5677  break;
5678 
5679  case CChemEq::MODIFIER:
5680  // it is a modifier
5682  break;
5683 
5684  default:
5685  fatalError();
5686  break;
5687  }
5688  }
5689 
5690  ++it;
5691  }
5692 
5693  if (it == endIt)
5694  {
5695  fatalError();
5696  }
5697  }
5698  else if (dynamic_cast<const CModelValue*>(object))
5699  {
5700  // it is a global parameter
5702  }
5703  else if (dynamic_cast<const CCompartment*>(object))
5704  {
5705  // it is a volume
5707  }
5708  else
5709  {
5710  fatalError()
5711  }
5712  }
5713  else
5714  {
5715  // it is a local parameter
5717  }
5718 
5719  pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
5720  ++parameterIndex;
5721  }
5722 }
5723 
5724 bool containsKey(const CCopasiVector < CChemEqElement >& list, const std::string& key)
5725 {
5728 
5729  for (; it != end; ++it)
5730  {
5731  if ((*it)->getMetaboliteKey() == key)
5732  {
5733  return true;
5734  }
5735  }
5736 
5737  return false;
5738 }
5739 
5740 void SBMLImporter::doMapping(CReaction* pCopasiReaction, const CEvaluationNodeCall* pCallNode)
5741 {
5742  // map the first argument of the call node to the first variable of the function of the reaction
5743  // and so on
5744  if (!pCallNode)
5745  {
5746  fatalError();
5747  }
5748 
5749  std::vector<CCopasiContainer*> listOfContainers;
5750  listOfContainers.push_back(this->mpCopasiModel);
5751 
5752  if (dynamic_cast<const CMassAction*>(pCopasiReaction->getFunction()))
5753  {
5754  const CEvaluationNodeObject* pChild = dynamic_cast<const CEvaluationNodeObject*>(pCallNode->getChild());
5755  std::string objectCN = pChild->getData();
5756  objectCN = objectCN.substr(1, objectCN.length() - 2);
5757  CCopasiObject* pObject = mpDataModel->ObjectFromName(listOfContainers, objectCN);
5758 
5759  if (!pObject)
5760  {
5761  fatalError();
5762  }
5763 
5764  if (pObject->isReference())
5765  {
5766  pObject = pObject->getObjectParent();
5767  }
5768 
5769  const std::string& objectKey = pObject->getKey();
5770 
5771  pCopasiReaction->setParameterMapping("k1", objectKey);
5772 
5773  const CCopasiVector<CChemEqElement>* metabolites = &pCopasiReaction->getChemEq().getSubstrates();
5774 
5775  size_t i, iMax = metabolites->size();
5776 
5777  size_t j, jMax;
5778 
5779  for (i = 0; i < iMax; ++i)
5780  for (j = 0, jMax = static_cast<int>(fabs((*metabolites)[i]->getMultiplicity())); j < jMax; j++)
5781  pCopasiReaction->addParameterMapping("substrate", (*metabolites)[i]->getMetaboliteKey());
5782 
5783  if (pCopasiReaction->isReversible())
5784  {
5785  pChild = dynamic_cast<const CEvaluationNodeObject*>(pChild->getSibling());
5786  std::string objectCN = pChild->getData();
5787  objectCN = objectCN.substr(1, objectCN.length() - 2);
5788  CCopasiObject* pObject = mpDataModel->ObjectFromName(listOfContainers, objectCN);
5789 
5790  if (!pObject)
5791  {
5792  fatalError();
5793  }
5794 
5795  if (pObject->isReference())
5796  {
5797  pObject = pObject->getObjectParent();
5798  }
5799 
5800  const std::string& objectKey = pObject->getKey();
5801 
5802  pCopasiReaction->setParameterMapping("k2", objectKey);
5803 
5804  const CCopasiVector<CChemEqElement>* metabolites = &pCopasiReaction->getChemEq().getProducts();
5805 
5806  iMax = metabolites->size();
5807 
5808  for (i = 0; i < iMax; ++i)
5809  for (j = 0, jMax = static_cast<int>(fabs((*metabolites)[i]->getMultiplicity())); j < jMax; j++)
5810  pCopasiReaction->addParameterMapping("product", (*metabolites)[i]->getMetaboliteKey());
5811  }
5812  }
5813  else
5814  {
5815  const CFunctionParameters & Variables = pCopasiReaction->getFunction()->getVariables();
5816  size_t i, iMax = Variables.size();
5817  const CEvaluationNodeObject* pChild = dynamic_cast<const CEvaluationNodeObject*>(pCallNode->getChild());
5818 
5819  for (i = 0; i < iMax; ++i)
5820  {
5821  if (!pChild)
5822  {
5823  fatalError();
5824  }
5825 
5826  std::string objectCN = pChild->getData();
5827  objectCN = objectCN.substr(1, objectCN.length() - 2);
5828  CCopasiObject* pObject = mpDataModel->ObjectFromName(listOfContainers, objectCN);
5829 
5830  if (!pObject)
5831  {
5832  fatalError();
5833  }
5834 
5835  if (pObject->isReference())
5836  {
5837  pObject = pObject->getObjectParent();
5838  }
5839 
5840  const std::string& objectKey = pObject->getKey();
5841 
5842  pCopasiReaction->setParameterMapping(i, objectKey);
5843 
5844  const CChemEq& eqn = pCopasiReaction->getChemEq();
5845 
5846  bool reversible = eqn.getReversibility();
5847 
5848  // We guess what the role of a variable of newly imported function is:
5849  if (Variables[i]->getUsage() == CFunctionParameter::VARIABLE)
5850  {
5852 
5853  if (pObject->getObjectType() == "Metabolite")
5854  {
5855  if (containsKey(eqn.getSubstrates(), objectKey))
5857 
5858  if (Role == CFunctionParameter::PARAMETER &&
5859  containsKey(eqn.getProducts(), objectKey))
5860  {
5861  // Fix for Bug 1882: if not reversible, only mark as product if it does
5862  // not appear in the list of modifiers
5863  if (!reversible && containsKey(eqn.getModifiers(), objectKey))
5865  else
5867  }
5868 
5869  // It is not a substrate and not a product therefore we must have a modifier
5870  if (Role == CFunctionParameter::PARAMETER)
5871  {
5873  }
5874  }
5875  else if (pObject->getObjectType() == "Model")
5876  {
5877  Role = CFunctionParameter::TIME;
5878  }
5879  else if (pObject->getObjectType() == "Compartment")
5880  {
5882  }
5883 
5884  const_cast< CFunctionParameter * >(Variables[i])->setUsage(Role);
5885  }
5886 
5887  pChild = dynamic_cast<const CEvaluationNodeObject*>(pChild->getSibling());
5888  }
5889  }
5890 }
5891 
5893 {
5894  // it is a simple function call if it is a CEvaluationNodeCall object and all
5895  // its arguments are object nodes.
5896  bool result = true;
5897 
5898  if (dynamic_cast<const CEvaluationNodeCall*>(pRootNode))
5899  {
5900  const CEvaluationNode* pChildNode = static_cast<const CEvaluationNode*>(pRootNode->getChild());
5901 
5902  // I guess it must have at least one child to qualify.
5903  if (!pChildNode)
5904  {
5905  result = false;
5906  }
5907 
5908  while (pChildNode)
5909  {
5910  if (!dynamic_cast<const CEvaluationNodeObject*>(pChildNode))
5911  {
5912  result = false;
5913  break;
5914  }
5915 
5916  pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
5917  }
5918  }
5919  else
5920  {
5921  result = false;
5922  }
5923 
5924  return result;
5925 }
5926 
5927 ConverterASTNode* SBMLImporter::isMultipliedByVolume(const ASTNode* node, const std::string& compartmentSBMLId)
5928 {
5929  ConverterASTNode* result = NULL;
5930 
5931  if (node->getType() == AST_TIMES || node->getType() == AST_DIVIDE)
5932  {
5933  ConverterASTNode* pTmpResultNode = new ConverterASTNode(node->getType());
5934  unsigned int i, iMax = node->getNumChildren();
5935  bool found = false;
5936 
5937  for (i = 0; i < iMax; ++i)
5938  {
5939  const ASTNode* child = node->getChild(i);
5940 
5941  if (node->getType() == AST_TIMES && child->getType() == AST_NAME && child->getName() == compartmentSBMLId)
5942  {
5943  found = true;
5944  }
5945  else if ((!found) && (child->getType() == AST_TIMES || child->getType() == AST_DIVIDE))
5946  {
5947  ASTNode* pSubResult = this->isMultipliedByVolume(child, compartmentSBMLId);
5948 
5949  if (pSubResult)
5950  {
5951  found = true;
5952 
5953  if (pSubResult->getNumChildren() > 1)
5954  {
5955  pTmpResultNode->addChild(pSubResult);
5956  }
5957  else if (pSubResult->getNumChildren() == 1)
5958  {
5959  pTmpResultNode->addChild(static_cast<ASTNode*>(static_cast<ConverterASTNode*>(pSubResult)->removeChild(0)));
5960  delete pSubResult;
5961  }
5962  else
5963  {
5964  delete pSubResult;
5965  }
5966  }
5967  else
5968  {
5969  pTmpResultNode->addChild(new ConverterASTNode(*child));
5970  }
5971  }
5972  else
5973  {
5974  // bail in case of hOSU
5975  if (child->getType() == AST_NAME && mSubstanceOnlySpecies.size() > 0)
5976  {
5977  std::map<Species*, Compartment*>::iterator it = this->mSubstanceOnlySpecies.begin();
5978  std::map<Species*, Compartment*>::iterator endIt = this->mSubstanceOnlySpecies.end();
5979 
5980  while (it != endIt)
5981  {
5982  if (it->first->getId() == child->getName()) return NULL;
5983 
5984  it++;
5985  }
5986  }
5987 
5988  pTmpResultNode->addChild(new ConverterASTNode(*child));
5989  }
5990  }
5991 
5992  if (found)
5993  {
5994  result = pTmpResultNode;
5995  }
5996  else
5997  {
5998  delete pTmpResultNode;
5999  }
6000  }
6001 
6002  return result;
6003 }
6004 
6005 CEvaluationNode* SBMLImporter::variables2objects(const CEvaluationNode* pOrigNode, const std::map<std::string, std::string>& replacementMap)
6006 {
6007  // TODO CRITICAL We need to use a node iterator
6008  CEvaluationNode* pResultNode = NULL;
6009 
6010  if (dynamic_cast<const CEvaluationNodeVariable*>(pOrigNode))
6011  {
6012  std::map<std::string , std::string>::const_iterator pos = replacementMap.find(pOrigNode->getData());
6013 
6014  if (pos == replacementMap.end()) fatalError();
6015 
6016  pResultNode = new CEvaluationNodeObject(CEvaluationNodeObject::CN, "<" + pos->second + ">");
6017  }
6018  else
6019  {
6020  pResultNode = CEvaluationNode::create(pOrigNode->getType(), pOrigNode->getData());
6021  const CEvaluationNode* pChildNode = static_cast<const CEvaluationNode*>(pOrigNode->getChild());
6022 
6023  while (pChildNode)
6024  {
6025  pResultNode->addChild(this->variables2objects(pChildNode, replacementMap));
6026  pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
6027  }
6028  }
6029 
6030  return pResultNode;
6031 }
6032 
6034 {
6035  std::vector<CCopasiContainer*> v;
6036  v.push_back(this->mpCopasiModel);
6037  CEvaluationNodeObject* pObjectNode = dynamic_cast<CEvaluationNodeObject*>(pCallNode->getChild());
6038  assert(pObjectNode);
6039  CCopasiObjectName objectName = CCopasiObjectName(pObjectNode->getData().substr(1, pObjectNode->getData().length() - 2));
6040  CCopasiObject* pObject = mpDataModel->ObjectFromName(v, objectName);
6041  assert(pObject);
6042 
6043  if (dynamic_cast<CCopasiParameter*>(pObject))
6044  {
6045  pObject->setObjectName("k1");
6046  pObjectNode->setData("<" + pObject->getCN() + ">");
6047  }
6048 
6049  pObjectNode = dynamic_cast<CEvaluationNodeObject*>(pObjectNode->getSibling());
6050 
6051  if (pObjectNode)
6052  {
6053  objectName = CCopasiObjectName(pObjectNode->getData().substr(1, pObjectNode->getData().length() - 2));
6054  pObject = mpDataModel->ObjectFromName(v, objectName);
6055  assert(pObject);
6056 
6057  if (dynamic_cast<CCopasiParameter*>(pObject))
6058  {
6059  pObject->setObjectName("k2");
6060  pObjectNode->setData("<" + pObject->getCN() + ">");
6061  }
6062  }
6063 }
6064 
6065 bool SBMLImporter::containsVolume(const ASTNode* pNode, const std::string& compartmentSBMLId)
6066 {
6067  bool result = false;
6068  unsigned int i, iMax = pNode->getNumChildren();
6069 
6070  for (i = 0; i < iMax; ++i)
6071  {
6072  if (pNode->getChild(i)->getType() == AST_NAME &&
6073  pNode->getChild(i)->getName() == compartmentSBMLId)
6074  {
6075  result = true;
6076  break;
6077  }
6078  }
6079 
6080  return result;
6081 }
6082 
6084 {mpImportHandler = pHandler;}
6085 
6087 {return mpImportHandler;}
6088 
6089 bool SBMLImporter::removeUnusedFunctions(CFunctionDB* pTmpFunctionDB, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
6090 {
6091  if (pTmpFunctionDB)
6092  {
6093  unsigned C_INT32 step = 0, totalSteps = 0;
6094  size_t hStep = C_INVALID_INDEX;
6095  size_t i, iMax = this->mpCopasiModel->getReactions().size();
6096 
6097  if (mpImportHandler)
6098  {
6099  step = 0;
6100  totalSteps =
6101  (unsigned C_INT32)(iMax + this->mpCopasiModel->getCompartments().size() + this->mpCopasiModel->getMetabolites().size() + this->mpCopasiModel->getModelValues().size());
6102  hStep = mpImportHandler->addItem("Searching used functions...",
6103  step,
6104  &totalSteps);
6105  }
6106 
6107  std::set<std::string> functionNameSet;
6108 
6109  for (i = 0; i < iMax; ++i)
6110  {
6111  const CFunction* pTree = this->mpCopasiModel->getReactions()[i]->getFunction();
6112 
6113  // this is a hack, but if we changed stoichiometries in reactions, we can no longer be sure that the
6114  // kinetic function fits the reaction, so in this case it is probably better to set all functions used
6115  // in directions to general
6116  // If we don't do this, the user will not be able to use the function in the reaction in the GUI if he
6117  // changed it once.
6118  // TODO only set those functions to general that really contain modified stoichiometries
6119  if (this->mConversionFactorFound)
6120  {
6121  // TODO another const_cast that should disappear
6122  // TODO because there is probably a better place to do this
6123  // TODO but this means rewriting the code that modifies the
6124  // TODO stoichiometries of reactions and put it all in one place.
6125  const_cast<CFunction*>(pTree)->setReversible(TriUnspecified);
6126  }
6127 
6128  if (functionNameSet.find(pTree->getObjectName()) == functionNameSet.end())
6129  {
6130  functionNameSet.insert(pTree->getObjectName());
6131  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
6132  }
6133 
6134  ++step;
6135 
6136  if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return false;
6137  }
6138 
6139  iMax = this->mpCopasiModel->getCompartments().size();
6140 
6141  for (i = 0; i < iMax; ++i)
6142  {
6143  CModelEntity* pME = this->mpCopasiModel->getCompartments()[i];
6144 
6145  if (pME->getStatus() != CModelEntity::FIXED)
6146  {
6147  const CEvaluationTree* pTree = pME->getExpressionPtr();
6148 
6149  if (pTree != NULL)
6150  {
6151  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
6152  }
6153  }
6154 
6155  if (pME->getStatus() != CModelEntity::ASSIGNMENT)
6156  {
6157  const CEvaluationTree* pTree = pME->getInitialExpressionPtr();
6158 
6159  if (pTree != NULL)
6160  {
6161  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
6162  }
6163  }
6164 
6165  ++step;
6166  }
6167 
6168  iMax = this->mpCopasiModel->getMetabolites().size();
6169 
6170  for (i = 0; i < iMax; ++i)
6171  {
6172  CModelEntity* pME = this->mpCopasiModel->getMetabolites()[i];
6173 
6174  if (pME->getStatus() != CModelEntity::FIXED)
6175  {
6176  const CEvaluationTree* pTree = pME->getExpressionPtr();
6177 
6178  if (pTree != NULL)
6179  {
6180  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
6181  }
6182  }
6183 
6184  if (pME->getStatus() != CModelEntity::ASSIGNMENT)
6185  {
6186  const CEvaluationTree* pTree = pME->getInitialExpressionPtr();
6187 
6188  if (pTree != NULL)
6189  {
6190  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
6191  }
6192  }
6193 
6194  ++step;
6195  }
6196 
6197  iMax = this->mpCopasiModel->getModelValues().size();
6198 
6199  for (i = 0; i < iMax; ++i)
6200  {
6201  CModelEntity* pME = this->mpCopasiModel->getModelValues()[i];
6202 
6203  if (pME->getStatus() != CModelEntity::FIXED)
6204  {
6205  const CEvaluationTree* pTree = pME->getExpressionPtr();
6206 
6207  if (pTree != NULL)
6208  {
6209  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
6210  }
6211  }
6212 
6213  if (pME->getStatus() != CModelEntity::ASSIGNMENT)
6214  {
6215  const CEvaluationTree* pTree = pME->getInitialExpressionPtr();
6216 
6217  if (pTree != NULL)
6218  {
6219  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
6220  }
6221  }
6222 
6223  ++step;
6224  }
6225 
6226  // find the used function in events
6227  iMax = this->mpCopasiModel->getEvents().size();
6228 
6229  for (i = 0; i < iMax; ++i)
6230  {
6231  CEvent* pEvent = this->mpCopasiModel->getEvents()[i];
6232  assert(pEvent != NULL);
6233  const CEvaluationTree* pTree = pEvent->getTriggerExpressionPtr();
6234  // an event has to have a trigger
6235  assert(pTree != NULL);
6236 
6237  if (pTree != NULL)
6238  {
6239  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
6240  }
6241 
6242  // handle the delay
6243  pTree = pEvent->getDelayExpressionPtr();
6244 
6245  if (pTree != NULL)
6246  {
6247  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
6248  }
6249 
6250  // handle all assignments
6251  size_t j, jMax = pEvent->getAssignments().size();
6252 
6253  for (j = 0; j < jMax; ++j)
6254  {
6255  CEventAssignment* pEventAssignment = pEvent->getAssignments()[j];
6256  assert(pEventAssignment != NULL);
6257 
6258  if (pEventAssignment != NULL)
6259  {
6260  pTree = pEventAssignment->getExpressionPtr();
6261  // each event assignment has to have an expression
6262  assert(pTree != NULL);
6263 
6264  if (pTree != NULL)
6265  {
6266  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
6267  }
6268  }
6269  }
6270  }
6271 
6273 
6274  if (mpImportHandler)
6275  {
6276  mpImportHandler->finishItem(hStep);
6277  step = 0;
6278  totalSteps = (unsigned C_INT32) pTmpFunctionDB->loadedFunctions().size();
6279  hStep = mpImportHandler->addItem("Removing unused functions...",
6280  step,
6281  &totalSteps);
6282  }
6283 
6284  // here we could have a dialog asking the user if unused functions should
6285  // be removed.
6286 
6289 
6290  for (; it != end; ++it)
6291  {
6292  CEvaluationTree * pTree = *it;
6293 
6294  if (functionNameSet.find(pTree->getObjectName()) == functionNameSet.end())
6295  {
6296  mUsedFunctions.erase(pTree->getObjectName());
6297  pFunctionDB->loadedFunctions().remove(pTree->getObjectName());
6298  // delete the entry from the copasi2sbmlmap.
6299  std::map<CCopasiObject*, SBase*>::iterator pos = copasi2sbmlmap.find(pTree);
6300  assert(pos != copasi2sbmlmap.end());
6301  copasi2sbmlmap.erase(pos);
6302  }
6303 
6304  ++step;
6305 
6306  if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return false;
6307  }
6308 
6309  if (mpImportHandler)
6310  {
6311  mpImportHandler->finishItem(hStep);
6312  }
6313  }
6314 
6315  return true;
6316 }
6317 
6318 void SBMLImporter::findFunctionCalls(const CEvaluationNode* pNode, std::set<std::string>& functionNameSet)
6319 {
6320  if (pNode)
6321  {
6324 
6325  while (treeIt != NULL)
6326  {
6327  if (CEvaluationNode::type((*treeIt).getType()) == CEvaluationNode::CALL)
6328  {
6329  // unQuote not necessary since getIndex in CCopasiVector takes care of this.
6330  CEvaluationTree* pTree = pFunctionDB->findFunction((*treeIt).getData());
6331 
6332  if (functionNameSet.find(pTree->getObjectName()) == functionNameSet.end())
6333  {
6334  functionNameSet.insert(pTree->getObjectName());
6335  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
6336  }
6337  }
6338 
6339  ++treeIt;
6340  }
6341  }
6342 }
6343 
6344 bool SBMLImporter::isStochasticModel(const Model* pSBMLModel)
6345 {
6346  bool stochastic = false;
6347  unsigned int i;
6348  const UnitDefinition* pUD = pSBMLModel->getUnitDefinition("substance");
6349 
6350  if (pUD)
6351  {
6352  stochastic = (pUD->getNumUnits() == 1 &&
6353  pUD->getUnit(0)->getKind() == UNIT_KIND_ITEM);
6354 
6355  for (i = 0; (stochastic == true) && (i < pSBMLModel->getNumReactions()); ++i)
6356  {
6357  stochastic = !pSBMLModel->getReaction(i)->getReversible();
6358  }
6359  }
6360 
6361  return stochastic;
6362 }
6363 
6364 void SBMLImporter::importSBMLRule(const Rule* sbmlRule, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap, Model* pSBMLModel)
6365 {
6366  // so far we only support assignment rules and rate rules
6367  int type = sbmlRule->getTypeCode();
6368 
6369  if (type == SBML_ASSIGNMENT_RULE)
6370  {
6371  const AssignmentRule* pAssignmentRule = dynamic_cast<const AssignmentRule*>(sbmlRule);
6372 
6373  if (pAssignmentRule && pAssignmentRule->isSetVariable())
6374  {
6375  this->importRule(pAssignmentRule, CModelEntity::ASSIGNMENT, copasi2sbmlmap, pSBMLModel);
6376  }
6377  else
6378  {
6379  fatalError();
6380  }
6381  }
6382  else if (type == SBML_RATE_RULE)
6383  {
6384  const RateRule* pRateRule = dynamic_cast<const RateRule*