COPASI API  4.16.103
CSBMLExporter.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) 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include <cmath>
16 
17 #define USE_LAYOUT 1
18 
19 #include "CSBMLExporter.h"
20 #include "SBMLUtils.h"
21 
22 #include "SBMLImporter.h"
24 #include "sbml/SBMLDocument.h"
25 #include "sbml/Compartment.h"
26 #if LIBSBML_VERSION >= 40100
27 #include "sbml/LocalParameter.h"
28 #endif // LIBSBML_VERSION
29 
30 #if LIBSBML_VERSION >= 50400
31 
32 #include <sbml/packages/layout/extension/LayoutModelPlugin.h>
33 #include <sbml/packages/layout/extension/LayoutExtension.h>
34 #include <sbml/conversion/ConversionProperties.h>
35 
36 #include "layout/CLDefaultStyles.h"
37 #include <sbml/packages/render/extension/RenderExtension.h>
38 #include <sbml/packages/render/extension/RenderListOfLayoutsPlugin.h>
39 #include <sbml/packages/render/sbml/GlobalRenderInformation.h>
40 
41 #define INIT_DEFAULTS(element) \
42  {\
43  element.initDefaults();\
44  }
45 
46 #else
47 
48 #define INIT_DEFAULTS(element) \
49  {\
50  }
51 
52 #endif // LIBSBML VERSION
53 
54 #include "sbml/Model.h"
55 #include "sbml/Species.h"
56 #include "sbml/Parameter.h"
57 #include "sbml/Reaction.h"
58 #include "sbml/KineticLaw.h"
59 #include "sbml/SBMLWriter.h"
60 #include "sbml/SpeciesReference.h"
61 #include "sbml/math/ASTNode.h"
62 #include "sbml/annotation/ModelHistory.h"
63 #include "sbml/annotation/CVTerm.h"
64 #include "sbml/SBMLErrorLog.h"
65 #include "sbml/SBMLError.h"
67 #include "SBMLIncompatibility.h"
68 #include "model/CCompartment.h"
69 #include "model/CModel.h"
70 #include "model/CEvent.h"
71 #include "model/CMetab.h"
72 #include "function/CExpression.h"
74 #include "model/CAnnotation.h"
75 #include "model/CReaction.h"
77 #include "model/CModelValue.h"
78 #include "function/CFunction.h"
79 #include "report/CKeyFactory.h"
80 #include "ConverterASTNode.h"
81 #include "utilities/CCopasiTree.h"
82 #include "model/CChemEqElement.h"
83 #include "utilities/CVersion.h"
84 #include "sbml/Trigger.h"
85 #include "sbml/Event.h"
86 #include "sbml/EventAssignment.h"
87 #include <sbml/xml/XMLInputStream.h>
89 #include "MIRIAM/CRDFUtilities.h"
91 #include "MIRIAM/CReference.h"
93 #include "MIRIAM/CConstants.h"
94 #include "MIRIAM/CCreator.h"
95 #include "MIRIAM/CModified.h"
96 #include "MIRIAM/CRDFPredicate.h"
97 #include "layout/CListOfLayouts.h"
99 #include "utilities/CVersion.h"
101 
102 // helper functions for function definitions
103 
104 std::string hasFunctionDefinitionForURI(SBMLDocument* pSBMLDocument,
105  const std::string& sNamespace,
106  const std::string& elementName,
107  const std::string& definition)
108 {
109  if (pSBMLDocument == NULL || pSBMLDocument->getModel() == NULL) return "";
110 
111  for (unsigned int i = 0; i < pSBMLDocument->getModel()->getNumFunctionDefinitions(); ++i)
112  {
113  FunctionDefinition* current = pSBMLDocument->getModel()->getFunctionDefinition(i);
114 
115  if (current == NULL) continue;
116 
117  if (!current->isSetAnnotation()) continue;
118 
119  const XMLNode* element = current->getAnnotation();
120 
121  if (element == NULL) continue;
122 
123  for (unsigned int i = 0 ; i < element->getNumChildren(); ++i)
124  {
125  const XMLNode& annot = element->getChild(i);
126 
127  if (annot.getURI() == sNamespace &&
128  annot.getName() == elementName &&
129  annot.getAttrValue("definition") == definition)
130  {
131  return current->getId();
132  }
133  }
134  }
135 
136  return "";
137 }
138 
139 std::string createFunctionDefinitonForURI(SBMLDocument* pSBMLDocument,
140  std::map<std::string, const SBase*>& idMap,
141  const char* id,
142  const std::string& sNamespace,
143  const std::string& elementName,
144  const std::string& definition,
145  const std::string& lambda)
146 {
147  if (pSBMLDocument == NULL || pSBMLDocument->getModel() == NULL) return id;
148 
149  std::string newId = CSBMLExporter::createUniqueId(idMap, id, false);
150 
151  FunctionDefinition *def = pSBMLDocument->getModel()->createFunctionDefinition();
152  def -> setId(newId);
153  def -> setMath(SBML_parseFormula(lambda.c_str()));
154 
155  std::stringstream annotation;
156  std::string annotElement = pSBMLDocument->getLevel() == 1 ? "annotations" : "annotation";
157  annotation << "<" << annotElement << "> <" << elementName
158  << " xmlns='" << sNamespace
159  << "' definition='" << definition
160  << "' /> </" << annotElement << ">";
161 
162  def->setAnnotation(annotation.str());
163 
164  return newId;
165 }
166 
167 std::string addRateOfIfItDoesNotExist(SBMLDocument* pSBMLDocument,
168  std::map<std::string, const SBase*>& idMap,
169  const char* id)
170 {
171  std::string newId = hasFunctionDefinitionForURI(pSBMLDocument,
172  "http://sbml.org/annotations/symbols",
173  "symbols",
174  "http://en.wikipedia.org/wiki/Derivative");
175 
176  if (!newId.empty()) return newId;
177 
179  pSBMLDocument,
180  idMap,
181  id,
182  "http://sbml.org/annotations/symbols",
183  "symbols",
184  "http://en.wikipedia.org/wiki/Derivative",
185  "lambda(a,NaN)"
186  );
187  return newId;
188 }
189 
190 std::string getUserDefinedFuctionForName(SBMLDocument* pSBMLDocument,
191  std::map<std::string, const SBase*>& idMap,
192  const char* id)
193 {
194  std::string newId;
195 
196  if (id == std::string("RNORMAL"))
197  {
198  newId = hasFunctionDefinitionForURI(pSBMLDocument,
199  "http://sbml.org/annotations/distribution",
200  "distribution",
201  "http://www.uncertml.org/distributions/normal");
202 
203  if (!newId.empty()) return newId;
204 
206  pSBMLDocument,
207  idMap,
208  id,
209  "http://sbml.org/annotations/distribution",
210  "distribution",
211  "http://www.uncertml.org/distributions/normal",
212  "lambda(m,s,m)"
213  );
214  return newId;
215  }
216  else if (id == std::string("RUNIFORM"))
217  {
218  newId = hasFunctionDefinitionForURI(pSBMLDocument,
219  "http://sbml.org/annotations/distribution",
220  "distribution",
221  "http://www.uncertml.org/distributions/uniform");
222 
223  if (!newId.empty()) return newId;
224 
226  pSBMLDocument,
227  idMap,
228  id,
229  "http://sbml.org/annotations/distribution",
230  "distribution",
231  "http://www.uncertml.org/distributions/uniform",
232  "lambda(a,b,(a+b)/2)"
233  );
234  return newId;
235  }
236  else if (id == std::string("RGAMMA"))
237  {
238  newId = hasFunctionDefinitionForURI(pSBMLDocument,
239  "http://sbml.org/annotations/distribution",
240  "distribution",
241  "http://www.uncertml.org/distributions/gamma");
242 
243  if (!newId.empty()) return newId;
244 
246  pSBMLDocument,
247  idMap,
248  id,
249  "http://sbml.org/annotations/distribution",
250  "distribution",
251  "http://www.uncertml.org/distributions/gamma",
252  "lambda(a,b,a*b)"
253  );
254  return newId;
255  }
256  else if (id == std::string("RPOISSON"))
257  {
258  newId = hasFunctionDefinitionForURI(pSBMLDocument,
259  "http://sbml.org/annotations/distribution",
260  "distribution",
261  "http://www.uncertml.org/distributions/poisson");
262 
263  if (!newId.empty()) return newId;
264 
266  pSBMLDocument,
267  idMap,
268  id,
269  "http://sbml.org/annotations/distribution",
270  "distribution",
271  "http://www.uncertml.org/distributions/poisson",
272  "lambda(mu,mu)"
273  );
274  return newId;
275  }
276  else if (id == std::string("MAX"))
277  {
278  newId = hasFunctionDefinitionForURI(pSBMLDocument,
279  "http://sbml.org/annotations/function",
280  "function",
281  "http://sbml.org/annotations/function/max");
282 
283  if (!newId.empty()) return newId;
284 
286  pSBMLDocument,
287  idMap,
288  id,
289  "http://sbml.org/annotations/function",
290  "function",
291  "http://sbml.org/annotations/function/max",
292  "lambda(a,b,piecewise(a,geq(a,b),b))"
293  );
294  return newId;
295  }
296  else if (id == std::string("MIN"))
297  {
298  newId = hasFunctionDefinitionForURI(pSBMLDocument,
299  "http://sbml.org/annotations/function",
300  "function",
301  "http://sbml.org/annotations/function/min");
302 
303  if (!newId.empty()) return newId;
304 
306  pSBMLDocument,
307  idMap,
308  id,
309  "http://sbml.org/annotations/function",
310  "function",
311  "http://sbml.org/annotations/function/min",
312  "lambda(a,b,piecewise(a,leq(a,b),b))"
313  );
314  return newId;
315  }
316  else if (id == std::string("rateOf"))
317  {
318  return addRateOfIfItDoesNotExist(pSBMLDocument, idMap, id);
319  }
320 
321  return id;
322 }
323 
324 #ifdef USE_SBMLUNIT
326 #endif // USE_SBMLUNIT
327 
328 CSBMLExporter::CSBMLExporter(): mpSBMLDocument(NULL), mSBMLLevel(2), mSBMLVersion(1), mIncompleteExport(false), mVariableVolumes(false), mpAvogadro(NULL), mAvogadroCreated(false), mMIRIAMWarning(false), mDocumentDisowned(false), mExportCOPASIMIRIAM(false), mExportedFunctions(2, 1)
329 {}
330 
332 {
333  if (this->mDocumentDisowned == false)
334  {
336  }
337 };
338 
339 /**
340  * Creates the units for the SBML model.
341  */
343 {
344  createLengthUnit(dataModel);
345  createAreaUnit(dataModel);
346  createVolumeUnit(dataModel);
347  createTimeUnit(dataModel);
348  createSubstanceUnit(dataModel);
349 }
350 
351 /**
352  * Creates the time unit for the SBML model.
353  */
355 {
356  if (dataModel.getModel() == NULL || this->mpSBMLDocument == NULL || this->mpSBMLDocument->getModel() == NULL) return;
357 
358  UnitDefinition uDef(this->mSBMLLevel, this->mSBMLVersion);
359  uDef.setName("time");
360  uDef.setId("time");
361  Unit unit(this->mSBMLLevel, this->mSBMLVersion);
362  INIT_DEFAULTS(unit);
363 
364  switch (dataModel.getModel()->getTimeUnitEnum())
365  {
366  case CModel::d:
367  unit.setKind(UNIT_KIND_SECOND);
368  unit.setExponent(1);
369  unit.setScale(0);
370  unit.setMultiplier(86400);
371  break;
372 
373  case CModel::h:
374  unit.setKind(UNIT_KIND_SECOND);
375  unit.setExponent(1);
376  unit.setScale(0);
377  unit.setMultiplier(3600);
378  break;
379 
380  case CModel::min:
381  unit.setKind(UNIT_KIND_SECOND);
382  unit.setExponent(1);
383  unit.setScale(0);
384  unit.setMultiplier(60);
385  break;
386 
387  case CModel::s:
388  unit.setKind(UNIT_KIND_SECOND);
389  unit.setExponent(1);
390  unit.setScale(0);
391  unit.setMultiplier(1);
392  break;
393 
394  case CModel::ms:
395  unit.setKind(UNIT_KIND_SECOND);
396  unit.setExponent(1);
397  unit.setScale(-3);
398  unit.setMultiplier(1);
399  break;
400 
401  case CModel::micros:
402  unit.setKind(UNIT_KIND_SECOND);
403  unit.setExponent(1);
404  unit.setScale(-6);
405  unit.setMultiplier(1);
406  break;
407 
408  case CModel::ns:
409  unit.setKind(UNIT_KIND_SECOND);
410  unit.setExponent(1);
411  unit.setScale(-9);
412  unit.setMultiplier(1);
413  break;
414 
415  case CModel::ps:
416  unit.setKind(UNIT_KIND_SECOND);
417  unit.setExponent(1);
418  unit.setScale(-12);
419  unit.setMultiplier(1);
420  break;
421 
422  case CModel::fs:
423  unit.setKind(UNIT_KIND_SECOND);
424  unit.setExponent(1);
425  unit.setScale(-15);
426  unit.setMultiplier(1);
427  break;
428 
430  unit.setKind(UNIT_KIND_DIMENSIONLESS);
431  unit.setExponent(1);
432  unit.setScale(0);
433  unit.setMultiplier(1);
434  break;
435 
436  default:
437  CCopasiMessage(CCopasiMessage::EXCEPTION, "SBMLExporter Error: Unknown copasi time unit.");
438  break;
439  }
440 
441  uDef.addUnit(&unit);
442  Model* pSBMLModel = this->mpSBMLDocument->getModel();
443  UnitDefinition* pUdef = pSBMLModel->getUnitDefinition("time");
444 
445  if (pUdef != NULL)
446  {
447  // check if it is the same unit as the existing one if there is one
448  // if yes, return, else replace the existing one
450  {
451  (*pUdef) = uDef;
452  }
453  }
454  else
455  {
456  // only add it if it is not the default unit definition anyway
457  // from SBML Level 3 on, there are no default units
458  // so we have to write it
459  if (this->mSBMLLevel > 2 || unit.getKind() != UNIT_KIND_SECOND || unit.getScale() != 0 || unit.getExponent() != 1 || unit.getMultiplier() != 1.0)
460  {
461  // set the unit definition
462  pSBMLModel->addUnitDefinition(&uDef);
463  }
464  }
465 
466 // if we write an SBML L3 document, we have to explicitely set the units on the model
467 #if LIBSBML_VERSION >= 40100
468 
469  if (this->mSBMLLevel > 2)
470  {
471  pSBMLModel->setTimeUnits(uDef.getId());
472  }
473 
474 #endif // LIBSBML_VERSION
475 }
476 
477 /**
478  * Creates the volume unit for the SBML model.
479  */
481 {
482  if (dataModel.getModel() == NULL || this->mpSBMLDocument == NULL || this->mpSBMLDocument->getModel() == NULL) return;
483 
484  UnitDefinition uDef(this->mSBMLLevel, this->mSBMLVersion);
485  uDef.setName("volume");
486  uDef.setId("volume");
487  Unit unit(this->mSBMLLevel, this->mSBMLVersion);
488  INIT_DEFAULTS(unit);
489 
490  switch (dataModel.getModel()->getVolumeUnitEnum())
491  {
492  case CModel::l:
493  unit.setKind(UNIT_KIND_LITRE);
494  unit.setExponent(1);
495  unit.setScale(0);
496  break;
497 
498  case CModel::ml:
499  unit.setKind(UNIT_KIND_LITRE);
500  unit.setExponent(1);
501  unit.setScale(-3);
502  break;
503 
504  case CModel::microl:
505  unit.setKind(UNIT_KIND_LITRE);
506  unit.setExponent(1);
507  unit.setScale(-6);
508  break;
509 
510  case CModel::nl:
511  unit.setKind(UNIT_KIND_LITRE);
512  unit.setExponent(1);
513  unit.setScale(-9);
514  break;
515 
516  case CModel::pl:
517  unit.setKind(UNIT_KIND_LITRE);
518  unit.setExponent(1);
519  unit.setScale(-12);
520  break;
521 
522  case CModel::fl:
523  unit.setKind(UNIT_KIND_LITRE);
524  unit.setExponent(1);
525  unit.setScale(-15);
526  break;
527 
528  case CModel::m3:
529  unit.setKind(UNIT_KIND_METRE);
530  unit.setExponent(3);
531  unit.setScale(0);
532  break;
533 
535  unit.setKind(UNIT_KIND_DIMENSIONLESS);
536  unit.setExponent(1);
537  unit.setScale(0);
538  break;
539 
540  default:
541  CCopasiMessage(CCopasiMessage::EXCEPTION, "SBMLExporter Error: Unknown copasi volume unit.");
542  break;
543  }
544 
545  unit.setMultiplier(1.0);
546  uDef.addUnit(&unit);
547  Model* pSBMLModel = this->mpSBMLDocument->getModel();
548  UnitDefinition* pUdef = pSBMLModel->getUnitDefinition("volume");
549 
550  if (pUdef != NULL)
551  {
552  // check if it is the same unit as the existing one if there is one
553  // if yes, return, else replace the existing one
555  {
556  (*pUdef) = uDef;
557  }
558  }
559  else
560  {
561  // only add it if it is not the default unit definition anyway
562  if (this->mSBMLLevel > 2 || unit.getKind() != UNIT_KIND_LITRE || unit.getScale() != 0 || unit.getExponent() != 1 || unit.getMultiplier() != 1.0)
563  {
564  // set the unit definition
565  pSBMLModel->addUnitDefinition(&uDef);
566  }
567  }
568 
569 // if we write an SBML L3 document, we have to explicitely set the units on the model
570 #if LIBSBML_VERSION >= 40100
571 
572  if (this->mSBMLLevel > 2)
573  {
574  pSBMLModel->setVolumeUnits(uDef.getId());
575  }
576 
577 #endif // LIBSBML_VERSION
578 }
579 
580 /**
581  * Creates the substance unit for the SBML model.
582  */
584 {
585  if (dataModel.getModel() == NULL || this->mpSBMLDocument == NULL || this->mpSBMLDocument->getModel() == NULL) return;
586 
587  UnitDefinition uDef(this->mSBMLLevel, this->mSBMLVersion);
588  uDef.setName("substance");
589  uDef.setId("substance");
590  Unit unit(this->mSBMLLevel, this->mSBMLVersion);
591  INIT_DEFAULTS(unit);
592 
593  switch (dataModel.getModel()->getQuantityUnitEnum())
594  {
595  case CModel::Mol:
596  unit.setKind(UNIT_KIND_MOLE);
597  unit.setExponent(1);
598  unit.setScale(0);
599  break;
600 
601  case CModel::mMol:
602  unit.setKind(UNIT_KIND_MOLE);
603  unit.setExponent(1);
604  unit.setScale(-3);
605  break;
606 
607  case CModel::microMol:
608  unit.setKind(UNIT_KIND_MOLE);
609  unit.setExponent(1);
610  unit.setScale(-6);
611  break;
612 
613  case CModel::nMol:
614  unit.setKind(UNIT_KIND_MOLE);
615  unit.setExponent(1);
616  unit.setScale(-9);
617  break;
618 
619  case CModel::pMol:
620  unit.setKind(UNIT_KIND_MOLE);
621  unit.setExponent(1);
622  unit.setScale(-12);
623  break;
624 
625  case CModel::fMol:
626  unit.setKind(UNIT_KIND_MOLE);
627  unit.setExponent(1);
628  unit.setScale(-15);
629  break;
630 
631  case CModel::number:
632  unit.setKind(UNIT_KIND_ITEM);
633  unit.setExponent(1);
634  unit.setScale(0);
635  break;
636 
638  unit.setKind(UNIT_KIND_DIMENSIONLESS);
639  unit.setExponent(1);
640  unit.setScale(0);
641  break;
642 
643  default:
644  CCopasiMessage(CCopasiMessage::EXCEPTION, "SBMLExporter Error: Unknown copasi quantity unit.");
645  break;
646  }
647 
648  unit.setMultiplier(1);
649  uDef.addUnit(&unit);
650  Model* pSBMLModel = this->mpSBMLDocument->getModel();
651  UnitDefinition* pUdef = pSBMLModel->getUnitDefinition("substance");
652 
653  if (pUdef != NULL)
654  {
655  // check if it is the same unit as the existing one if there is one
656  // if yes, return, else replace the existing one
658  {
659  (*pUdef) = uDef;
660  }
661  }
662  else
663  {
664  // only add it if it is not the default unit definition anyway
665  if (this->mSBMLLevel > 2 || unit.getKind() != UNIT_KIND_MOLE || unit.getScale() != 0 || unit.getExponent() != 1 || unit.getMultiplier() != 1.0)
666  {
667  // set the unit definition
668  pSBMLModel->addUnitDefinition(&uDef);
669  }
670  }
671 
672 // if we write an SBML L3 document, we have to explicitely set the units on the model
673 #if LIBSBML_VERSION >= 40100
674 
675  if (this->mSBMLLevel > 2)
676  {
677  pSBMLModel->setSubstanceUnits(uDef.getId());
678  // here we also set the extent unit to the same unit as the substance unit
679  // because COPASI does not know about different extent units
680  pSBMLModel->setExtentUnits(uDef.getId());
681  }
682 
683 #endif // LIBSBML_VERSION
684 }
685 
686 /**
687  * Creates the length unit for the SBML model.
688  */
690 {
691  if (dataModel.getModel() == NULL || this->mpSBMLDocument == NULL || this->mpSBMLDocument->getModel() == NULL) return;
692 
693  UnitDefinition uDef(this->mSBMLLevel, this->mSBMLVersion);
694  uDef.setName("length");
695  uDef.setId("length");
696  Unit unit(this->mSBMLLevel, this->mSBMLVersion);
697  INIT_DEFAULTS(unit);
698 
699  switch (dataModel.getModel()->getLengthUnitEnum())
700  {
701  case CModel::m:
702  unit.setKind(UNIT_KIND_METRE);
703  unit.setExponent(1);
704  unit.setScale(0);
705  break;
706 
707  case CModel::dm:
708  unit.setKind(UNIT_KIND_METRE);
709  unit.setExponent(1);
710  unit.setScale(-1);
711  break;
712 
713  case CModel::cm:
714  unit.setKind(UNIT_KIND_METRE);
715  unit.setExponent(1);
716  unit.setScale(-2);
717  break;
718 
719  case CModel::mm:
720  unit.setKind(UNIT_KIND_METRE);
721  unit.setExponent(1);
722  unit.setScale(-3);
723  break;
724 
725  case CModel::microm:
726  unit.setKind(UNIT_KIND_METRE);
727  unit.setExponent(1);
728  unit.setScale(-6);
729  break;
730 
731  case CModel::nm:
732  unit.setKind(UNIT_KIND_METRE);
733  unit.setExponent(1);
734  unit.setScale(-9);
735  break;
736 
737  case CModel::pm:
738  unit.setKind(UNIT_KIND_METRE);
739  unit.setExponent(1);
740  unit.setScale(-12);
741  break;
742 
743  case CModel::fm:
744  unit.setKind(UNIT_KIND_METRE);
745  unit.setExponent(1);
746  unit.setScale(-15);
747  break;
748 
750  unit.setKind(UNIT_KIND_DIMENSIONLESS);
751  unit.setExponent(1);
752  unit.setScale(0);
753  break;
754 
755  default:
756  CCopasiMessage(CCopasiMessage::EXCEPTION, "SBMLExporter Error: Unknown copasi length unit.");
757  break;
758  }
759 
760  unit.setMultiplier(1.0);
761  uDef.addUnit(&unit);
762  Model* pSBMLModel = this->mpSBMLDocument->getModel();
763  UnitDefinition* pUdef = pSBMLModel->getUnitDefinition("length");
764 
765  if (pUdef != NULL)
766  {
767  // check if it is the same unit as the existing one if there is one
768  // if yes, return, else replace the existing one
770  {
771  (*pUdef) = uDef;
772  }
773  }
774  else
775  {
776  // only add it if it is not the default unit definition anyway
777  if (this->mSBMLLevel > 2 || unit.getKind() != UNIT_KIND_METRE || unit.getScale() != 0 || unit.getExponent() != 1 || unit.getMultiplier() != 1.0)
778  {
779  // set the unit definition
780  pSBMLModel->addUnitDefinition(&uDef);
781  }
782  }
783 
784 // if we write an SBML L3 document, we have to explicitely set the units on the model
785 #if LIBSBML_VERSION >= 40100
786 
787  if (this->mSBMLLevel > 2)
788  {
789  pSBMLModel->setLengthUnits(uDef.getId());
790  }
791 
792 #endif // LIBSBML_VERSION
793 }
794 
795 /**
796  * Creates the area unit for the SBML model.
797  */
799 {
800  if (dataModel.getModel() == NULL || this->mpSBMLDocument == NULL || this->mpSBMLDocument->getModel() == NULL) return;
801 
802  UnitDefinition uDef(this->mSBMLLevel, this->mSBMLVersion);
803  uDef.setName("area");
804  uDef.setId("area");
805  Unit unit(this->mSBMLLevel, this->mSBMLVersion);
806  INIT_DEFAULTS(unit);
807 
808  switch (dataModel.getModel()->getAreaUnitEnum())
809  {
810  case CModel::m:
811  unit.setKind(UNIT_KIND_METRE);
812  unit.setExponent(2);
813  unit.setScale(0);
814  break;
815 
816  case CModel::dm:
817  unit.setKind(UNIT_KIND_METRE);
818  unit.setExponent(2);
819  unit.setScale(-1);
820  break;
821 
822  case CModel::cm:
823  unit.setKind(UNIT_KIND_METRE);
824  unit.setExponent(2);
825  unit.setScale(-2);
826  break;
827 
828  case CModel::mm:
829  unit.setKind(UNIT_KIND_METRE);
830  unit.setExponent(2);
831  unit.setScale(-3);
832  break;
833 
834  case CModel::microm:
835  unit.setKind(UNIT_KIND_METRE);
836  unit.setExponent(2);
837  unit.setScale(-6);
838  break;
839 
840  case CModel::nm:
841  unit.setKind(UNIT_KIND_METRE);
842  unit.setExponent(2);
843  unit.setScale(-9);
844  break;
845 
846  case CModel::pm:
847  unit.setKind(UNIT_KIND_METRE);
848  unit.setExponent(2);
849  unit.setScale(-12);
850  break;
851 
852  case CModel::fm:
853  unit.setKind(UNIT_KIND_METRE);
854  unit.setExponent(2);
855  unit.setScale(-15);
856  break;
857 
859  unit.setKind(UNIT_KIND_DIMENSIONLESS);
860  unit.setExponent(1);
861  unit.setScale(0);
862  break;
863 
864  default:
865  CCopasiMessage(CCopasiMessage::EXCEPTION, "SBMLExporter Error: Unknown copasi area unit.");
866  break;
867  }
868 
869  unit.setMultiplier(1.0);
870  uDef.addUnit(&unit);
871  Model* pSBMLModel = this->mpSBMLDocument->getModel();
872  UnitDefinition* pUdef = pSBMLModel->getUnitDefinition("area");
873 
874  if (pUdef != NULL)
875  {
876  // check if it is the same unit as the existing one if there is one
877  // if yes, return, else replace the existing one
879  {
880  (*pUdef) = uDef;
881  }
882  }
883  else
884  {
885  // only add it if it is not the default unit definition anyway
886  if (this->mSBMLLevel > 2 ||
887  unit.getKind() != UNIT_KIND_METRE ||
888  unit.getScale() != 0 ||
889  unit.getExponent() != 2 ||
890  unit.getMultiplier() != 1.0)
891  {
892  // set the unit definition
893  pSBMLModel->addUnitDefinition(&uDef);
894  }
895  }
896 
897 // if we write an SBML L3 document, we have to explicitely set the units on the model
898 #if LIBSBML_VERSION >= 40100
899 
900  if (this->mSBMLLevel > 2)
901  {
902  pSBMLModel->setAreaUnits(uDef.getId());
903  }
904 
905 #endif // LIBSBML_VERSION
906 }
907 
908 /**
909  * Creates the compartments for the model.
910  */
912 {
913  // make sure the SBML Document already exists and that it has a Model set
914  if (dataModel.getModel() == NULL || this->mpSBMLDocument == NULL || this->mpSBMLDocument->getModel() == NULL) return;
915 
917 
918  while (it != endit)
919  {
920  createCompartment(**it);
921  ++it;
922  }
923 }
924 
925 /**
926  * Creates the compartment for the given COPASI compartment.
927  */
929 {
930  Compartment* pSBMLCompartment = NULL;
931  std::string sbmlId = compartment.getSBMLId();
932 
933  if (!sbmlId.empty())
934  {
935  pSBMLCompartment = this->mpSBMLDocument->getModel()->getCompartment(sbmlId);
936 
937  if (pSBMLCompartment == NULL)
938  {
939  pSBMLCompartment = this->mpSBMLDocument->getModel()->createCompartment();
940  this->mCOPASI2SBMLMap[&compartment] = pSBMLCompartment;
941  pSBMLCompartment->setId(sbmlId);
942  }
943  }
944  else
945  {
946  pSBMLCompartment = this->mpSBMLDocument->getModel()->createCompartment();
947  this->mCOPASI2SBMLMap[&compartment] = pSBMLCompartment;
948  sbmlId = CSBMLExporter::createUniqueId(this->mIdMap, compartment.getObjectName(), false);
949  compartment.setSBMLId(sbmlId);
950  pSBMLCompartment->setId(sbmlId);
951  }
952 
953  INIT_DEFAULTS((*pSBMLCompartment));
954 
955  this->mIdMap.insert(std::pair<const std::string, const SBase*>(sbmlId, pSBMLCompartment));
956  this->mHandledSBMLObjects.insert(pSBMLCompartment);
957 
958  // don't call setName on level 1 objects because this will also
959  // change the id
960  if (this->mpSBMLDocument->getLevel() > 1)
961  {
962  pSBMLCompartment->setName(compartment.getObjectName().c_str());
963  }
964 
965  pSBMLCompartment->setSpatialDimensions((unsigned int)compartment.getDimensionality());
966  double value = compartment.getInitialValue();
967 
968  // if the value is NaN, unset the initial volume
969  if (!isnan(value))
970  {
971  pSBMLCompartment->setVolume(value);
972  }
973  else
974  {
975  pSBMLCompartment->unsetVolume();
976  }
977 
978  // fill assignment set
979  // a compartment can either have an assignment or an initial assignment but not
980  // both
981  CModelEntity::Status status = compartment.getStatus();
982 
983  if (status == CModelEntity::ASSIGNMENT)
984  {
985  if (compartment.getDimensionality() != 0)
986  {
987  this->mAssignmentVector.push_back(&compartment);
988  pSBMLCompartment->setConstant(false);
989  removeInitialAssignment(pSBMLCompartment->getId());
990  }
991  else
992  {
993  fatalError();
994  }
995  }
996  else if (status == CModelEntity::ODE)
997  {
998  if (compartment.getDimensionality() != 0)
999  {
1000  this->mODEVector.push_back(&compartment);
1001  pSBMLCompartment->setConstant(false);
1002 
1003  if (compartment.getInitialExpression() != "")
1004  {
1005  this->mInitialAssignmentVector.push_back(&compartment);
1006  }
1007  else
1008  {
1009  removeInitialAssignment(pSBMLCompartment->getId());
1010  }
1011  }
1012  else
1013  {
1014  fatalError();
1015  }
1016  }
1017  else
1018  {
1019  // do this explicitly since we might handle an existing object from an
1020  // earlier import that had it attribute set already
1021  if (this->mSBMLLevel != 1)
1022  {
1023  pSBMLCompartment->setConstant(true);
1024  }
1025  else
1026  {
1027  // Level 1 does not know the constant flag and libsbmnl does not drop
1028  // it automatically
1029  pSBMLCompartment->setConstant(false);
1030  }
1031 
1032  removeRule(pSBMLCompartment->getId());
1033 
1034  // fill initial assignment set
1035  if (compartment.getInitialExpression() != "")
1036  {
1037  if (compartment.getDimensionality() != 0)
1038  {
1039  this->mInitialAssignmentVector.push_back(&compartment);
1040  }
1041  else
1042  {
1043  fatalError();
1044  }
1045  }
1046  else
1047  {
1048  removeInitialAssignment(pSBMLCompartment->getId());
1049  }
1050  }
1051 
1052  if (pSBMLCompartment != NULL)
1053  {
1054  CSBMLExporter::setSBMLNotes(pSBMLCompartment, &compartment);
1055  }
1056 
1057  if (pSBMLCompartment != NULL && mSBMLLevel == 3)
1058  {
1059  pSBMLCompartment->setUnits("volume");
1060  }
1061 
1062  CSBMLExporter::updateMIRIAMAnnotation(&compartment, pSBMLCompartment, this->mMetaIdMap);
1063 }
1064 
1065 /**
1066  * Creates the compartments for the model.
1067  */
1069 {
1070  // make sure the SBML Document already exists and that it has a Model set
1071  if (dataModel.getModel() == NULL || this->mpSBMLDocument == NULL || this->mpSBMLDocument->getModel() == NULL) return;
1072 
1073  if (this->mSBMLLevel > 2 || (this->mSBMLLevel == 2 && this->mSBMLVersion >= 3))
1074  {
1076  }
1077 
1078  CCopasiVector<CMetab>::const_iterator it = dataModel.getModel()->getMetabolites().begin(), endit = dataModel.getModel()->getMetabolites().end();
1079  this->mSpatialSizeUnitsSpecies.clear();
1080 
1081  while (it != endit)
1082  {
1083  createMetabolite(**it);
1084  ++it;
1085  }
1086 
1087  if (!this->mSpatialSizeUnitsSpecies.empty())
1088  {
1089  std::ostringstream os;
1090  std::set<std::string>::const_iterator sit = this->mSpatialSizeUnitsSpecies.begin(), sendit = this->mSpatialSizeUnitsSpecies.end();
1091 
1092  while (sit != sendit)
1093  {
1094  os << *sit << ", ";
1095  ++sit;
1096  }
1097 
1098  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 84, os.str().substr(0, os.str().size() - 2).c_str());
1099  }
1100 }
1101 
1102 /**
1103  * Creates the species for the given COPASI metabolite.
1104  */
1106 {
1107  Species* pSBMLSpecies = NULL;
1108  std::string sbmlId = metab.getSBMLId();
1109 
1110  if (!sbmlId.empty())
1111  {
1112  pSBMLSpecies = this->mpSBMLDocument->getModel()->getSpecies(sbmlId);
1113 
1114  if (pSBMLSpecies == NULL)
1115  {
1116  pSBMLSpecies = this->mpSBMLDocument->getModel()->createSpecies();
1117  this->mCOPASI2SBMLMap[&metab] = pSBMLSpecies;
1118  pSBMLSpecies->setId(sbmlId);
1119  }
1120  else
1121  {
1122 #if LIBSBML_VERSION >= 40100
1123 
1124  if (this->mSBMLLevel > 2)
1125  {
1126  pSBMLSpecies->unsetConversionFactor();
1127  // we have to remove the conversionFactor because on import we multiplied the stoichiometries
1128  // with this factor
1129  }
1130 
1131 #endif // LIBSBML_VERSION
1132 
1133  // clear the spatialSizeUnits attribute if there is any
1134  if (this->mSBMLLevel > 2 || (this->mSBMLLevel == 2 && this->mSBMLVersion >= 3))
1135  {
1136  if (pSBMLSpecies->isSetSpatialSizeUnits())
1137  {
1138  pSBMLSpecies->unsetSpatialSizeUnits();
1139  this->mSpatialSizeUnitsSpecies.insert(pSBMLSpecies->getId());
1140  }
1141  }
1142  }
1143  }
1144  else
1145  {
1146  pSBMLSpecies = this->mpSBMLDocument->getModel()->createSpecies();
1147  this->mCOPASI2SBMLMap[&metab] = pSBMLSpecies;
1148  sbmlId = CSBMLExporter::createUniqueId(this->mIdMap, metab.getObjectDisplayName(), false);
1149  metab.setSBMLId(sbmlId);
1150  pSBMLSpecies->setId(sbmlId);
1151  }
1152 
1153  INIT_DEFAULTS((*pSBMLSpecies));
1154 
1155  this->mIdMap.insert(std::pair<const std::string, const SBase*>(sbmlId, pSBMLSpecies));
1156  this->mHandledSBMLObjects.insert(pSBMLSpecies);
1157 
1158  // don't call setName on level 1 objects because this will also
1159  // change the id
1160  if (this->mpSBMLDocument->getLevel() > 1)
1161  {
1162  pSBMLSpecies->setName(metab.getObjectName().c_str());
1163  }
1164 
1165  //const Compartment* pSBMLCompartment = dynamic_cast<const Compartment*>(this->mCOPASI2SBMLMap[const_cast<CCompartment*>(metab.getCompartment())]);
1166  const Compartment* pSBMLCompartment = this->mpSBMLDocument->getModel()->getCompartment(metab.getCompartment()->getSBMLId());
1167  assert(pSBMLCompartment);
1168  pSBMLSpecies->setCompartment(pSBMLCompartment->getId());
1169 
1170  if (this->mVariableVolumes == true)
1171  {
1172  pSBMLSpecies->setHasOnlySubstanceUnits(true);
1173  }
1174 
1175  double value = metab.getInitialConcentration();
1176 
1177  // if the value is NaN, unset the initial amount
1178  if (!isnan(value))
1179  {
1180  // if we set the concentration on a species that had the amount
1181  // set, the meaning of the model is different when the user changed
1182  // to size of the corresponding compartment
1183  // So if the amount had been set, we try to keep this.
1184  // we also have to set the initial amount if the model has variable
1185  // volumes since those models export all species with the
1186  // hasOnlySubstanceUnits flag set to true
1187  //
1188  // libsbml 4 does not set the initial concentration on a species any more
1189  // so we also have to set the initial amount
1190  if (pSBMLSpecies->isSetInitialAmount() || this->mVariableVolumes == true || pSBMLSpecies->getLevel() == 1)
1191  {
1192  pSBMLSpecies->setInitialAmount(value * metab.getCompartment()->getInitialValue());
1193  }
1194  else
1195  {
1196  pSBMLSpecies->setInitialConcentration(value);
1197  }
1198  }
1199  else
1200  {
1201  pSBMLSpecies->unsetInitialConcentration();
1202  pSBMLSpecies->unsetInitialAmount();
1203  }
1204 
1205  CModelEntity::Status status = metab.getStatus();
1206 
1207  // a species can either have an assignment or an initial assignment but not
1208  // both
1209  if (status == CModelEntity::ASSIGNMENT)
1210  {
1211  this->mAssignmentVector.push_back(&metab);
1212  pSBMLSpecies->setConstant(false);
1213  pSBMLSpecies->setBoundaryCondition(true);
1214  removeInitialAssignment(pSBMLSpecies->getId());
1215  }
1216  else if (status == CModelEntity::ODE)
1217  {
1218  this->mODEVector.push_back(&metab);
1219  pSBMLSpecies->setConstant(false);
1220  pSBMLSpecies->setBoundaryCondition(true);
1221 
1222  if (metab.getInitialExpression() != "")
1223  {
1224  this->mInitialAssignmentVector.push_back(&metab);
1225  }
1226  else
1227  {
1228  removeInitialAssignment(pSBMLSpecies->getId());
1229  }
1230  }
1231  else if (status == CModelEntity::FIXED)
1232  {
1233  // do this explicitly since we might handle an existing object from an
1234  // earlier import that had it attribute set already
1235  if (this->mSBMLLevel != 1)
1236  {
1237  pSBMLSpecies->setConstant(true);
1238  }
1239  else
1240  {
1241  // Level 1 doesn't have the constant attribute and libSBML 3.2.0 does
1242  // not drop it automatically
1243  pSBMLSpecies->setConstant(false);
1244  }
1245 
1246  pSBMLSpecies->setBoundaryCondition(true);
1247  removeRule(pSBMLSpecies->getId());
1248 
1249  // fill initial assignment set
1250  if (metab.getInitialExpression() != "")
1251  {
1252  this->mInitialAssignmentVector.push_back(&metab);
1253  }
1254  else
1255  {
1256  removeInitialAssignment(pSBMLSpecies->getId());
1257  }
1258  }
1259  else if (status == CModelEntity::REACTIONS)
1260  {
1261  pSBMLSpecies->setConstant(false);
1262  pSBMLSpecies->setBoundaryCondition(false);
1263 
1264  if (metab.getInitialExpression() != "")
1265  {
1266  this->mInitialAssignmentVector.push_back(&metab);
1267  }
1268  }
1269 
1270  if (pSBMLSpecies != NULL)
1271  {
1272  CSBMLExporter::setSBMLNotes(pSBMLSpecies, &metab);
1273  }
1274 
1275  if (pSBMLSpecies != NULL && mSBMLLevel == 3)
1276  {
1277  // for l3 the unit needs to be set explicitly to substance,
1278  pSBMLSpecies->setSubstanceUnits("substance");
1279  }
1280 
1281  CSBMLExporter::updateMIRIAMAnnotation(&metab, pSBMLSpecies, this->mMetaIdMap);
1282 }
1283 
1284 /**
1285  * Creates the parameters for the model.
1286  */
1288 {
1289  // make sure the SBML Document already exists and that it has a Model set
1290  if (dataModel.getModel() == NULL || this->mpSBMLDocument == NULL || this->mpSBMLDocument->getModel() == NULL) return;
1291 
1293 
1294  while (it != endit)
1295  {
1296  createParameter(**it);
1297  ++it;
1298  }
1299 }
1300 
1301 /**
1302  * Creates the parameter for the given COPASI parameter.
1303  */
1305 {
1306  Parameter* pParameter = NULL;
1307  std::string sbmlId = modelValue.getSBMLId();
1308 
1309  if (!sbmlId.empty())
1310  {
1311  pParameter = this->mpSBMLDocument->getModel()->getParameter(sbmlId);
1312 
1313  if (pParameter == NULL)
1314  {
1315  pParameter = this->mpSBMLDocument->getModel()->createParameter();
1316  this->mCOPASI2SBMLMap[&modelValue] = pParameter;
1317  pParameter->setId(sbmlId);
1318  }
1319  }
1320  else
1321  {
1322  pParameter = this->mpSBMLDocument->getModel()->createParameter();
1323  this->mCOPASI2SBMLMap[&modelValue] = pParameter;
1324  sbmlId = CSBMLExporter::createUniqueId(this->mIdMap, modelValue.getObjectName(), false);
1325  modelValue.setSBMLId(sbmlId);
1326  pParameter->setId(sbmlId);
1327  }
1328 
1329  INIT_DEFAULTS((*pParameter));
1330 
1331  this->mIdMap.insert(std::pair<const std::string, const SBase*>(sbmlId, pParameter));
1332  this->mHandledSBMLObjects.insert(pParameter);
1333 
1334  // don't call setName on level 1 objects because this will also
1335  // change the id
1336  if (this->mpSBMLDocument->getLevel() > 1)
1337  {
1338  pParameter->setName(modelValue.getObjectName());
1339  }
1340 
1341  double value = modelValue.getInitialValue();
1342 
1343  // if the value is NaN, unset the parameters value
1344  if (!isnan(value))
1345  {
1346  pParameter->setValue(value);
1347  }
1348  else
1349  {
1350  pParameter->unsetValue();
1351  }
1352 
1353  // fill assignment set
1354  // a parameter can either have an assignment or an initial assignment but not
1355  // both
1356  CModelEntity::Status status = modelValue.getStatus();
1357 
1358  if (status == CModelEntity::ASSIGNMENT)
1359  {
1360  this->mAssignmentVector.push_back(&modelValue);
1361  pParameter->setConstant(false);
1362  removeInitialAssignment(pParameter->getId());
1363  }
1364  else if (status == CModelEntity::ODE)
1365  {
1366  this->mODEVector.push_back(&modelValue);
1367  pParameter->setConstant(false);
1368 
1369  if (modelValue.getInitialExpression() != "")
1370  {
1371  this->mInitialAssignmentVector.push_back(&modelValue);
1372  }
1373  else
1374  {
1375  removeInitialAssignment(pParameter->getId());
1376  }
1377  }
1378  else
1379  {
1380  // do this explicitly since we might handle an existing object from an
1381  // earlier import that had it attribute set already
1382  if (this->mSBMLLevel != 1)
1383  {
1384  pParameter->setConstant(true);
1385  }
1386  else
1387  {
1388  // Level 1 does not have a constant flag
1389  // and libsbml does not drop it automatically
1390  pParameter->setConstant(false);
1391  }
1392 
1393  // fill initial assignment set
1394  removeRule(pParameter->getId());
1395 
1396  if (modelValue.getInitialExpression() != "")
1397  {
1398  this->mInitialAssignmentVector.push_back(&modelValue);
1399  }
1400  else
1401  {
1402  removeInitialAssignment(pParameter->getId());
1403  }
1404  }
1405 
1406  if (pParameter != NULL)
1407  {
1408  CSBMLExporter::setSBMLNotes(pParameter, &modelValue);
1409  }
1410 
1411  CSBMLExporter::updateMIRIAMAnnotation(&modelValue, pParameter, this->mMetaIdMap);
1412 }
1413 
1414 /**
1415  * Creates the reactions for the model.
1416  */
1418 {
1419  // make sure the SBML Document already exists and that it has a Model set
1420  if (dataModel.getModel() == NULL || this->mpSBMLDocument == NULL || this->mpSBMLDocument->getModel() == NULL) return;
1421 
1422  CCopasiVectorNS<CReaction>::const_iterator it = dataModel.getModel()->getReactions().begin(), endit = dataModel.getModel()->getReactions().end();
1423 
1424  while (it != endit)
1425  {
1426  createReaction(**it, dataModel);
1427  ++it;
1428  }
1429 }
1430 
1431 /**
1432  * Creates the reaction for the given COPASI reaction.
1433  */
1435 {
1436  // make sure the mCOPASI2SBMLMap is already updated
1437  Reaction* pSBMLReaction = NULL;
1438 
1439  // if the reaction has nothing set but the name, we don't do anything
1440  // This is mandated since SBML Level 2 Version 2. A reactions has to have either
1441  // a substrate or a product.
1442  if (reaction.getChemEq().getSubstrates().size() == 0 &&
1443  reaction.getChemEq().getProducts().size() == 0) return;
1444 
1445  std::string sbmlId = reaction.getSBMLId();
1446 
1447  if (!sbmlId.empty())
1448  {
1449  pSBMLReaction = this->mpSBMLDocument->getModel()->getReaction(sbmlId);
1450 
1451  if (pSBMLReaction == NULL)
1452  {
1453  /* create a new reaction object */
1454  pSBMLReaction = this->mpSBMLDocument->getModel()->createReaction();
1455  mCOPASI2SBMLMap[&reaction] = pSBMLReaction;
1456  pSBMLReaction->setId(sbmlId);
1457  }
1458 
1459  // maybe this id was assigned by assignSBMLIdstoReactions and there is no object associated with it
1460  // If this is the case, we associate the object here.
1461  std::map<std::string, const SBase*>::const_iterator pos = this->mIdMap.find(sbmlId);
1462  assert(pos != this->mIdMap.end());
1463 
1464  if (pos->second == NULL)
1465  {
1466  this->mIdMap[sbmlId] = pSBMLReaction;
1467  }
1468  }
1469  else
1470  {
1471  pSBMLReaction = this->mpSBMLDocument->getModel()->createReaction();
1472  this->mCOPASI2SBMLMap[&reaction] = pSBMLReaction;
1473  sbmlId = CSBMLExporter::createUniqueId(this->mIdMap, reaction.getObjectName(), false);
1474  reaction.setSBMLId(sbmlId);
1475  pSBMLReaction->setId(sbmlId);
1476  }
1477 
1478  INIT_DEFAULTS((*pSBMLReaction));
1479 
1480  this->mIdMap.insert(std::pair<const std::string, const SBase*>(sbmlId, pSBMLReaction));
1481  this->mHandledSBMLObjects.insert(pSBMLReaction);
1482 
1483  // don't call setName on level 1 objects because this will also
1484  // change the id
1485  if (this->mpSBMLDocument->getLevel() > 1)
1486  {
1487  pSBMLReaction->setName(reaction.getObjectName().c_str());
1488  }
1489 
1490  pSBMLReaction->setReversible(reaction.isReversible());
1491  const CChemEq& chemicalEquation = reaction.getChemEq();
1492  /* Add all substrates */
1493  unsigned int counter;
1494  std::set<std::string> usedReferences;
1495 
1496  for (counter = 0; counter < chemicalEquation.getSubstrates().size(); counter++)
1497  {
1498  const CChemEqElement* element = chemicalEquation.getSubstrates()[counter];
1499  const CMetab* pMetabolite = element->getMetabolite();
1500  assert(pMetabolite);
1501  SpeciesReference* sRef = NULL;
1502 
1503  if (!(sRef = pSBMLReaction->getReactant(pMetabolite->getSBMLId())))
1504  {
1505  sRef = pSBMLReaction->createReactant();
1506  sRef->setSpecies(pMetabolite->getSBMLId().c_str());
1507  }
1508 
1509  INIT_DEFAULTS((*sRef));
1510 #if LIBSBML_VERSION > 40100
1511 
1512  if (this->mSBMLLevel > 2)
1513  sRef->setConstant(true);
1514 
1515 #endif
1516  sRef->setStoichiometry(element->getMultiplicity());
1517  sRef->setDenominator(1);
1518  usedReferences.insert(sRef->getSpecies());
1519  }
1520 
1521  ListOf* l = pSBMLReaction->getListOfReactants();
1522 
1523  for (counter = l->size(); counter > 0; --counter)
1524  {
1525  if (usedReferences.find(static_cast<SimpleSpeciesReference*>(l->get(counter - 1))->getSpecies()) == usedReferences.end())
1526  {
1527  l->remove(counter - 1);
1528  }
1529  }
1530 
1531  /* Add all products */
1532  usedReferences.clear();
1533 
1534  for (counter = 0; counter < chemicalEquation.getProducts().size(); counter++)
1535  {
1536  const CChemEqElement* element = chemicalEquation.getProducts()[counter];
1537  const CMetab* pMetabolite = element->getMetabolite();
1538  assert(pMetabolite);
1539  SpeciesReference* sRef = NULL;
1540 
1541  if (!(sRef = pSBMLReaction->getProduct(pMetabolite->getSBMLId())))
1542  {
1543  sRef = pSBMLReaction->createProduct();
1544  sRef->setSpecies(pMetabolite->getSBMLId().c_str());
1545  }
1546 
1547  INIT_DEFAULTS((*sRef));
1548 #if LIBSBML_VERSION > 40100
1549 
1550  if (this->mSBMLLevel > 2)
1551  sRef->setConstant(true);
1552 
1553 #endif
1554 
1555  sRef->setStoichiometry(element->getMultiplicity());
1556  sRef->setDenominator(1);
1557  usedReferences.insert(sRef->getSpecies());
1558  }
1559 
1560  l = pSBMLReaction->getListOfProducts();
1561 
1562  for (counter = l->size(); counter > 0; --counter)
1563  {
1564  if (usedReferences.find(static_cast<SimpleSpeciesReference*>(l->get(counter - 1))->getSpecies()) == usedReferences.end())
1565  {
1566  l->remove(counter - 1);
1567  }
1568  }
1569 
1570  /* Add all modifiers */
1571  usedReferences.clear();
1572 
1573  for (counter = 0; counter < chemicalEquation.getModifiers().size(); counter++)
1574  {
1575  const CChemEqElement* element = chemicalEquation.getModifiers()[counter];
1576  const CMetab* pMetabolite = element->getMetabolite(); assert(pMetabolite);
1577  ModifierSpeciesReference* sRef = NULL;
1578 
1579  if (!(sRef = pSBMLReaction->getModifier(pMetabolite->getSBMLId())))
1580  {
1581  if (pSBMLReaction->getLevel() > 1)
1582  {
1583  sRef = pSBMLReaction->createModifier();
1584  assert(sRef != NULL);
1585  sRef->setSpecies(pMetabolite->getSBMLId().c_str());
1586  }
1587  }
1588 
1589  if (pSBMLReaction->getLevel() > 1)
1590  {
1591  usedReferences.insert(sRef->getSpecies());
1592  }
1593  }
1594 
1595  l = pSBMLReaction->getListOfModifiers();
1596 
1597  for (counter = l->size(); counter > 0; --counter)
1598  {
1599  if (usedReferences.find(static_cast<SimpleSpeciesReference*>(l->get(counter - 1))->getSpecies()) == usedReferences.end())
1600  {
1601  l->remove(counter - 1);
1602  }
1603  }
1604 
1605  /* create the kinetic law */
1606 
1607  /* if there is one on COPASI */
1609  {
1610  // make sure the creaated kinetic law has the same level and version as the reaction
1611  // we intend to add it to, otherwise the setKineticLaw function of libsbml 4 will fail
1612  // to set the kinetic law
1613  KineticLaw* pKineticLaw = this->createKineticLaw(reaction, dataModel, pSBMLReaction->getLevel(), pSBMLReaction->getVersion());
1614 
1615  if (pKineticLaw != NULL)
1616  {
1617  pSBMLReaction->setKineticLaw(pKineticLaw);
1618  delete pKineticLaw;
1619  }
1620  else
1621  {
1622  // create an error message and abort
1623  // maybe only abort if incomplete export is not wanted
1624  if (this->mIncompleteExport != true)
1625  {
1627  }
1628  else
1629  {
1630  pSBMLReaction->unsetKineticLaw();
1631  }
1632  }
1633  }
1634  else
1635  {
1636  pSBMLReaction->unsetKineticLaw();
1637  }
1638 
1639  if (pSBMLReaction != NULL)
1640  {
1641  CSBMLExporter::setSBMLNotes(pSBMLReaction, &reaction);
1642  }
1643 
1644  CSBMLExporter::updateMIRIAMAnnotation(&reaction, pSBMLReaction, this->mMetaIdMap);
1645 }
1646 
1647 /**
1648  * Creates the initial assignments for the model.
1649  */
1651 {
1652  // make sure the mInitialAssignmentVector has been filled already
1653 
1654  // create the initial assignments
1655  size_t i, iMax = this->mInitialAssignmentVector.size();
1656  const CModelEntity* pME = NULL;
1657 
1658  for (i = 0; i < iMax; ++i)
1659  {
1660  pME = mInitialAssignmentVector[i];
1661  assert(pME != NULL);
1662 
1663  if (pME != NULL)
1664  {
1665  createInitialAssignment(*pME, dataModel);
1666  }
1667  }
1668 }
1669 
1670 /**
1671  * Creates the initial assignment for the given COPASI model entity.
1672  */
1674 {
1675  // check the expression
1676  std::vector<SBMLIncompatibility> result;
1678  , dataModel
1679  , this->mSBMLLevel
1680  , this->mSBMLVersion
1681  , result
1682  , std::string("initial expression for object named \"" + modelEntity.getObjectName() + "\"").c_str()
1683  , true
1684  , &mInitialValueMap);
1685 
1686  // collect directly used functions
1687  if (result.empty())
1688  {
1689  std::set<std::string> directlyUsedFunctionNames;
1690  CSBMLExporter::findDirectlyUsedFunctions(modelEntity.getInitialExpressionPtr()->getRoot(), directlyUsedFunctionNames);
1691  std::set<CFunction*> usedFunctions = CSBMLExporter::createFunctionSetFromFunctionNames(directlyUsedFunctionNames, CCopasiRootContainer::getFunctionList());
1692 
1693 #if defined _MSC_VER && _MSC_VER < 1201 // 1200 Identifies Visual C++ 6.0
1694  {
1695  std::set<CFunction*>::const_iterator it = usedFunctions.begin();
1696  std::set<CFunction*>::const_iterator end = usedFunctions.end();
1697 
1698  for (; it != end; ++it)
1699  this->mUsedFunctions.insert(*it);
1700  }
1701 #else
1702  this->mUsedFunctions.insert(usedFunctions.begin(), usedFunctions.end());
1703 #endif
1704 
1705  // create the actual initial assignment
1706  InitialAssignment* pInitialAssignment = this->mpSBMLDocument->getModel()->getInitialAssignment(modelEntity.getSBMLId());
1707 
1708  if (pInitialAssignment == NULL)
1709  {
1710  pInitialAssignment = this->mpSBMLDocument->getModel()->createInitialAssignment();
1711  pInitialAssignment->setSymbol(modelEntity.getSBMLId());
1712  }
1713 
1714  // set the math
1715  this->mHandledSBMLObjects.insert(pInitialAssignment);
1716  const CEvaluationNode* pOrigNode = modelEntity.getInitialExpressionPtr()->getRoot();
1717 
1719  {
1720  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 70, "initial assignment", modelEntity.getObjectType().c_str(), modelEntity.getObjectName().c_str());
1721  }
1722 
1723  // the next few lines replace references to species depending on whether
1724  // it is a reference to an amount or a reference to a concentration.
1725  // Other factors that influence this replacement are if the model
1726  // contains variable volumes or if the quantity units are set to CModel::number
1727  pOrigNode = this->replaceSpeciesReferences(pOrigNode, dataModel);
1728  assert(pOrigNode != NULL);
1729 
1730  // check if the rule is for an amount species
1731  // if this is the case, we have to multiply the expression by the volume
1732  // of the compartment the species is in
1733  // create the actual rule
1734  const CMetab* pMetab = dynamic_cast<const CMetab*>(&modelEntity);
1735 
1736  if (pMetab != NULL)
1737  {
1738  std::map<const CCopasiObject*, SBase*>::const_iterator pos = this->mCOPASI2SBMLMap.find(&modelEntity);
1739  assert(pos != this->mCOPASI2SBMLMap.end());
1740 
1741  if (dynamic_cast<const Species*>(pos->second)->getHasOnlySubstanceUnits() == true)
1742  {
1743  const CCompartment* pCompartment = pMetab->getCompartment();
1744 
1745  if (pCompartment->getDimensionality() != 0)
1746  {
1747  CEvaluationNode* pNode = CSBMLExporter::multiplyByObject(pOrigNode, pCompartment->getInitialValueReference());
1748  assert(pNode != NULL);
1749 
1750  if (pNode != NULL)
1751  {
1752  delete pOrigNode;
1753  pOrigNode = pNode;
1754  }
1755  }
1756  }
1757  }
1758 
1759  ASTNode* pNode = this->convertToASTNode(pOrigNode, dataModel);
1760  delete pOrigNode;
1761  // convert local parameters that are referenced in an expression to global
1762  // parameters
1763  this->replace_local_parameters(pNode, dataModel);
1764 
1765  if (pNode != NULL)
1766  {
1767  pInitialAssignment->setMath(pNode);
1768  delete pNode;
1769  }
1770  else
1771  {
1772  if (this->mIncompleteExport == true)
1773  {
1774  // remove the initial assignment from the SBML model
1775  unsigned int i = 0, iMax = this->mpSBMLDocument->getModel()->getNumInitialAssignments();
1776 
1777  while (i < iMax)
1778  {
1779  if (this->mpSBMLDocument->getModel()->getInitialAssignment(i)->getSymbol() == modelEntity.getSBMLId())
1780  {
1781  this->mpSBMLDocument->getModel()->getListOfInitialAssignments()->remove(i);
1782  break;
1783  }
1784 
1785  ++i;
1786  }
1787  }
1788  else
1789  {
1790  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 60, "initial assignment", modelEntity.getObjectType().c_str(), modelEntity.getObjectName().c_str());
1791  }
1792  }
1793  }
1794  else
1795  {
1796  this->mIncompatibilities.insert(this->mIncompatibilities.end(), result.begin(), result.end());
1797 
1798  if (!this->mIncompleteExport)
1799  {
1800  this->outputIncompatibilities();
1801  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 60, "initial assignment", modelEntity.getObjectType().c_str(), modelEntity.getObjectName().c_str());
1802  }
1803  }
1804 }
1805 
1806 /*
1807  * This function replaces all references to initial values, volumes and concentrations with the
1808  * name of the global parameter that will be created.
1809  */
1810 std::string convertExpression(const std::string& expression, const std::map<const std::string, Parameter*>& initialValueMap)
1811 {
1812  if (initialValueMap.empty())
1813  return expression;
1814 
1815  std::string result = expression;
1816  std::map<const std::string, Parameter*>::const_iterator it;
1817 
1818  for (it = initialValueMap.begin(); it != initialValueMap.end(); ++it)
1819  {
1820  size_t length = it->first.length();
1821  size_t start;
1822 
1823  while ((start = result.find(it->first)) != std::string::npos)
1824  {
1825  result.replace(start, length, it->second->getId());
1826  }
1827  }
1828 
1829  return result;
1830 }
1831 
1832 /**
1833  * Creates the rules for the model.
1834  */
1836 {
1837  // make sure the mAssignmentVector has been filled already
1838  // make sure the mODEVector has been filled already
1839  // order the rules for Level 1 export
1840  // rules in Level 2 are not ordered.
1841  //if (this->mSBMLLevel == 1 || (this->mSBMLLevel == 2 && this->mSBMLVersion == 1))
1842  // {
1843  std::vector<CModelEntity*> orderedAssignmentRules = orderRules(dataModel);
1844  //}
1845  //
1846  // remove all rules from the SBML model
1847  // and put them in a map
1848  std::map<std::string, Rule*> ruleMap;
1849  Model* pSBMLModel = this->mpSBMLDocument->getModel();
1850  assert(pSBMLModel != NULL);
1851 
1852  while (pSBMLModel->getNumRules() != 0)
1853  {
1854  Rule* pRule = pSBMLModel->getRule(0);
1855  assert(pRule != NULL);
1856 
1857  // ignore algebraic rules
1858  if (pRule->getTypeCode() == SBML_ASSIGNMENT_RULE)
1859  {
1860  AssignmentRule* pARule = dynamic_cast<AssignmentRule*>(pRule);
1861  assert(pARule != NULL);
1862  ruleMap[pARule->getVariable()] = pARule;
1863  }
1864  else if (pRule->getTypeCode() == SBML_RATE_RULE)
1865  {
1866  RateRule* pRRule = dynamic_cast<RateRule*>(pRule);
1867  assert(pRRule != NULL);
1868  ruleMap[pRRule->getVariable()] = pRRule;
1869  }
1870 
1871  pSBMLModel->getListOfRules()->remove(0);
1872  }
1873 
1874  // create the rules
1875  const CModelEntity* pME = NULL;
1876  size_t i, iMax = orderedAssignmentRules.size();
1877 
1878  for (i = 0; i < iMax; ++i)
1879  {
1880  pME = orderedAssignmentRules[i];
1881  assert(pME != NULL);
1882 
1883  if (pME != NULL)
1884  {
1885  std::map<std::string, Rule*>::const_iterator pos = ruleMap.find(pME->getSBMLId());
1886  Rule* pOldRule = NULL;
1887 
1888  if (pos != ruleMap.end())
1889  {
1890  pOldRule = pos->second;
1891  assert(pOldRule != NULL);
1892  ruleMap.erase(pos->first);
1893  }
1894 
1895  createRule(*pME, dataModel, pOldRule);
1896  }
1897  }
1898 
1899  iMax = mODEVector.size();
1900 
1901  for (i = 0; i < iMax; ++i)
1902  {
1903  pME = mODEVector[i];
1904  assert(pME != NULL);
1905 
1906  if (pME != NULL)
1907  {
1908  std::map<std::string, Rule*>::const_iterator pos = ruleMap.find(pME->getSBMLId());
1909  Rule* pOldRule = NULL;
1910 
1911  if (pos != ruleMap.end())
1912  {
1913  pOldRule = pos->second;
1914  assert(pOldRule != NULL);
1915  ruleMap.erase(pos->first);
1916  }
1917 
1918  createRule(*pME, dataModel, pOldRule);
1919  }
1920  }
1921 
1922  // delete the remaining rules
1923  std::map<std::string, Rule*>::iterator mapIt = ruleMap.begin(), mapEndit = ruleMap.end();
1924 
1925  while (mapIt != mapEndit)
1926  {
1927  delete mapIt->second;
1928  ++mapIt;
1929  }
1930 
1931  ruleMap.clear();
1932 }
1933 
1934 /**
1935  * Creates the rule for the given COPASI model entity.
1936  */
1937 void CSBMLExporter::createRule(const CModelEntity& modelEntity, CCopasiDataModel& dataModel, Rule* pOldRule)
1938 {
1939  // check the expression
1940  std::vector<SBMLIncompatibility> result;
1942  , dataModel
1943  , this->mSBMLLevel
1944  , this->mSBMLVersion
1945  , result
1946  , std::string("rule for object named \"" + modelEntity.getObjectName() + "\"").c_str()
1947  , false
1948  , &mInitialValueMap);
1949 
1950  // collect directly used functions
1951  if (result.empty())
1952  {
1953  std::set<std::string> directlyUsedFunctionNames;
1954  CSBMLExporter::findDirectlyUsedFunctions(modelEntity.getExpressionPtr()->getRoot(), directlyUsedFunctionNames);
1955  std::set<CFunction*> usedFunctions = CSBMLExporter::createFunctionSetFromFunctionNames(directlyUsedFunctionNames, CCopasiRootContainer::getFunctionList());
1956 
1957 # if defined _MSC_VER && _MSC_VER < 1201 // 1200 Identifies Visual C++ 6.0
1958  {
1959  std::set<CFunction*>::const_iterator it = usedFunctions.begin();
1960  std::set<CFunction*>::const_iterator end = usedFunctions.end();
1961 
1962  for (; it != end; ++it)
1963  this->mUsedFunctions.insert(*it);
1964  }
1965 #else
1966  this->mUsedFunctions.insert(usedFunctions.begin(), usedFunctions.end());
1967 #endif
1968 
1969  // create the actual rule
1970  //Rule* pRule = this->mpSBMLDocument->getModel()->getRule(modelEntity.getSBMLId());
1971  const CMetab* pMetab = dynamic_cast<const CMetab*>(&modelEntity);
1972 
1973  if (pOldRule == NULL)
1974  {
1975  if (modelEntity.getStatus() == CModelEntity::ASSIGNMENT)
1976  {
1977  pOldRule = this->mpSBMLDocument->getModel()->createAssignmentRule();
1978  }
1979  else
1980  {
1981  if (pMetab != NULL)
1982  {
1983  // check if the compartment is fixed
1984  if (pMetab->getCompartment()->getStatus() != CModelEntity::FIXED)
1985  {
1986  CCopasiMessage(CCopasiMessage::ERROR, MCSBML + 52, pMetab->getObjectName().c_str());
1987  }
1988  }
1989 
1990  pOldRule = this->mpSBMLDocument->getModel()->createRateRule();
1991  }
1992 
1993  pOldRule->setVariable(modelEntity.getSBMLId());
1994  }
1995  else
1996  {
1997  // Read the rule to the model
1998  this->mpSBMLDocument->getModel()->getListOfRules()->appendAndOwn(pOldRule);
1999  }
2000 
2001  // set the math
2002 
2003  const std::string& changedExpression = convertExpression(modelEntity.getExpression(), mInitialValueMap);
2004  CEvaluationTree tree;
2005  tree.setInfix(changedExpression);
2006  const CEvaluationNode* pOrigNode = tree.getRoot();
2007 
2009  {
2010  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 70, "assignment", modelEntity.getObjectType().c_str(), modelEntity.getObjectName().c_str());
2011  }
2012 
2013  // the next few lines replace references to species depending on whether
2014  // it is a reference to an amount or a reference to a concentration.
2015  // Other factors that influence this replacement are if the model
2016  // contains variable volumes or if the quantity units are set to CModel::number
2017  pOrigNode = this->replaceSpeciesReferences(pOrigNode, dataModel);
2018  assert(pOrigNode != NULL);
2019 
2020  // check if the rule is for an amount species
2021  // if this is the case, we have to multiply the expression by the volume
2022  // of the compartment the species is in
2023  if (pMetab != NULL)
2024  {
2025  std::map<const CCopasiObject*, SBase*>::const_iterator pos = this->mCOPASI2SBMLMap.find(&modelEntity);
2026  assert(pos != this->mCOPASI2SBMLMap.end());
2027 
2028  if (dynamic_cast<const Species*>(pos->second)->getHasOnlySubstanceUnits() == true)
2029  {
2030  const CCompartment* pCompartment = pMetab->getCompartment();
2031 
2032  if (pCompartment->getDimensionality() != 0)
2033  {
2034  CEvaluationNode* pNode = CSBMLExporter::multiplyByObject(pOrigNode, pCompartment->getValueReference());
2035  assert(pNode != NULL);
2036 
2037  if (pNode != NULL)
2038  {
2039  delete pOrigNode;
2040  pOrigNode = pNode;
2041  }
2042  }
2043  }
2044  }
2045 
2046  ASTNode* pNode = this->convertToASTNode(pOrigNode, dataModel);
2047  this->replace_local_parameters(pNode, dataModel);
2048  delete pOrigNode;
2049 
2050  if (pNode != NULL)
2051  {
2052  pOldRule->setMath(pNode);
2053  delete pNode;
2054  }
2055  else
2056  {
2057  if (this->mIncompleteExport != true)
2058  {
2059  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 60, "rule", modelEntity.getObjectType().c_str(), modelEntity.getObjectName().c_str());
2060  }
2061  }
2062  }
2063  else
2064  {
2065  this->mIncompatibilities.insert(this->mIncompatibilities.end(), result.begin(), result.end());
2066 
2067  if (!this->mIncompleteExport)
2068  {
2069  this->outputIncompatibilities();
2070  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 60, "rule", modelEntity.getObjectType().c_str(), modelEntity.getObjectName().c_str());
2071  }
2072  }
2073 }
2074 
2075 const std::map<std::string, const SBase*> CSBMLExporter::createIdMap(const Model& sbmlModel)
2076 {
2077  // go through all function definitions, compartments, species, reactions,
2078  // parameters and events and create a set with all used ids.
2079  std::map<std::string, const SBase*> idMap;
2080  unsigned int i, iMax;
2081 
2082  if (sbmlModel.isSetId())
2083  {
2084  idMap.insert(std::pair<const std::string, const SBase*>(sbmlModel.getId(), &sbmlModel));
2085  }
2086 
2087  if (sbmlModel.getListOfFunctionDefinitions()->isSetId())
2088  {
2089  idMap.insert(std::pair<const std::string, const SBase*>(sbmlModel.getListOfFunctionDefinitions()->getId(), sbmlModel.getListOfFunctionDefinitions()));
2090  }
2091 
2092  iMax = sbmlModel.getNumFunctionDefinitions();
2093 
2094  for (i = 0; i < iMax; ++i)
2095  {
2096  idMap.insert(std::pair<const std::string, const SBase*>(sbmlModel.getFunctionDefinition(i)->getId(), sbmlModel.getFunctionDefinition(i)));
2097  }
2098 
2099  if (sbmlModel.getListOfCompartments()->isSetId())
2100  {
2101  idMap.insert(std::pair<const std::string, const SBase*>(sbmlModel.getListOfCompartments()->getId(), sbmlModel.getListOfCompartments()));
2102  }
2103 
2104  iMax = sbmlModel.getNumCompartments();
2105 
2106  for (i = 0; i < iMax; ++i)
2107  {
2108  idMap.insert(std::pair<const std::string, const SBase*>(sbmlModel.getCompartment(i)->getId(), sbmlModel.getCompartment(i)));
2109  }
2110 
2111  if (sbmlModel.getListOfSpecies()->isSetId())
2112  {
2113  idMap.insert(std::pair<const std::string, const SBase*>(sbmlModel.getListOfSpecies()->getId(), sbmlModel.getListOfSpecies()));
2114  }
2115 
2116  iMax = sbmlModel.getNumSpecies();
2117 
2118  for (i = 0; i < iMax; ++i)
2119  {
2120  idMap.insert(std::pair<const std::string, const SBase*>(sbmlModel.getSpecies(i)->getId(), sbmlModel.getSpecies(i)));
2121  }
2122 
2123  if (sbmlModel.getListOfParameters()->isSetId())
2124  {
2125  idMap.insert(std::pair<const std::string, const SBase*>(sbmlModel.getListOfParameters()->getId(), sbmlModel.getListOfParameters()));
2126  }
2127 
2128  iMax = sbmlModel.getNumParameters();
2129 
2130  for (i = 0; i < iMax; ++i)
2131  {
2132  idMap.insert(std::pair<const std::string, const SBase*>(sbmlModel.getParameter(i)->getId(), sbmlModel.getParameter(i)));
2133  }
2134 
2135  if (sbmlModel.getListOfReactions()->isSetId())
2136  {
2137  idMap.insert(std::pair<const std::string, const SBase*>(sbmlModel.getListOfReactions()->getId(), sbmlModel.getListOfReactions()));
2138  }
2139 
2140  iMax = sbmlModel.getNumReactions();
2141 
2142  for (i = 0; i < iMax; ++i)
2143  {
2144  const Reaction* pReaction = sbmlModel.getReaction(i);
2145 
2146  if (pReaction != NULL)
2147  {
2148  idMap.insert(std::pair<const std::string, const SBase*>(pReaction->getId(), pReaction));
2149 
2150  if (pReaction->getListOfReactants()->isSetId())
2151  {
2152  idMap.insert(std::pair<const std::string, const SBase*>(pReaction->getListOfReactants()->getId(), pReaction->getListOfReactants()));
2153  }
2154 
2155  unsigned int j, jMax = pReaction->getNumReactants();
2156  const SpeciesReference* pSRef = NULL;
2157 
2158  for (j = 0; j < jMax; ++j)
2159  {
2160  pSRef = pReaction->getReactant(j);
2161 
2162  if (pSRef->isSetId())
2163  {
2164  idMap.insert(std::pair<const std::string, const SBase*>(pSRef->getId(), pSRef));
2165  }
2166  }
2167 
2168  if (pReaction->getListOfProducts()->isSetId())
2169  {
2170  idMap.insert(std::pair<const std::string, const SBase*>(pReaction->getListOfProducts()->getId(), pReaction->getListOfProducts()));
2171  }
2172 
2173  jMax = pReaction->getNumProducts();
2174 
2175  for (j = 0; j < jMax; ++j)
2176  {
2177  pSRef = pReaction->getProduct(j);
2178 
2179  if (pSRef->isSetId())
2180  {
2181  idMap.insert(std::pair<const std::string, const SBase*>(pSRef->getId(), pSRef));
2182  }
2183  }
2184 
2185  if (pReaction->getListOfModifiers()->isSetId())
2186  {
2187  idMap.insert(std::pair<const std::string, const SBase*>(pReaction->getListOfModifiers()->getId(), pReaction->getListOfModifiers()));
2188  }
2189 
2190  jMax = pReaction->getNumModifiers();
2191  const ModifierSpeciesReference* pModSRef;
2192 
2193  for (j = 0; j < jMax; ++j)
2194  {
2195  pModSRef = pReaction->getModifier(j);
2196 
2197  if (pModSRef->isSetId())
2198  {
2199  idMap.insert(std::pair<const std::string, const SBase*>(pModSRef->getId(), pModSRef));
2200  }
2201  }
2202  }
2203  }
2204 
2205  if (sbmlModel.getListOfEvents()->isSetId())
2206  {
2207  idMap.insert(std::pair<const std::string, const SBase*>(sbmlModel.getListOfEvents()->getId(), sbmlModel.getListOfEvents()));
2208  }
2209 
2210  iMax = sbmlModel.getNumEvents();
2211 
2212  for (i = 0; i < iMax; ++i)
2213  {
2214  idMap.insert(std::pair<const std::string, const SBase*>(sbmlModel.getEvent(i)->getId(), sbmlModel.getEvent(i)));
2215  }
2216 
2217  // if COPASI is compiled with layout, we have to add those ids as well
2218  const LayoutModelPlugin* lmPlugin = (LayoutModelPlugin*)sbmlModel.getPlugin("layout");
2219 
2220  if (lmPlugin != NULL)
2221  {
2222  if (lmPlugin->getListOfLayouts()->isSetId())
2223  {
2224  idMap.insert(std::pair<const std::string, const SBase*>(lmPlugin->getListOfLayouts()->getId(), lmPlugin->getListOfLayouts()));
2225  }
2226 
2227  iMax = lmPlugin->getListOfLayouts()->size();
2228 
2229  for (i = 0; i < iMax; ++i)
2230  {
2231  const Layout* pLayout = lmPlugin->getLayout(i);
2232 
2233  if (pLayout != NULL)
2234  {
2235  if (pLayout->isSetId())
2236  {
2237  idMap.insert(std::pair<const std::string, const SBase*>(pLayout->getId(), pLayout));
2238  }
2239  }
2240  }
2241  }
2242 
2243  return idMap;
2244 }
2245 
2246 /**
2247  * Create a unique id for an SBML object.
2248  * I can't just take the COPASI key of the object since this might conflict
2249  * with an already existing SBML id which came from the sbmlid attribute in a
2250  * COPASI file or directly by importing an SBML file.
2251  */
2252 const std::string CSBMLExporter::createUniqueId(const std::map < std::string, const SBase* > & idMap,
2253  const std::string& prefix,
2254  bool addIndexForFirst,
2255  const std::string & separator)
2256 {
2257  std::string Prefix = nameToSbmlId(prefix);
2258 
2259  unsigned int i = 0;
2260  std::ostringstream Id;
2261 
2262  if (addIndexForFirst)
2263  {
2264  Id << Prefix << separator << i++;
2265  }
2266  else
2267  {
2268  Id << Prefix;
2269  }
2270 
2271  while (idMap.find(Id.str()) != idMap.end())
2272  {
2273  Id.str("");
2274  Id << Prefix << separator << i++;
2275  }
2276 
2277  return Id.str();
2278 }
2279 
2280 /**
2281  * Checks all assignments (initial and transient) for references to objects
2282  * that can not be exported to SBML.
2283  */
2284 void CSBMLExporter::checkForUnsupportedObjectReferences(const CCopasiDataModel& dataModel, unsigned int sbmlLevel, unsigned int sbmlVersion, std::vector<SBMLIncompatibility>& result)
2285 {
2286  // check all metabolites,parameters and compartments
2287  const CModel* pModel = dataModel.getModel();
2288  assert(pModel);
2289 
2290  if (pModel == NULL) return;
2291 
2292  size_t i, iMax = this->mAssignmentVector.size();
2293  const CModelEntity* pME;
2294 
2295  for (i = 0; i < iMax; ++i)
2296  {
2297  pME = this->mAssignmentVector[i];
2298  assert(pME != NULL);
2299 
2300  if (pME != NULL)
2301  {
2302  checkForUnsupportedObjectReferences(*pME->getExpressionPtr(), dataModel, sbmlLevel, sbmlVersion, result, false, &mInitialValueMap);
2303  }
2304  }
2305 
2306  // check ode rules
2307  iMax = this->mODEVector.size();
2308 
2309  for (i = 0; i < iMax; ++i)
2310  {
2311  pME = this->mODEVector[i];
2312  assert(pME != NULL);
2313 
2314  if (pME != NULL)
2315  {
2316  checkForUnsupportedObjectReferences(*pME->getExpressionPtr(), dataModel, sbmlLevel, sbmlVersion, result, false, &mInitialValueMap);
2317  }
2318  }
2319 
2320  // check initial assignments
2321  iMax = this->mInitialAssignmentVector.size();
2322 
2323  for (i = 0; i < iMax; ++i)
2324  {
2325  pME = this->mInitialAssignmentVector[i];
2326  assert(pME != NULL);
2327 
2328  if (pME != NULL)
2329  {
2330  checkForUnsupportedObjectReferences(*pME->getInitialExpressionPtr(), dataModel, sbmlLevel, sbmlVersion, result, false, &mInitialValueMap);
2331  }
2332  }
2333 }
2334 
2335 std::string getAnnotationStringFor(const CCopasiObject* pObjectParent)
2336 {
2337  std::stringstream str;
2338  str << "<initialValue xmlns='http://copasi.org/initialValue' ";
2339  str << "parent='" << ((CModelEntity*)pObjectParent)->getSBMLId() << "' />";
2340  return str.str();
2341 }
2342 
2343 /*
2344  * Adds the given object to the initialMap, and creates a new parameter
2345  * for it.
2346  */
2347 void addToInitialValueMap(std::map<const std::string, Parameter*>* initialMap
2348  , const CCopasiObject* pObject
2349  , const CCopasiObject* pObjectParent
2350  , int sbmlLevel
2351  , int sbmlVersion)
2352 {
2353  if (initialMap == NULL || pObject == NULL || pObjectParent == NULL)
2354  return;
2355 
2356  const std::string& cn = pObject->getCN();
2357 
2358  if ((*initialMap)[cn] != NULL)
2359  {
2360  // already have the initial assignment no need to add another one
2361  return;
2362  }
2363 
2364  Parameter* initial = new Parameter(sbmlLevel, sbmlVersion);
2365  initial->setAnnotation(getAnnotationStringFor(pObjectParent));
2366  initial->initDefaults();
2367  initial->setId(pObjectParent->getKey());
2368  initial->setName("Initial for " + pObjectParent->getObjectName());
2369 
2370  if (pObject->isValueDbl())
2371  initial->setValue(*((const C_FLOAT64*) pObject->getValuePointer()));
2372 
2373  (*initialMap) [cn] = initial;
2374 }
2375 
2377  const CEvaluationTree& expr
2378  , const CCopasiDataModel& dataModel
2379  , unsigned int sbmlLevel
2380  , unsigned int sbmlVersion
2381  , std::vector<SBMLIncompatibility>& result
2382  , bool initialExpression /*= false*/
2383  , std::map<const std::string, Parameter*>* initialMap /*= NULL*/)
2384 {
2385  // SBML Level 1 and Level 2 Version 1 can have references to transient values of
2386  // compartments, metabolites and global parameters as well as time (model)
2387 
2388  // SBML Level 2 Version 2 and above can have the same references as level 1 plus references to
2389  // reaction fluxes.
2390 
2391  // Starting with COPASI Build 30 references to local parameters are allowed
2392  // as well because they are converted to global parameters
2393  const std::vector<CEvaluationNode*>& objectNodes = expr.getNodeList();
2394  size_t j, jMax = objectNodes.size();
2395 
2396  for (j = 0; j < jMax; ++j)
2397  {
2398  if (CEvaluationNode::type(objectNodes[j]->getType()) == CEvaluationNode::OBJECT)
2399  {
2400  const CEvaluationNodeObject* pObjectNode = dynamic_cast<const CEvaluationNodeObject*>(objectNodes[j]);
2401  assert(pObjectNode);
2402 
2403  if (pObjectNode == NULL) continue;
2404 
2405  std::vector<CCopasiContainer*> containers;
2406  containers.push_back(const_cast<CModel*>(dataModel.getModel()));
2407  const CCopasiObject* pObject =
2408  dataModel.ObjectFromName(containers, pObjectNode->getObjectCN());
2409  assert(pObject);
2410 
2411  if (pObject->isReference())
2412  {
2413  const CCopasiObject* pObjectParent = pObject->getObjectParent();
2414  assert(pObjectParent);
2415  std::string typeString = pObjectParent->getObjectType();
2416 
2417  if (typeString == "Compartment")
2418  {
2419  // must be a reference to the (transient) or initial volume
2420  if (initialExpression == true)
2421  {
2422  if (pObject->getObjectName() != "InitialVolume")
2423  {
2424  result.push_back(SBMLIncompatibility(1, pObject->getObjectName().c_str(), "compartment", pObjectParent->getObjectName().c_str()));
2425  }
2426  }
2427  else
2428  {
2429  if (pObject->getObjectName() == "InitialVolume"
2430  && initialMap != NULL
2431  && (sbmlLevel > 2 || (sbmlLevel == 2 && sbmlVersion > 1)))
2432  {
2433  // add new parameter for the initial value
2434  addToInitialValueMap(initialMap, pObject, pObjectParent, sbmlLevel, sbmlVersion);
2435  }
2436  else if (pObject->getObjectName() != "Volume"
2437  && pObject->getObjectName() != "Rate")
2438  {
2439  result.push_back(SBMLIncompatibility(1, pObject->getObjectName().c_str(), "compartment", pObjectParent->getObjectName().c_str()));
2440  }
2441  }
2442  }
2443  else if (typeString == "Metabolite")
2444  {
2445  // must be a reference to the transient or initial concentration
2446  if (initialExpression == true)
2447  {
2448  if (pObject->getObjectName() != "InitialConcentration" && pObject->getObjectName() != "InitialParticleNumber")
2449  {
2450  result.push_back(SBMLIncompatibility(1, pObject->getObjectName().c_str(), "metabolite", pObjectParent->getObjectName().c_str()));
2451  }
2452  }
2453  else
2454  {
2455  if (pObject->getObjectName() == "InitialConcentration"
2456  && initialMap != NULL
2457  && (sbmlLevel > 2 || (sbmlLevel == 2 && sbmlVersion > 1)))
2458  {
2459  // add new parameter for the initial value
2460  addToInitialValueMap(initialMap, pObject, pObjectParent, sbmlLevel, sbmlVersion);
2461  }
2462  else if (pObject->getObjectName() != "Concentration"
2463  && pObject->getObjectName() != "ParticleNumber"
2464  && pObject->getObjectName() != "Rate")
2465  {
2466 
2467  result.push_back(SBMLIncompatibility(1, pObject->getObjectName().c_str(), "metabolite", pObjectParent->getObjectName().c_str()));
2468  }
2469  }
2470  }
2471  else if (typeString == "ModelValue")
2472  {
2473  // must be a reference to the transient or initial value
2474  if (initialExpression == true)
2475  {
2476  if (pObject->getObjectName() != "InitialValue")
2477  {
2478  result.push_back(SBMLIncompatibility(1, pObject->getObjectName().c_str(), "parameter", pObjectParent->getObjectName().c_str()));
2479  }
2480  }
2481  else
2482  {
2483  if (pObject->getObjectName() == "InitialValue"
2484  && initialMap != NULL
2485  && (sbmlLevel > 2 || (sbmlLevel == 2 && sbmlVersion > 1)))
2486  {
2487  // add new parameter for the initial value
2488  addToInitialValueMap(initialMap, pObject, pObjectParent, sbmlLevel, sbmlVersion);
2489  }
2490  else if (pObject->getObjectName() != "Value"
2491  && pObject->getObjectName() != "Rate")
2492  {
2493  result.push_back(SBMLIncompatibility(1, pObject->getObjectName().c_str(), "parameter", pObjectParent->getObjectName().c_str()));
2494  }
2495  }
2496  }
2497  else if (typeString == "Model")
2498  {
2499  // must be a reference to the model time
2500  if (pObject->getObjectName() != "Time")
2501  {
2502  result.push_back(SBMLIncompatibility(1, pObject->getObjectName().c_str(), "model", pObjectParent->getObjectName().c_str()));
2503  }
2504  }
2505  else if (typeString == "Parameter")
2506  {
2507  // must be a local parameter
2508  if (pObject->getObjectName() != "Value")
2509  {
2510  result.push_back(SBMLIncompatibility(1, pObject->getObjectName().c_str(), "local parameter", pObjectParent->getObjectName().c_str()));
2511  }
2512  }
2513  else
2514  {
2515  if (sbmlLevel == 1 || (sbmlLevel == 2 && sbmlVersion == 1))
2516  {
2517  result.push_back(SBMLIncompatibility(10, pObject->getObjectName().c_str(), typeString.c_str(), pObjectParent->getObjectName().c_str()));
2518  }
2519  else
2520  {
2521  if (typeString == "Reaction")
2522  {
2523  if (pObject->getObjectName() != "Flux")
2524  {
2525  result.push_back(SBMLIncompatibility(1, pObject->getObjectName().c_str(), "reaction", pObjectParent->getObjectName().c_str()));
2526  }
2527  }
2528  }
2529  }
2530  }
2531  else
2532  {
2533  // normally local parameters are referenced via their value and
2534  // not directly, so this code should never be called
2535  const CCopasiParameter* pLocalParameter = dynamic_cast<const CCopasiParameter*>(pObject);
2536 
2537  if (pLocalParameter == NULL)
2538  {
2539  result.push_back(SBMLIncompatibility(1, "value", pObject->getObjectType().c_str() , pObject->getObjectName().c_str()));
2540  }
2541  }
2542  }
2543  }
2544 }
2545 
2546 /**
2547  * Checks whether the given data model can be exported to SBML Level1
2548  * If it can be exported, the result vector will be empty, otherwise it will
2549  * contain a number of messages that specify why it can't be exported.
2550  */
2551 void CSBMLExporter::isModelSBMLL1Compatible(const CCopasiDataModel& dataModel, std::vector<SBMLIncompatibility>& result)
2552 {
2553  // check for piecewise functions
2554  checkForPiecewiseFunctions(dataModel, result);
2555 
2556  // check for initial assignments
2557  checkForInitialAssignments(dataModel, result);
2558 
2559  // check for event
2560  checkForEvents(dataModel, result);
2561 }
2562 
2563 /**
2564  * Checks whether the given data model can be exported to SBML Level2 Version1.
2565  * If it can be exported, the result vector will be empty, otherwise it will
2566  * contain a number of messages that specify why it can't be exported.
2567  */
2568 void CSBMLExporter::isModelSBMLL2V1Compatible(const CCopasiDataModel& dataModel, std::vector<SBMLIncompatibility>& result)
2569 {
2570  // check for initial assignments
2571  CSBMLExporter::checkForInitialAssignments(dataModel, result);
2572 }
2573 
2574 /**
2575  * Checks whether the given data model can be exported to SBML Level2 Version3.
2576  * If it can be exported, the result vector will be empty, otherwise it will
2577  * contain a number of messages that specify why it can't be exported.
2578  */
2579 void CSBMLExporter::isModelSBMLL2V3Compatible(const CCopasiDataModel& dataModel, std::vector<SBMLIncompatibility>& result)
2580 {
2581  // if there is an SBML model, which means the model was imported from SBML,
2582  // we have to check for spatial size units on species in the SBML model
2583  // because those are not allowed in SBMLLl2V3 and above
2584  check_for_spatial_size_units(dataModel, result);
2585 }
2586 
2587 /**
2588  * Go through all species in the model and check if the corresponding species
2589  * in the SBML model has the spatialSizeUnits attribute set.
2590  * This attribute is not supported in SBML L2V3 and above, so we have to get
2591  * rid of this attribute when we export to a level equal to or higher than
2592  * L2V3.
2593  * If the attribute has the same value as the compartments units, we can just
2594  * delete it without changing the model, otherwise we have to give a
2595  * corresponding warning.
2596  */
2597 void CSBMLExporter::check_for_spatial_size_units(const CCopasiDataModel& dataModel, std::vector<SBMLIncompatibility>& result)
2598 {
2599  const SBMLDocument* pSBMLDocument = const_cast<const SBMLDocument*>(const_cast<CCopasiDataModel&>(dataModel).getCurrentSBMLDocument());
2600 
2601  if (pSBMLDocument != NULL)
2602  {
2603  // check all species in the model if they have a spatial size attribute set
2604  // and if it is identical to the unit of the compartment the species is in
2605  const CModel* pModel = dataModel.getModel();
2606 
2607  if (pModel != NULL)
2608  {
2609  CCopasiVector<CMetab>::const_iterator it = pModel->getMetabolites().begin(), endit = pModel->getMetabolites().end();
2610  std::set<std::string> badSpecies;
2611  const std::map<CCopasiObject*, SBase*>& copasi2sbmlmap = const_cast<CCopasiDataModel&>(dataModel).getCopasi2SBMLMap();
2612  std::map<CCopasiObject*, SBase*>::const_iterator pos;
2613  const Species* pSBMLSpecies = NULL;
2614  std::string spatialSizeUnits;
2615 
2616  while (it != endit)
2617  {
2618  pos = copasi2sbmlmap.find(*it);
2619 
2620  if (pos != copasi2sbmlmap.end())
2621  {
2622  // check for the spatial size units attribute
2623  pSBMLSpecies = dynamic_cast<const Species*>(pos->second);
2624  assert(pSBMLSpecies != NULL);
2625 
2626  if (pSBMLSpecies == NULL) continue;
2627 
2628  if (pSBMLSpecies->isSetSpatialSizeUnits())
2629  {
2630  spatialSizeUnits = pSBMLSpecies->getSpatialSizeUnits();
2631  // check if the units are the same as the one on the species
2632  // compartment
2633  assert(pSBMLDocument->getModel() != NULL);
2634  const Compartment* pCompartment = pSBMLDocument->getModel()->getCompartment(pSBMLSpecies->getCompartment());
2635  assert(pCompartment != NULL);
2636 
2637  if (pCompartment != NULL)
2638  {
2639  UnitDefinition* pUDef1 = NULL;
2640  UnitDefinition* pUDef2 = NULL;
2641 
2642  if (pCompartment->isSetUnits())
2643  {
2644  assert(pSBMLDocument->getModel() != NULL);
2645 
2646  if (pSBMLDocument->getModel() != NULL)
2647  {
2648  pUDef1 = SBMLImporter::getSBMLUnitDefinitionForId(pCompartment->getUnits(), pSBMLDocument->getModel());
2649  }
2650  }
2651  else
2652  {
2653  // the compartment has the default units associated with the
2654  // symbol length , area or volume depending on the spatial size
2655  // of the compartment
2656  assert(pSBMLDocument->getModel() != NULL);
2657 
2658  if (pSBMLDocument->getModel() != NULL)
2659  {
2660  switch (pCompartment->getSpatialDimensions())
2661  {
2662  case 0:
2663  // the species is not allowed to have a
2664  // spatialDimensionsUnit attribute
2665  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 83 , pSBMLSpecies->getId().c_str());
2666  break;
2667 
2668  case 1:
2669  pUDef1 = SBMLImporter::getSBMLUnitDefinitionForId("length", pSBMLDocument->getModel());
2670  break;
2671 
2672  case 2:
2673  pUDef1 = SBMLImporter::getSBMLUnitDefinitionForId("area", pSBMLDocument->getModel());
2674  break;
2675 
2676  case 3:
2677  pUDef1 = SBMLImporter::getSBMLUnitDefinitionForId("volume", pSBMLDocument->getModel());
2678  break;
2679 
2680  default:
2681  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 82 , pCompartment->getId().c_str());
2682  break;
2683  }
2684  }
2685  }
2686 
2687  if (pUDef1 != NULL && pUDef2 != NULL)
2688  {
2689  // compare the two unit definitions
2690  if (!SBMLImporter::areSBMLUnitDefinitionsIdentical(pUDef1, pUDef2))
2691  {
2692  // add the species to bad species
2693  badSpecies.insert(pSBMLSpecies->getId());
2694  }
2695  }
2696 
2697  // delete the unit definitions
2698  if (pUDef1 != NULL)
2699  {
2700  delete pUDef1;
2701  }
2702 
2703  if (pUDef2 != NULL)
2704  {
2705  delete pUDef2;
2706  }
2707  }
2708  }
2709  }
2710 
2711  ++it;
2712  }
2713 
2714  if (!badSpecies.empty())
2715  {
2716  // create the incompatibility message
2717  std::ostringstream os;
2718  std::set<std::string>::const_iterator sit = badSpecies.begin(), sendit = badSpecies.end();
2719 
2720  while (sit != sendit)
2721  {
2722  os << *sit << ", ";
2723  ++sit;
2724  }
2725 
2726  result.push_back(SBMLIncompatibility(2, os.str().substr(0, os.str().size() - 2).c_str()));
2727  }
2728  }
2729  }
2730 }
2731 
2732 /**
2733  * Checks whether the model contains a metabolite that is defined by an ODE
2734  * expression and that is located in a variable volume. Since COPASI
2735  * interprets the expression differently from SBML, we can not correctly
2736  * export this yet. See Bug 903.
2737  */
2738 void CSBMLExporter::checkForODESpeciesInNonfixedCompartment(const CCopasiDataModel& dataModel, std::vector<SBMLIncompatibility> result)
2739 {
2740  const CModel* pModel = dataModel.getModel();
2741  const CCopasiVector<CMetab>& metabolites = pModel->getMetabolites();
2742  CCopasiVector<CMetab>::const_iterator it = metabolites.begin(), endit = metabolites.end();
2743 
2744  while (it != endit)
2745  {
2746  if ((*it)->getStatus() == CModelValue::ODE)
2747  {
2748  const CCompartment* pCompartment = (*it)->getCompartment();
2749  assert(pCompartment != NULL);
2750 
2751  if (pCompartment->getStatus() != CModelValue::FIXED)
2752  {
2753  result.push_back(SBMLIncompatibility(3, (*it)->getObjectName().c_str(), pCompartment->getObjectName().c_str()));
2754  }
2755  }
2756 
2757  ++it;
2758  }
2759 }
2760 
2761 /**
2762  * Checks whether the rule in the given model entity can be exported to
2763  * the specified version of SBML.
2764  * If it can be exported, the result vector will be empty, otherwise it will
2765  * contain a number of messages that specify why it can't be exported.
2766  */
2768  , const CCopasiDataModel& dataModel
2769  , int sbmlLevel
2770  , int sbmlVersion
2771  , std::vector<SBMLIncompatibility>& result
2772  , const std::string& objectDescription
2773  , bool initialExpression /* = false */
2774  , std::map<const std::string, Parameter*>* initialMap /* = NULL */)
2775 {
2776  checkForUnsupportedObjectReferences(expr, dataModel, sbmlLevel, sbmlVersion, result, initialExpression, initialMap);
2777  std::set<CEvaluationNodeFunction::SubType> unsupportedFunctionTypes = CSBMLExporter::createUnsupportedFunctionTypeSet(sbmlLevel);
2778  checkForUnsupportedFunctionCalls(*expr.getRoot(), unsupportedFunctionTypes, result, objectDescription);
2779 }
2780 
2781 /**
2782  * This static methods checks, whether the model uses any function calls
2783  * that can not be expressed in SBML like the random distribution
2784  * functions.
2785  */
2787  unsigned int sbmlLevel, unsigned int /*sbmlVersion*/, std::vector<SBMLIncompatibility>& result)
2788 {
2789  // Fill the set of unsupported functions depending on the level and
2790  // version
2791  std::set<CEvaluationNodeFunction::SubType> unsupportedFunctionTypes = CSBMLExporter::createUnsupportedFunctionTypeSet(sbmlLevel);
2792  // check all metabolites,parameters and compartments
2793  // make sure the list of assignments and initial assignments is filled
2794  // before this function is called
2795  size_t i, iMax = mAssignmentVector.size();
2796  const CModelEntity* pME = NULL;
2797 
2798  for (i = 0; i < iMax; ++i)
2799  {
2800  pME = mAssignmentVector[i];
2801  assert(pME != NULL);
2802 
2803  if (pME != NULL)
2804  {
2805  checkForUnsupportedFunctionCalls(*pME->getExpressionPtr()->getRoot(), unsupportedFunctionTypes, result, std::string(pME->getObjectType() + " called \"" + pME->getObjectName() + "\""));
2806  }
2807  }
2808 
2809  // check ode rules
2810  iMax = mODEVector.size();
2811 
2812  for (i = 0; i < iMax; ++i)
2813  {
2814  pME = mODEVector[i];
2815  assert(pME != NULL);
2816 
2817  if (pME != NULL)
2818  {
2819  checkForUnsupportedFunctionCalls(*pME->getExpressionPtr()->getRoot(), unsupportedFunctionTypes, result, std::string(pME->getObjectType() + " called \"" + pME->getObjectName() + "\""));
2820  }
2821  }
2822 
2823  // check initial assignments
2824  iMax = mInitialAssignmentVector.size();
2825 
2826  for (i = 0; i < iMax; ++i)
2827  {
2828  pME = mInitialAssignmentVector[i];
2829  assert(pME != NULL);
2830 
2831  if (pME != NULL)
2832  {
2833  checkForUnsupportedFunctionCalls(*pME->getInitialExpressionPtr()->getRoot(), unsupportedFunctionTypes, result, std::string(pME->getObjectType() + " called \"" + pME->getObjectName() + "\""));
2834  }
2835  }
2836 
2837  // if we already have a list of used functions, we can go through this here
2838  // make sure the list of used function has been filled before this is
2839  // called
2840  iMax = mUsedFunctions.size();
2841  std::set<CFunction*>::iterator it = this->mUsedFunctions.begin(), endit = this->mUsedFunctions.end();
2842 
2843  while (it != endit)
2844  {
2845  assert(*it != NULL);
2846 
2847  if (*it != NULL)
2848  {
2849  checkForUnsupportedFunctionCalls(*(*it)->getRoot(), unsupportedFunctionTypes, result, std::string("function called \"" + (*it)->getObjectName() + "\""));
2850  }
2851 
2852  ++it;
2853  }
2854 }
2855 
2856 /**
2857  * This static methods checks recursively, whether the given CEvaluationNode constains any function calls
2858  * that can not be expressed in SBML like the random distribution
2859  * functions.
2860  */
2862  const std::set<CEvaluationNodeFunction::SubType>& unsupportedFunctions,
2863  std::vector<SBMLIncompatibility>& result, const std::string& objectDescription)
2864 {
2866  {
2868 
2869  if (unsupportedFunctions.find(subtype) != unsupportedFunctions.end())
2870  {
2871  result.push_back(SBMLIncompatibility(2, node.getData().c_str(), objectDescription.c_str()));
2872  }
2873  }
2874 
2875  const CEvaluationNode* pChild = dynamic_cast<const CEvaluationNode*>(node.getChild());
2876 
2877  while (pChild != NULL)
2878  {
2879  checkForUnsupportedFunctionCalls(*pChild, unsupportedFunctions, result, objectDescription);
2880  pChild = dynamic_cast<const CEvaluationNode*>(pChild->getSibling());
2881  }
2882 }
2883 
2884 /**
2885  * This method checks whether the given model contains any initial assignments.
2886  */
2887 void CSBMLExporter::checkForInitialAssignments(const CCopasiDataModel& dataModel, std::vector<SBMLIncompatibility>& result)
2888 {
2889  const CModel* pModel = dataModel.getModel();
2890 
2891  if (pModel != NULL)
2892  {
2893  // check for rules
2894  const CCopasiVectorNS<CCompartment>& compartments = pModel->getCompartments();
2895  CCopasiVectorNS<CCompartment>::const_iterator compIt = compartments.begin(), compEndit = compartments.end();
2896 
2897  while (compIt != compEndit)
2898  {
2899  if ((*compIt)->getInitialExpression() != "")
2900  {
2901  result.push_back(SBMLIncompatibility(5, "Compartment", (*compIt)->getObjectName().c_str()));
2902  }
2903 
2904  ++compIt;
2905  }
2906 
2907  const CCopasiVector<CMetab>& metabs = pModel->getMetabolites();
2908 
2909  CCopasiVector<CMetab>::const_iterator metabIt = metabs.begin(), metabEndit = metabs.end();
2910 
2911  while (metabIt != metabEndit)
2912  {
2913  if ((*metabIt)->getInitialExpression() != "")
2914  {
2915  result.push_back(SBMLIncompatibility(5, "Metabolite", (*metabIt)->getObjectName().c_str()));
2916  }
2917 
2918  ++metabIt;
2919  }
2920 
2921  const CCopasiVectorN<CModelValue>& parameters = pModel->getModelValues();
2922 
2923  CCopasiVectorN<CModelValue>::const_iterator mvIt = parameters.begin(), mvEndit = parameters.end();
2924 
2925  while (mvIt != mvEndit)
2926  {
2927  if ((*mvIt)->getInitialExpression() != "")
2928  {
2929  result.push_back(SBMLIncompatibility(5, "Parameter", (*mvIt)->getObjectName().c_str()));
2930  }
2931 
2932  ++mvIt;
2933  }
2934  }
2935 }
2936 
2937 /**
2938  * Create all function definitions.
2939  */
2941 {
2942  this->mExportedFunctions.clear(true);
2943  this->mFunctionMap.clear();
2944  // make sure the list of used functions is filled before this is
2945  // called
2946  // make sure the mCOPASI2SBMLMap is up to date
2947 
2948  // find all indirectly called functions
2949  std::vector<CFunction*> usedFunctions = findUsedFunctions(this->mUsedFunctions, CCopasiRootContainer::getFunctionList());
2950  this->mUsedFunctions.clear();
2951  this->mUsedFunctions.insert(usedFunctions.begin(), usedFunctions.end());
2952 
2953  // remove the function calls that have been created by copasi solely for use
2954  // as a kinetic law term if it is no longer needed
2955  Model* pModel = this->mpSBMLDocument->getModel();
2956  assert(pModel);
2957 
2958  if (pModel == NULL) fatalError();
2959 
2960  std::map<SBase*, const CCopasiObject*> sbml2copasiMap;
2961  std::map<const CCopasiObject*, SBase*>::iterator mapIt = this->mCOPASI2SBMLMap.begin();
2962  std::map<const CCopasiObject*, SBase*>::iterator mapEndit = this->mCOPASI2SBMLMap.end();
2963 
2964  while (mapIt != mapEndit)
2965  {
2966  sbml2copasiMap.insert(std::make_pair(mapIt->second, mapIt->first));
2967  ++mapIt;
2968  }
2969 
2970  std::set<CFunction*> unusedFunctions;
2971  unsigned int i = 0, iMax = pModel->getNumFunctionDefinitions();
2972  std::map<SBase*, const CCopasiObject*>::iterator mapPos;
2973  std::set<std::string> toRemove;
2974 
2975  while (i < iMax)
2976  {
2977  FunctionDefinition* pFunDef = pModel->getFunctionDefinition(i);
2978 
2979  if (pFunDef != NULL)
2980  {
2981  mapPos = sbml2copasiMap.find(pFunDef);
2982 
2983  if (mapPos != sbml2copasiMap.end())
2984  {
2985  CFunction* pFun = dynamic_cast<CFunction*>(const_cast<CCopasiObject*>(mapPos->second));
2986 
2987  if (pFun != NULL && this->mUsedFunctions.find(pFun) == this->mUsedFunctions.end())
2988  {
2989  // the function exists in the model, but it is not used in any
2990  // expression
2991  if (pFun->getObjectName().find("function_4_") == 0 || pFun->getObjectName().find("Function for ") == 0)
2992  {
2993  // store the function definition that is to be removed
2994  toRemove.insert(pFunDef->getId());
2995  }
2996  else
2997  {
2998  // those need to be stored in a separate list since we also
2999  // need to find indirectly called functions
3000  unusedFunctions.insert(pFun);
3001  }
3002  }
3003  }
3004  }
3005 
3006  ++i;
3007  }
3008 
3009  // find all indirectly called functions for the unused functions
3010  std::vector<CFunction*> functionsVect = findUsedFunctions(unusedFunctions, CCopasiRootContainer::getFunctionList());
3011  usedFunctions.insert(usedFunctions.end(), functionsVect.begin(), functionsVect.end());
3012  // reset the used functions set
3013  this->mUsedFunctions.clear();
3014  this->mUsedFunctions.insert(usedFunctions.begin(), usedFunctions.end());
3015 
3016  // now we remove the function definitions from the SBML model
3017  std::set<std::string>::iterator toRemoveIt = toRemove.begin();
3018  std::set<std::string>::iterator toRemoveEndit = toRemove.end();
3019 
3020  while (toRemoveIt != toRemoveEndit)
3021  {
3022  FunctionDefinition* pFunDef = pModel->getFunctionDefinition(*toRemoveIt);
3023 
3024  if (pFunDef != NULL)
3025  {
3026  mapPos = sbml2copasiMap.find(pFunDef);
3027 
3028  if (mapPos != sbml2copasiMap.end() &&
3029  this->mUsedFunctions.find(dynamic_cast<CFunction*>(const_cast<CCopasiObject*>(mapPos->second))) == this->mUsedFunctions.end())
3030  {
3031  pModel->getListOfFunctionDefinitions()->remove(*toRemoveIt);
3032  }
3033  }
3034 
3035  ++toRemoveIt;
3036  }
3037 
3038  // order the remaining function definitions
3039  // remove duplicates from the vector, always keep the last one
3040  std::vector<CFunction*>::reverse_iterator reverseIt = usedFunctions.rbegin(), reverseEndit = usedFunctions.rend();
3041  functionsVect.clear();
3042 
3043  while (reverseIt != reverseEndit)
3044  {
3045  if (std::find(functionsVect.begin(), functionsVect.end(), *reverseIt) == functionsVect.end())
3046  {
3047  functionsVect.insert(functionsVect.begin(), *reverseIt);
3048  }
3049 
3050  ++reverseIt;
3051  }
3052 
3053  // remove all existing function definitions from the list
3054  // and remove their entries from the copasi2sbmlmap
3055  iMax = pModel->getNumFunctionDefinitions();
3056 
3057  for (unsigned int i = 0; i < iMax; ++i)
3058  {
3059  mapIt = this->mCOPASI2SBMLMap.begin();
3060  mapEndit = this->mCOPASI2SBMLMap.end();
3061 
3062  while (mapIt != mapEndit)
3063  {
3064  if (mapIt->second == pModel->getFunctionDefinition(i))
3065  {
3066  this->mCOPASI2SBMLMap.erase(mapIt->first);
3067  break;
3068  }
3069 
3070  ++mapIt;
3071  }
3072  }
3073 
3074  // don't want to remove the new custom function definitions if present
3075  // pModel->getListOfFunctionDefinitions()->clear(true);
3076  for (int i = ((int)pModel->getNumFunctionDefinitions()) - 1; i >= 0; --i)
3077  {
3078  FunctionDefinition* current = pModel->getFunctionDefinition(i);
3079 
3080  if (current != NULL &&
3081  current->isSetAnnotation() &&
3082  current->getAnnotation()->getNumChildren() == 1 &&
3083  current->getAnnotation()->getChild(0).getURI().find("http://sbml.org/annotations/") != std::string::npos)
3084  continue;
3085 
3086  // if we are here it is not one of the custom ones, so delete it
3087  FunctionDefinition* toDelete = pModel->removeFunctionDefinition(i);
3088  delete toDelete;
3089  }
3090 
3091  std::vector<CFunction*>::iterator it = functionsVect.begin(), endit = functionsVect.end();
3092 
3093  while (it != endit)
3094  {
3095  if (*it != NULL)
3096  {
3097  createFunctionDefinition(**it, dataModel);
3098  }
3099 
3100  ++it;
3101  }
3102 }
3103 
3104 /**
3105  * Create the SBML function definition from the given COPASI function.
3106  */
3108 {
3109  // check the expression
3110  std::map<const CCopasiObject*, SBase*>::iterator pos = this->mCOPASI2SBMLMap.find(&function);
3111  FunctionDefinition* pFunDef = NULL;
3112 
3113  // the entry could be NULL when exporting to L2 first and then to L3 because the
3114  // SBMLExporter::collectIds will add a NULL entry into the id table which will be
3115  // tranferred to the COPASI2SBMLMap.
3116  if (pos != this->mCOPASI2SBMLMap.end() && pos->second != NULL)
3117  {
3118  pFunDef = dynamic_cast<FunctionDefinition*>(pos->second);
3119  assert(pFunDef);
3120  this->mpSBMLDocument->getModel()->getListOfFunctionDefinitions()->appendAndOwn(pFunDef);
3121  }
3122  else
3123  {
3124  if (this->mpSBMLDocument->getLevel() == 1)
3125  {
3126  pFunDef = new FunctionDefinition(2, 4);
3127  }
3128  else
3129  {
3130  pFunDef = new FunctionDefinition(this->mpSBMLDocument->getLevel(), this->mpSBMLDocument->getVersion());
3131  }
3132 
3133  this->mFunctionMap[pFunDef] = &function;
3134  pFunDef->setName(function.getObjectName());
3135  std::string id = function.getSBMLId();
3136 
3137  if (id.empty())
3138  {
3139  id = function.getObjectName();
3140 
3141  if (CSBMLExporter::isValidSId(id))
3142  {
3143  if (this->mIdMap.find(id) != this->mIdMap.end())
3144  {
3145  id = CSBMLExporter::createUniqueId(this->mIdMap, id, true);
3146  }
3147  }
3148  else
3149  {
3150  id = CSBMLExporter::createUniqueId(this->mIdMap, function.getObjectName(), false);
3151  }
3152  }
3153 
3154  this->mIdMap.insert(std::pair<const std::string, const SBase*>(id, pFunDef));
3155  pFunDef->setId(id);
3156  function.setSBMLId(id);
3157  this->mCOPASI2SBMLMap[&function] = pFunDef;
3158  }
3159 
3160  if (function.getRoot() == NULL)
3161  {
3162  std::string errorMessage = std::string("Can not export function '");
3163  errorMessage += function.getObjectName();
3164  errorMessage += std::string("'. Function does not have a valid root node.");
3165  CCopasiMessage(CCopasiMessage::EXCEPTION, errorMessage.c_str());
3166  }
3167  else
3168  {
3169  std::vector<SBMLIncompatibility> result;
3170  CSBMLExporter::isExpressionSBMLCompatible(function, dataModel, this->mSBMLLevel, this->mSBMLVersion, result, std::string("function with name \"" + function.getObjectName() + "\"").c_str());
3171 
3172  if (result.empty())
3173  {
3174  ASTNode* pFunNode = this->convertToASTNode(function.getRoot(), dataModel);
3175  // go through the AST tree and replace all function call nodes with with a call to the sbml id
3176  ASTNode* pLambda = new ASTNode(AST_LAMBDA);
3177  // add the parameters to the function definition
3178  const CFunctionParameters& funParams = function.getVariables();
3179  size_t i, iMax = funParams.size();
3180  ASTNode* pParamNode = NULL;
3181 
3182  for (i = 0; i < iMax; ++i)
3183  {
3184  pParamNode = new ASTNode(AST_NAME);
3185  pParamNode->setName(funParams[i]->getObjectName().c_str());
3186  pLambda->addChild(pParamNode);
3187  }
3188 
3189  pLambda->addChild(pFunNode);
3190 
3191  if (pFunDef->setMath(pLambda) != LIBSBML_OPERATION_SUCCESS)
3192  {
3194  "The function '%s' uses an invalid expression '%s' that cannot be exported to SBML, the function definition is not exported."
3195  , function.getObjectName().c_str()
3196  , function.getInfix().c_str());
3197  }
3198 
3199  this->mExportedFunctions.appendAndOwn(pFunDef);
3200  delete pLambda;
3201  }
3202  else
3203  {
3204  this->mIncompatibilities.insert(this->mIncompatibilities.end(), result.begin(), result.end());
3205  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 28, function.getObjectName().c_str());
3206  }
3207  }
3208 
3209  if (pFunDef != NULL)
3210  {
3211  CSBMLExporter::setSBMLNotes(pFunDef, &function);
3212  }
3213 
3214  CSBMLExporter::updateMIRIAMAnnotation(&function, pFunDef, this->mMetaIdMap);
3215 }
3216 
3217 /*
3218  * Some parameters and initial assignments that are only needed for the export
3219  * are tagged so as not to be removed. This function removes that tag and is to
3220  * be called after the export has been completed.
3221  */
3222 void removeStickyTagFromElements(SBMLDocument *pSBMLDocument)
3223 {
3224  if (pSBMLDocument == NULL || pSBMLDocument->getModel() == NULL)
3225  return;
3226 
3227  // reset sticky parameters
3228  for (unsigned int i = 0; i < pSBMLDocument->getModel()->getNumParameters(); ++i)
3229  pSBMLDocument->getModel()->getParameter(i)->setUserData(NULL);
3230 
3231  // reset sticky initial assignments
3232  for (unsigned int i = 0; i < pSBMLDocument->getModel()->getNumInitialAssignments(); ++i)
3233  pSBMLDocument->getModel()->getInitialAssignment(i)->setUserData(NULL);
3234 }
3235 
3236 /**
3237  * Export the model to SBML.
3238  * The SBML model is returned as a string. In case of an error, an
3239  * empty string is returned.
3240  */
3241 const std::string CSBMLExporter::exportModelToString(CCopasiDataModel& dataModel, unsigned int sbmlLevel, unsigned int sbmlVersion)
3242 {
3243  this->mSBMLLevel = sbmlLevel;
3244  this->mSBMLVersion = sbmlVersion;
3245  mHandledSBMLObjects.clear();
3246  createSBMLDocument(dataModel);
3247 
3248  if (this->mpSBMLDocument && this->mpSBMLDocument->getModel())
3249  {
3250  LayoutModelPlugin* lmPlugin = static_cast<LayoutModelPlugin*>(this->mpSBMLDocument->getModel()->getPlugin("layout"));
3251 
3252  if (lmPlugin != NULL && sbmlLevel > 1)
3253  {
3254  // the layout can only be exported for Level 2 or Level 3
3255  dataModel.getListOfLayouts()->exportToSBML(lmPlugin->getListOfLayouts(),
3256  this->mCOPASI2SBMLMap, mIdMap, this->mpSBMLDocument->getLevel(), this->mpSBMLDocument->getVersion());
3257 
3258 #if LIBSBML_VERSION >= 50400
3259 
3260  // also ensure ther is one global render information object
3261  if (lmPlugin->getNumLayouts() > 0 && getNumDefaultStyles() > 0)
3262  {
3263  RenderListOfLayoutsPlugin* plugin = static_cast<RenderListOfLayoutsPlugin *>(lmPlugin->getListOfLayouts()->getPlugin("render"));
3264 
3265  if (plugin != NULL)
3266  {
3267  if (plugin->getNumGlobalRenderInformationObjects() == 0)
3268  {
3269  GlobalRenderInformation* info = plugin->createGlobalRenderInformation();
3270  getDefaultStyle(0)->toSBML(info, this->mpSBMLDocument->getLevel(), this->mpSBMLDocument->getVersion());
3271  }
3272  }
3273  }
3274 
3275 #endif // LIBSBML_VERSION >= 50400
3276  }
3277  }
3278 
3279 #ifdef USE_SBMLUNIT
3280 
3281  if (this->mpSBMLDocument != NULL)
3282  {
3283  CSBMLunitInterface uif(this->mpSBMLDocument->getModel(), true);
3284  uif.determineUnits();
3285 
3286  // check if there were conflicts
3287  if (uif.getStatistics().all[5] == 0)
3288  {
3289  // check if there are unresolved parameter units left
3290  if ((uif.getStatistics().local[0] != 0 || uif.getStatistics().global[0] != 0) && uif.getStatistics().numbers[0] != 0)
3291  {
3292  // try with heuristics
3293  CSBMLunitInterface uif2(this->mpSBMLDocument->getModel(), true);
3294  uif2.setAssumeDimensionlessOne(true);
3295  uif2.determineUnits();
3296 
3297  // check again if there have been conflicts
3298  if (uif2.getStatistics().all[5] == 0)
3299  {
3300  // done use result from uif2
3301  // there were no conflicts, so we write the units that could be determined
3302  uif2.writeBackToModel();
3303  std::cerr << "undetermined: " << uif2.getStatistics().global[0] + uif2.getStatistics().local[0] << std::endl;
3304  }
3305  else
3306  {
3307  // TODO create appropriate warning
3308  std::cerr << "Warning. " << uif2.getStatistics().all[5] << " conflicts found." << std::endl;
3309  }
3310  }
3311  else
3312  {
3313  // no conflicts, so we can write the parameters that could be determined to the model
3314  uif.writeBackToModel();
3315  std::cerr << "undetermined: " << uif.getStatistics().global[0] + uif.getStatistics().local[0] << std::endl;
3316  }
3317  }
3318  else
3319  {
3320  // TODO create a warning
3321  std::cerr << "Warning. " << uif.getStatistics().all[5] << " conflicts found." << std::endl;
3322  }
3323  }
3324 
3325 #endif // USE_SBMLUNIT
3326 
3327  // export the model to a string
3328  if (this->mpSBMLDocument == NULL) return std::string();
3329 
3331  SBMLWriter* writer = new SBMLWriter();
3332 
3333  writer->setProgramName("COPASI");
3334  writer->setProgramVersion(CVersion::VERSION.getVersion().c_str());
3335 
3336  char* d = writer->writeToString(this->mpSBMLDocument);
3337  std::string returnValue = d;
3338 
3339  if (d) free(d);
3340 
3341  pdelete(writer);
3342 
3343  // mark some temporary objects as no longer needed
3345 
3346  // the workaround is no longer necessary, as the libSBML > 5.11.4 will include
3347  // a converter that not only converts to l1v1 but also inlines the compartment
3348  // sizes and changes pow implementation as needed by Gepasi
3349 #if LIBSBML_VERSION <= 51104
3350 
3351  // actually most of the work seems to be done by
3352  // libsbml already, so the following method doesn't have to
3353  // do much and most of the code is obsolete, at least for now
3354  if (sbmlLevel == 1 && sbmlVersion == 1)
3355  {
3356  // if there is an exception in the routine, we hand it up to the calling routine
3357  convert_to_l1v1(returnValue);
3358  }
3359 
3360 #endif
3361 
3362  return returnValue;
3363 }
3364 
3365 /*
3366  * Adds all created parameters to the model and sets them as initial assignments
3367  * to the originating elements.
3368  */
3369 void addInitialAssignmentsToModel(SBMLDocument* doc
3370  , std::map<const std::string, Parameter*>& initialValueMap
3371  , const CCopasiDataModel &dataModel)
3372 {
3373  if (doc == NULL || doc->getModel() == NULL || initialValueMap.size() == 0)
3374  return;
3375 
3376  std::map<const std::string, Parameter*>::const_iterator it;
3377 
3378  for (it = initialValueMap.begin(); it != initialValueMap.end(); ++it)
3379  {
3380  Parameter* param = it->second;
3381  // add parameter to model
3382  int result = doc->getModel()->addParameter(param);
3383  doc->getModel()->getParameter(param->getId())->setUserData((void*)"1");
3384 
3385  const CCopasiObject *obj = static_cast<const CCopasiObject *>(dataModel.getObject(it->first));
3386  const std::string &sbmlId = (static_cast<const CModelEntity*>(obj->getObjectParent()))->getSBMLId();
3387 
3388  // create an initial assignment for the newly created initial quantity
3389  // to synchronize it with the original element.
3390  InitialAssignment* ia = doc->getModel()->createInitialAssignment();
3391  ia->setSymbol(it->second->getId());
3392  ia->setMath(SBML_parseFormula(sbmlId.c_str()));
3393  ia->setUserData((void*)"1");
3394 
3395  delete param;
3396  }
3397 }
3398 
3400 {
3401  // reset warnings for missing entries in modelhistory
3405 
3406  const SBMLDocument* pOldSBMLDocument = dataModel.getCurrentSBMLDocument();
3407  const CModel* pModel = dataModel.getModel();
3408  assert(pModel != NULL);
3409 
3410  if (pOldSBMLDocument == NULL)
3411  {
3412  this->mpSBMLDocument = new SBMLDocument(this->mSBMLLevel, this->mSBMLVersion);
3413  }
3414  else
3415  {
3416  this->mpSBMLDocument = dynamic_cast<SBMLDocument*>(pOldSBMLDocument->clone());
3417  }
3418 
3419  if (this->mpSBMLDocument == NULL) fatalError();
3420 
3421 #if LIBSBML_VERSION >= 50000
3422  const std::string uri = (this->mSBMLLevel < 3 ? LayoutExtension::getXmlnsL2() : LayoutExtension::getXmlnsL3V1V1());
3423  this->mpSBMLDocument->enablePackage(uri, "layout", true);
3424 
3425  if (this->mSBMLLevel > 2)
3426  this->mpSBMLDocument->setPackageRequired("layout", false);
3427 
3428  const std::string renderuri = (this->mSBMLLevel < 3 ? RenderExtension::getXmlnsL2() : RenderExtension::getXmlnsL3V1V1());
3429  this->mpSBMLDocument->enablePackage(renderuri, "render", true);
3430 
3431  if (this->mSBMLLevel > 2)
3432  this->mpSBMLDocument->setPackageRequired("render", false);
3433 
3434 #endif
3435 
3436  if (this->mpSBMLDocument->getModel() == NULL)
3437  {
3438  // there might still be COPASI object that have an SBML id, we have to
3439  // collect those as well
3440  CSBMLExporter::collectIds(dataModel, this->mIdMap);
3441  std::string id = pModel->getSBMLId();
3442 
3443  if (id.empty())
3444  {
3445  id = CSBMLExporter::createUniqueId(this->mIdMap, pModel->getObjectName(), false);
3446  this->mIdMap.insert(std::pair<const std::string, const SBase*>(id, this->mpSBMLDocument->getModel()));
3447  }
3448 
3449  this->mpSBMLDocument->createModel(id);
3450  }
3451  else
3452  {
3453  SBMLUtils::collectIds(this->mpSBMLDocument->getModel(), this->mIdMap, this->mMetaIdMap);
3454  }
3455 
3456  this->mFunctionIdMap.clear();
3457  // update the copasi2sbmlmap
3458  updateCOPASI2SBMLMap(dataModel);
3459  // update the comments on the model
3460  this->mpSBMLDocument->getModel()->setName(pModel->getObjectName());
3461 
3462  if (pModel != NULL && this->mpSBMLDocument->getModel() != NULL)
3463  {
3464  CSBMLExporter::setSBMLNotes(this->mpSBMLDocument->getModel(), pModel);
3465  }
3466 
3467  // update the MIRIAM annotation on the model
3468  CSBMLExporter::updateMIRIAMAnnotation(pModel, this->mpSBMLDocument->getModel(), this->mMetaIdMap);
3469 #if LIBSBML_VERSION >= 40100
3470 
3471  if (this->mSBMLLevel > 2)
3472  {
3473  // we have to remove the conversionFactor because on import we multiplied the stoichiometries
3474  // with this factor
3475  this->mpSBMLDocument->getModel()->unsetConversionFactor();
3476  }
3477 
3478 #endif // LIBSBML_VERSION
3479 
3480  // create units, compartments, species, parameters, reactions, initial
3481  // assignment, assignments, (event) and function definitions
3482  createUnits(dataModel);
3483  // try to find a parameter that represents Avogadros number
3484  findAvogadro(dataModel);
3485 
3486  // check if there is an assignment to a volume
3487  // in that case we has to export all species with the
3488  // hasOnlySubstanceUnits flag set to true
3489  this->mVariableVolumes = this->hasVolumeAssignment(dataModel);
3490 
3491  createCompartments(dataModel);
3492  createMetabolites(dataModel);
3493  createParameters(dataModel);
3494 
3495  // TODO here we should check the expressions but the function definitions
3496  // TODO have not been created yet, so we can't
3497  // TODO Maybe we don't have to do all checks here.
3498  // TODO it would probably be enough to check if the expressions are valid and
3499  // TODO to check the function calls later on when the functions are created.
3500 
3501  // only export initial assignments for Level 2 Version 2 and above
3502  if (this->mSBMLLevel != 1 && !(this->mSBMLLevel == 2 && this->mSBMLVersion == 1))
3503  {
3504  createInitialAssignments(dataModel);
3505  }
3506  else
3507  {
3509  }
3510 
3511  // since the flux of a reaction can be referenced in an assignment,
3512  // we have to make sure that asll reactions do have SBML Ids prior to creating the rules
3513  // and events
3514  assignSBMLIdsToReactions(dataModel.getModel());
3515 
3516  createRules(dataModel);
3517  createEvents(dataModel);
3518 
3519  // we have to export the reactions after all other entities that have
3520  // expressions since an expression could contain a local parameter
3521  // If it did, it would have to be converted to a global parameter and that
3522  // has to be taken into consideration when creating the reactions
3523  createReactions(dataModel);
3524  // find all used functions
3525  createFunctionDefinitions(dataModel);
3526 
3527  if (this->mSBMLLevel == 1)
3528  {
3529  // do SBML Level1 Voodoo
3530  // e.g. replace some of the unsupported nodes with workarounds
3531  convertToLevel1();
3532  // delete all function definitions since L1 does not have functions
3533  Model* pSBMLModel = this->mpSBMLDocument->getModel();
3534  assert(pSBMLModel != NULL);
3535  int i = pSBMLModel->getListOfFunctionDefinitions()->size();
3536 
3537  while (i > 0)
3538  {
3539  --i;
3540  // fix for bug 1086
3541  // remove the object from the COPASI2SBML Map
3542  std::map<const CCopasiObject*, SBase*>::iterator it = this->mCOPASI2SBMLMap.begin(), endit = this->mCOPASI2SBMLMap.end();
3543 
3544  while (it != endit)
3545  {
3546  if (it->second == pSBMLModel->getFunctionDefinition(i))
3547  {
3548  this->mCOPASI2SBMLMap.erase(it);
3549  break;
3550  }
3551 
3552  ++it;
3553  }
3554 
3555  delete pSBMLModel->getListOfFunctionDefinitions()->remove(i);
3556  }
3557  }
3558  else
3559  {
3560  // add all function definitions to the model
3561  Model* pSBMLModel = this->mpSBMLDocument->getModel();
3562  assert(pSBMLModel != NULL);
3563 
3564  if (pSBMLModel != NULL)
3565  {
3566  unsigned int i = 0, iMax = this->mExportedFunctions.size();
3567 #if LIBSBML_VERSION >= 40100
3568  int result;
3569 #endif // LIBSBML_VERSION >= 40100
3570  FunctionDefinition* pFunDef = NULL;
3571  std::map<const FunctionDefinition*, const CCopasiObject*>::const_iterator funPos;
3572 
3573  while (i < iMax)
3574  {
3575  pFunDef = this->mExportedFunctions.get(i);
3576  assert(pFunDef != NULL);
3577  // match the namespaces to those of the model
3578  pFunDef->setSBMLNamespaces(pSBMLModel->getSBMLNamespaces());
3579 // add methods only return a value starting with libsbml 4
3580 #if LIBSBML_VERSION >= 40100
3581  result = pSBMLModel->addFunctionDefinition(pFunDef);
3582  assert(result == LIBSBML_OPERATION_SUCCESS);
3583 #else
3584  pSBMLModel->addFunctionDefinition(pFunDef);
3585  assert(pSBMLModel->getFunctionDefinition(pFunDef->getId()) != NULL);
3586 #endif // LIBSBML_VERSION >= 40100
3587  // now we need to add the newly created FunctionDefinition to the copasi2sbml map
3588  funPos = this->mFunctionMap.find(pFunDef);
3589  assert(funPos != this->mFunctionMap.end());
3590  pFunDef = pSBMLModel->getFunctionDefinition(pFunDef->getId());
3591  assert(pFunDef != NULL);
3592  this->mCOPASI2SBMLMap[funPos->second] = pFunDef;
3593  ++i;
3594  }
3595 
3596  //the assert below no longer holds, as we export uniform / normal as functiondefinitions
3597  //assert(pSBMLModel->getNumFunctionDefinitions() == this->mExportedFunctions.size());
3598  }
3599  }
3600 
3601  // In case the document is just a clone of an old model, we have to set the level and version
3602  if (this->mpSBMLDocument->getLevel() != this->mSBMLLevel || this->mpSBMLDocument->getVersion() != this->mSBMLVersion)
3603  {
3604 #if LIBSBML_VERSION >= 50400
3605 
3606  SBMLNamespaces targetNs(mSBMLLevel, mSBMLVersion); // use reference, as the ns gets cloned
3607 
3608  if (mSBMLLevel == 1 && mSBMLVersion == 1)
3609  {
3610  // l1v1 export tailor made to suit GEPASI, with pow -> ^ notation and
3611  // inlining of compartment sizes in kinetics.
3612 
3613  ConversionProperties prop(&targetNs);
3614  prop.addOption("convertToL1V1", true,
3615  "convert the document to SBML Level 1 Version 1");
3616  prop.addOption("changePow", true,
3617  "change pow expressions to the (^) hat notation");
3618  prop.addOption("inlineCompartmentSizes", true,
3619  "if true, occurrances of compartment ids in expressions will be replaced with their initial size");
3620 
3621  mpSBMLDocument->convert(prop);
3622  }
3623  else
3624  {
3625 
3626  ConversionProperties prop(&targetNs);
3627  prop.addOption("strict", false);
3628  prop.addOption("setLevelAndVersion", true);
3629  prop.addOption("ignorePackages", true);
3630 
3631  mpSBMLDocument->convert(prop);
3632  }
3633 
3634 #else
3635  this->mpSBMLDocument->setLevelAndVersion(this->mSBMLLevel, this->mSBMLVersion);
3636 #endif
3637  }
3638 
3639  if (this->mpSBMLDocument->getLevel() != this->mSBMLLevel || this->mpSBMLDocument->getVersion() != this->mSBMLVersion)
3640  {
3641  unsigned int i, iMax = this->mpSBMLDocument->getNumErrors();
3642  std::string message("libSBML could not convert the model to the requested level and/or version. The reason(s) given were:");
3643 
3644  for (i = 0; i < iMax ; ++i)
3645  {
3646  const XMLError* pSBMLError = this->mpSBMLDocument->getError(i);
3647  message += "\n";
3648  message += pSBMLError->getMessage();
3649  }
3650 
3651  CCopasiMessage(CCopasiMessage::EXCEPTION, message.c_str());
3652  }
3653 
3654  // remove mpAvogadro from the model again
3655  if (this->mAvogadroCreated == true)
3656  {
3657  std::map<const CCopasiObject*, SBase*>::iterator pos = this->mCOPASI2SBMLMap.find(this->mpAvogadro);
3658  this->mCOPASI2SBMLMap.erase(pos);
3659  dataModel.getModel()->removeModelValue(this->mpAvogadro->getKey(), true);
3660  }
3661 
3662  this->outputIncompatibilities();
3663 
3664  // if the model is incompatible with SBML export, throw an exception
3665  if (!this->mIncompatibilities.empty() && !this->mIncompleteExport)
3666  {
3667  CCopasiMessage(CCopasiMessage::EXCEPTION, "Model incompatible with chosen version and/or level of SBML.");
3668  }
3669 
3670  // if initial assignments were used we will have to add them to the document
3672  mInitialValueMap.clear();
3673 
3674  // delete all temporary function definitions
3675  // and the function map
3676  this->mExportedFunctions.clear(true);
3677  this->mFunctionMap.clear();
3678 }
3679 
3680 const std::set<CEvaluationNodeFunction::SubType> CSBMLExporter::createUnsupportedFunctionTypeSet(unsigned int sbmlLevel)
3681 {
3682  std::set<CEvaluationNodeFunction::SubType> unsupportedFunctionTypes;
3683 
3684  if (sbmlLevel == 1)
3685  {
3686  // check what to do with all the functions that are not supported by
3687  // SBML Level 1
3688  // Also take care of other unsupported nodes in L1
3689  // (inf,nan,pi,exponentiale)
3690  // most functions that are not supported by Level 1, but Level 2
3691  // can not be converted to something that can be exported.
3692  // Roundtripping of those will not succeed, but the results should stay
3693  // the same. The only function that can not be correctly converted
3694  // is ARCCOTH since we would need a piecewise for that which is not
3695  // supported in Level 1
3696  unsupportedFunctionTypes.insert(CEvaluationNodeFunction::ARCCOTH);
3697  //unsupportedFunctionTypes.insert(CEvaluationNodeFunction::RNORMAL);
3698  //unsupportedFunctionTypes.insert(CEvaluationNodeFunction::RUNIFORM);
3699  unsupportedFunctionTypes.insert(CEvaluationNodeFunction::MAX);
3700  unsupportedFunctionTypes.insert(CEvaluationNodeFunction::MIN);
3701  }
3702  else
3703  {
3704  //unsupportedFunctionTypes.insert(CEvaluationNodeFunction::RNORMAL);
3705  //unsupportedFunctionTypes.insert(CEvaluationNodeFunction::RUNIFORM);
3706  //unsupportedFunctionTypes.insert(CEvaluationNodeFunction::MAX);
3707  //unsupportedFunctionTypes.insert(CEvaluationNodeFunction::MIN);
3708  }
3709 
3710  return unsupportedFunctionTypes;
3711 }
3712 
3713 /**
3714  * Export the model to SBML.
3715  * The model is written to the file given by filename.
3716  * If the export fails, false is returned.
3717  */
3718 bool CSBMLExporter::exportModel(CCopasiDataModel& dataModel, const std::string& filename, unsigned int sbmlLevel, unsigned int sbmlVersion, bool overwrite)
3719 {
3720  bool success = true;
3721  /* create a string that represents the SBMLDocument */
3722  std::string str = this->exportModelToString(dataModel, sbmlLevel, sbmlVersion);
3723 
3724  if (!str.empty())
3725  {
3726  /* check if the file already exists.
3727  If yes, write if overwrite is true,
3728  else create an appropriate CCopasiMessage. */
3729  std::ifstream testInfile(CLocaleString::fromUtf8(filename).c_str(), std::ios::in);
3730 
3731  if (testInfile && !overwrite)
3732  {
3733  // create a CCopasiMessage with the appropriate error
3734  CCopasiMessage(CCopasiMessage::ERROR, MCDirEntry + 1, filename.c_str());
3735  return false;
3736  }
3737 
3738  /* write the document to a file */
3739  std::ofstream outfile(CLocaleString::fromUtf8(filename).c_str(), std::ios::out | std::ios::trunc);
3740  outfile << str;
3741  outfile.close();
3742  }
3743  else
3744  {
3745  /* if no SBMLDocument could be created return false */
3746  success = false;
3747  }
3748 
3749  return success;
3750 }
3751 
3752 /**
3753  * Checks whether the given data model can be exported to a certain version of SBML.
3754  * If it can be exported, the result vector will be empty, otherwise it will
3755  * contain a number of messages that specify why it can't be exported.
3756  */
3757 const std::vector<SBMLIncompatibility> CSBMLExporter::isModelSBMLCompatible(const CCopasiDataModel& dataModel, int sbmlLevel, int sbmlVersion)
3758 {
3759  const CModel* pModel = dataModel.getModel();
3760  std::vector<SBMLIncompatibility> result;
3761 
3762  if (pModel == NULL) return result;
3763 
3764  // general checks
3765  // check if there is a species with an ode rule that is in a nonfixed
3766  // compartment
3767  checkForODESpeciesInNonfixedCompartment(dataModel, result);
3768 
3769  // check if the model contains references to model entities that can not be
3770  // represented in SBML like the initial value of something as opposed to the
3771  // transient value
3772  // check if the model contains calls to functions that are not supported in
3773  // the given version of SBML
3774  CModelEntity::Status status;
3775  const CExpression* pExpression = NULL;
3776  std::set<std::string> usedFunctionNames;
3779 
3780  while (compIt != compEndit)
3781  {
3782  status = (*compIt)->getStatus();
3783 
3784  if (status == CModelEntity::ODE || status == CModelEntity::ASSIGNMENT)
3785  {
3786  pExpression = (*compIt)->getExpressionPtr();
3787  assert(pExpression != NULL);
3788 
3789  if (pExpression != NULL)
3790  {
3791  // check for unsupported object references and unsupported function
3792  // calls
3793  usedFunctionNames.clear();
3794  CSBMLExporter::isExpressionSBMLCompatible(*pExpression, dataModel, sbmlLevel, sbmlVersion, result, std::string("compartment with name \"" + (*compIt)->getObjectName() + "\"").c_str());
3795  CSBMLExporter::findDirectlyUsedFunctions(pExpression->getRoot(), usedFunctionNames);
3796  }
3797  }
3798 
3799  pExpression = (*compIt)->getInitialExpressionPtr();
3800 
3801  if (pExpression != NULL)
3802  {
3803  // check for unsupported object references and unsupported function
3804  // calls
3805  CSBMLExporter::isExpressionSBMLCompatible(*pExpression, dataModel, sbmlLevel, sbmlVersion, result, std::string("initial expression for compartment named \"" + (*compIt)->getObjectName() + "\"").c_str());
3806  usedFunctionNames.clear();
3807  CSBMLExporter::findDirectlyUsedFunctions(pExpression->getRoot(), usedFunctionNames);
3808  }
3809 
3810  ++compIt;
3811  }
3812 
3814  CCopasiVector<CMetab>::const_iterator metabEndit = pModel->getMetabolites().end();
3815 
3816  while (metabIt != metabEndit)
3817  {
3818  status = (*metabIt)->getStatus();
3819 
3820  if (status == CModelEntity::ODE || status == CModelEntity::ASSIGNMENT)
3821  {
3822  pExpression = (*metabIt)->getExpressionPtr();
3823  assert(pExpression != NULL);
3824 
3825  if (pExpression != NULL)
3826  {
3827  // check for unsupported object references and unsupported function
3828  // calls
3829  CSBMLExporter::isExpressionSBMLCompatible(*pExpression, dataModel, sbmlLevel, sbmlVersion, result, std::string("rule for species named \"" + (*metabIt)->getObjectName() + "\"").c_str());
3830  usedFunctionNames.clear();
3831  CSBMLExporter::findDirectlyUsedFunctions(pExpression->getRoot(), usedFunctionNames);
3832  }
3833  }
3834 
3835  pExpression = (*metabIt)->getInitialExpressionPtr();
3836 
3837  if (pExpression != NULL)
3838  {
3839  // check for unsupported object references and unsupported function
3840  // calls
3841  CSBMLExporter::isExpressionSBMLCompatible(*pExpression, dataModel, sbmlLevel, sbmlVersion, result, std::string("initial species for metabolite named \"" + (*metabIt)->getObjectName() + "\"").c_str());
3842  usedFunctionNames.clear();
3843  CSBMLExporter::findDirectlyUsedFunctions(pExpression->getRoot(), usedFunctionNames);
3844  }
3845 
3846  ++metabIt;
3847  }
3848 
3851 
3852  while (mvIt != mvEndit)
3853  {
3854  status = (*mvIt)->getStatus();
3855 
3856  if (status == CModelEntity::ODE || status == CModelEntity::ASSIGNMENT)
3857  {
3858  pExpression = (*mvIt)->getExpressionPtr();
3859  assert(pExpression != NULL);
3860 
3861  if (pExpression != NULL)
3862  {
3863  // check for unsupported object references and unsupported function
3864  // calls
3865  CSBMLExporter::isExpressionSBMLCompatible(*pExpression, dataModel, sbmlLevel, sbmlVersion, result, std::string("rule for global parameter named \"" + (*mvIt)->getObjectName() + "\"").c_str());
3866  usedFunctionNames.clear();
3867  CSBMLExporter::findDirectlyUsedFunctions(pExpression->getRoot(), usedFunctionNames);
3868  }
3869  }
3870 
3871  pExpression = (*mvIt)->getInitialExpressionPtr();
3872 
3873  if (pExpression != NULL)
3874  {
3875  // check for unsupported object references and unsupported function
3876  // calls
3877  CSBMLExporter::isExpressionSBMLCompatible(*pExpression, dataModel, sbmlLevel, sbmlVersion, result, std::string("initial expression for global parameter named \"" + (*mvIt)->getObjectName() + "\"").c_str());
3878  usedFunctionNames.clear();
3879  CSBMLExporter::findDirectlyUsedFunctions(pExpression->getRoot(), usedFunctionNames);
3880  }
3881 
3882  ++mvIt;
3883  }
3884 
3885  // since all kinetic laws in COPASI are simple function calls and the
3886  // arguments to the call can only be parameters or species, the expressions
3887  // in kinetic laws are always automatically valid for SBML export
3888  // However the functions that are called must not be valid
3889 
3890  // check all functions that are used if they contain invalid function calls
3891 
3892  // level dependent checks
3893  CCopasiVectorN<CEvent>::const_iterator eventIt, eventEndit;
3894 
3895  switch (sbmlLevel)
3896  {
3897  case 2:
3898  // check all events
3899  eventIt = pModel->getEvents().begin();
3900  eventEndit = pModel->getEvents().end();
3901 
3902  while (eventIt != eventEndit)
3903  {
3904  CSBMLExporter::isEventSBMLCompatible(*eventIt, dataModel, sbmlLevel, sbmlVersion, result);
3905  ++eventIt;
3906  }
3907 
3908  switch (sbmlVersion)
3909  {
3910  case 1:
3911  CSBMLExporter::isModelSBMLL2V1Compatible(dataModel, result);
3912  break;
3913 
3914  case 2:
3915  case 3:
3916  CSBMLExporter::isModelSBMLL2V3Compatible(dataModel, result);
3917  break;
3918 
3919  default:
3920  result.push_back(SBMLIncompatibility(6, sbmlLevel, sbmlVersion));
3921  break;
3922  }
3923 
3924  break;
3925 
3926  case 1:
3927  CSBMLExporter::isModelSBMLL1Compatible(dataModel, result);
3928  break;
3929 
3930  default:
3931  result.push_back(SBMLIncompatibility(6, sbmlLevel, sbmlVersion));
3932  break;
3933  }
3934 
3935  return result;
3936 }
3937 
3939 {
3940  if (this->mSBMLLevel == 1)
3941  {
3943  return;
3944  }
3945 
3946  // bail early
3947  if (dataModel.getModel() == NULL)
3948  return;
3949 
3950  Model* pSBMLModel = this->mpSBMLDocument->getModel();
3951 
3952  if (pSBMLModel == NULL)
3953  return;
3954 
3955  // remove all events from the SBML model and put them in a set
3956  std::set<Event*> eventSet;
3957 
3958  while (pSBMLModel->getNumEvents() != 0)
3959  {
3960  Event* pEvent = pSBMLModel->getListOfEvents()->remove(pSBMLModel->getNumEvents() - 1);
3961  assert(pEvent != NULL);
3962  eventSet.insert(pEvent);
3963  }
3964 
3965  const CCopasiVectorN<CEvent>& events = dataModel.getModel()->getEvents();
3966 
3967  CCopasiVectorN<CEvent>::const_iterator it = events.begin(), endit = events.end();
3968 
3969  Event* pSBMLEvent = NULL;
3970 
3971  std::map<const CCopasiObject*, SBase*>::const_iterator pos;
3972 
3973  while (it != endit)
3974  {
3975  // find the old event if there is one
3976  pos = this->mCOPASI2SBMLMap.find(*it);
3977 
3978  if (pos != this->mCOPASI2SBMLMap.end())
3979  {
3980  pSBMLEvent = dynamic_cast<Event*>(pos->second);
3981  }
3982  else
3983  {
3984  pSBMLEvent = NULL;
3985  }
3986 
3987  // remove the event from the set
3988  if (pSBMLEvent) eventSet.erase(pSBMLEvent);
3989 
3990  this->createEvent(**it, pSBMLEvent, dataModel);
3991  ++it;
3992  }
3993 
3994  // make sure events that have already been deleted do not get
3995  // exported
3996  std::set<Event*>::iterator it2 = eventSet.begin(), endit2 = eventSet.end();
3997 
3998  while (it2 != endit2)
3999  {
4000  delete *it2;
4001  ++it2;
4002  }
4003 }
4004 
4005 void CSBMLExporter::createEvent(CEvent& event, Event* pSBMLEvent, CCopasiDataModel& dataModel)
4006 {
4007  // once events are functional, we create them here
4008  // don't forget to call replaceSpeciesReferences
4009  // on all mathematical expressions
4010 
4011  // first check if we already have a corresponding event from an earlier
4012  // import or export
4013  // if none is found, we create a new event object
4014  if (pSBMLEvent == NULL)
4015  {
4016  pSBMLEvent = this->mpSBMLDocument->getModel()->createEvent();
4017  }
4018  else
4019  {
4020  // Read the event to the model
4021  if (this->mpSBMLDocument->getModel()->getListOfEvents()->appendAndOwn(pSBMLEvent) != LIBSBML_OPERATION_SUCCESS)
4022  {
4023  // the existing event is not valid and cannot be added to the model, delete and create a new one instead
4024  delete pSBMLEvent;
4025  pSBMLEvent = this->mpSBMLDocument->getModel()->createEvent();
4026  }
4027  }
4028 
4029  // add the object to the COPASI2SBMLMap
4030  this->mCOPASI2SBMLMap[&event] = pSBMLEvent;
4031 
4032  // a) we set the name and the id
4033  if (!pSBMLEvent->isSetId())
4034  {
4035  std::string sbmlId = CSBMLExporter::createUniqueId(this->mIdMap, event.getObjectName(), false);
4036  this->mIdMap.insert(std::pair<const std::string, const SBase*>(sbmlId, pSBMLEvent));
4037  pSBMLEvent->setId(sbmlId);
4038  event.setSBMLId(sbmlId);
4039  }
4040 
4041  // The attribute 'useValuesfromTriggerTime' is mandatory in L3V1
4042  if (this->mSBMLLevel > 2)
4043  {
4044  pSBMLEvent->setUseValuesFromTriggerTime(event.getDelayAssignment());
4045  }
4046 
4047 #if LIBSBML_VERSION >= 40200
4048 
4049  // if the event came from an SBML L3 model, we have to make sure that the
4050  // priority is not set when we export the event, otherwise the exported model
4051  // and the COPASI model will not behave the same.
4052  if (this->mSBMLLevel > 2)
4053  {
4054  pSBMLEvent->setPriority(NULL);
4055  }
4056 
4057 #endif // LIBSBML_VERSION >= 40200
4058 
4059  // if the name doesn't consists of whitespace only, we set the name
4060  if (event.getObjectName().find_first_not_of("\t\r\n ") != std::string::npos)
4061  {
4062  pSBMLEvent->setName(event.getObjectName());
4063  }
4064 
4065  // b) we check the trigger expression and export the trigger
4066  const CExpression* pExpression = event.getTriggerExpressionPtr();
4067  assert(pExpression != NULL);
4068  std::vector<SBMLIncompatibility> result;
4070  , dataModel
4071  , this->mSBMLLevel
4072  , this->mSBMLVersion
4073  , result
4074  , std::string("event trigger for event with id\"+" + event.getSBMLId() + "\"")
4075  , false
4076  , &mInitialValueMap);
4077 
4078  // collect directly used functions
4079  if (result.empty())
4080  {
4081  std::set<std::string> directlyUsedFunctionNames;
4082  CSBMLExporter::findDirectlyUsedFunctions(pExpression->getRoot(), directlyUsedFunctionNames);
4083  std::set<CFunction*> usedFunctions = CSBMLExporter::createFunctionSetFromFunctionNames(directlyUsedFunctionNames, CCopasiRootContainer::getFunctionList());
4084 
4085 # if defined _MSC_VER && _MSC_VER < 1201 // 1200 Identifies Visual C++ 6.0
4086  {
4087  std::set<CFunction*>::const_iterator it = usedFunctions.begin();
4088  std::set<CFunction*>::const_iterator end = usedFunctions.end();
4089 
4090  for (; it != end; ++it)
4091  this->mUsedFunctions.insert(*it);
4092  }
4093 #else
4094  this->mUsedFunctions.insert(usedFunctions.begin(), usedFunctions.end());
4095 #endif
4096  }
4097  else
4098  {
4099  this->mIncompatibilities.insert(this->mIncompatibilities.end(), result.begin(), result.end());
4100 
4101  if (!this->mIncompleteExport)
4102  {
4103  this->outputIncompatibilities();
4104  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 60, "trigger expression", "event", event.getObjectName().c_str());
4105  }
4106  }
4107 
4108  const std::string& changedExpression = convertExpression(pExpression->getInfix(), mInitialValueMap);
4109 
4110  CEvaluationTree tree;
4111 
4112  tree.setInfix(changedExpression);
4113 
4114  const CEvaluationNode* pOrigNode = tree.getRoot();
4115 
4116  //const CEvaluationNode* pOrigNode = pExpression->getRoot();
4117 
4119  {
4120  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 70, "trigger", "event", event.getObjectName().c_str());
4121  }
4122 
4123  // the next few lines replace references to species depending on whether
4124  // it is a reference to an amount or a reference to a concentration.
4125  // Other factors that influence this replacement are if the model
4126  // contains variable volumes or if the quantity units are set to CModel::number
4127  pOrigNode = this->replaceSpeciesReferences(pOrigNode, dataModel);
4128  assert(pOrigNode != NULL);
4129  ASTNode* pNode = this->convertToASTNode(pOrigNode, dataModel);
4130  delete pOrigNode;
4131  // convert local parameters that are referenced in an expression to global
4132  // parameters
4133  this->replace_local_parameters(pNode, dataModel);
4134 
4135  if (pNode != NULL)
4136  {
4137  Trigger* pTrigger = pSBMLEvent->createTrigger();
4138  pTrigger->setMath(pNode);
4139 #if LIBSBML_VERSION >= 40200
4140 
4141  // we need to make sure that the initial value of the trigger is set to true
4142  // as well as the persistent flag
4143  // for L3 event because otherwise the exported model will not do the same as
4144  // the model in COPASI
4145  if (this->mSBMLLevel > 2)
4146  {
4147  pTrigger->setInitialValue(true);
4148  pTrigger->setPersistent(true);
4149  }
4150 
4151 #endif // LIBSBML_VERSION >= 40200
4152  delete pNode;
4153  }
4154  else
4155  {
4156  if (this->mIncompleteExport != true)
4157  {
4158  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 60, "trigger", "event", event.getObjectName().c_str());
4159  }
4160  }
4161 
4162  // c) we check the delay expression and export it if it is set
4163  pExpression = event.getDelayExpressionPtr();
4164  const std::string delayInfix = pExpression == NULL ? "" : pExpression->getInfix();
4165 
4166  if (!delayInfix.empty())
4167  {
4168  result.clear();
4170  , dataModel
4171  , this->mSBMLLevel
4172  , this->mSBMLVersion
4173  , result
4174  , std::string("event delay for event with id\"+" + event.getSBMLId() + "\"").c_str()
4175  , false
4176  , &mInitialValueMap);
4177 
4178  // collect directly used functions
4179  if (result.empty())
4180  {
4181  std::set<std::string> directlyUsedFunctionNames;
4182  CSBMLExporter::findDirectlyUsedFunctions(pExpression->getRoot(), directlyUsedFunctionNames);
4183  std::set<CFunction*> usedFunctions = CSBMLExporter::createFunctionSetFromFunctionNames(directlyUsedFunctionNames, CCopasiRootContainer::getFunctionList());
4184 
4185 # if defined _MSC_VER && _MSC_VER < 1201 // 1200 Identifies Visual C++ 6.0
4186  {
4187  it = usedFunctions.begin();
4188  usedFunctions.end();
4189 
4190  for (; it != end; ++it)
4191  this->mUsedFunctions.insert(*it);
4192  }
4193 #else
4194  this->mUsedFunctions.insert(usedFunctions.begin(), usedFunctions.end());
4195 #endif
4196  }
4197  else
4198  {
4199  this->mIncompatibilities.insert(this->mIncompatibilities.end(), result.begin(), result.end());
4200 
4201  if (!this->mIncompleteExport)
4202  {
4203  this->outputIncompatibilities();
4204  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 60, "delay expression", "event", event.getObjectName().c_str());
4205  }
4206  }
4207 
4208  //pOrigNode = pExpression->getRoot();
4209  const std::string& changedExpression = convertExpression(pExpression->getInfix(), mInitialValueMap);
4210  CEvaluationTree tree;
4211  tree.setInfix(changedExpression);
4212  const CEvaluationNode* pOrigNode = tree.getRoot();
4213 
4215  {
4216  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 70, "delay", "event", event.getObjectName().c_str());
4217  }
4218 
4219  // the next few lines replace references to species depending on whether
4220  // it is a reference to an amount or a reference to a concentration.
4221  // Other factors that influence this replacement are if the model
4222  // contains variable volumes or if the quantity units are set to CModel::number
4223  pOrigNode = this->replaceSpeciesReferences(pOrigNode, dataModel);
4224  assert(pOrigNode != NULL);
4225  pNode = this->convertToASTNode(pOrigNode, dataModel);
4226  delete pOrigNode;
4227  // convert local parameters that are referenced in an expression to global
4228  // parameters
4229  this->replace_local_parameters(pNode, dataModel);
4230 
4231  if (pNode != NULL)
4232  {
4233  Delay* pDelay = new Delay(this->mSBMLLevel, this->mSBMLVersion);
4234  pDelay->setMath(pNode);
4235  pSBMLEvent->setDelay(pDelay);
4236  delete pNode;
4237  delete pDelay;
4238 
4239  if (this->mSBMLLevel == 2 && this->mSBMLVersion >= 4)
4240  {
4241  pSBMLEvent->setUseValuesFromTriggerTime(event.getDelayAssignment());
4242  }
4243  // if the delayAssignment is set to false, we have a problem
4244  // because the model semantic changes on export
4245  else if (event.getDelayAssignment() == false &&
4246  this->mSBMLLevel <= 2 &&
4247  this->mSBMLVersion <= 3)
4248  {
4249  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 77, event.getObjectName().c_str(), this->mSBMLLevel, this->mSBMLVersion);
4250  }
4251  }
4252  else
4253  {
4254  if (this->mIncompleteExport != true)
4255  {
4256  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 60, "delay", "event", event.getObjectName().c_str());
4257  }
4258  }
4259  }
4260 
4261  // export all event assignments
4262  this->exportEventAssignments(event, pSBMLEvent, dataModel);
4263 
4264  // check if the event has event assignments, if not, we have to delete
4265  // the event and create an error message (for SBML Levels < 3)
4266  if (pSBMLEvent->getNumEventAssignments() == 0 && mpSBMLDocument->getLevel() < 3)
4267  {
4268  unsigned int i, iMax = this->mpSBMLDocument->getModel()->getNumEvents();
4269 
4270  // the following code does not check if the corrct pointer is found
4271  // but this should not be a problem
4272  for (i = 0; i < iMax; ++i)
4273  {
4274  if (this->mpSBMLDocument->getModel()->getEvent(i) == pSBMLEvent)
4275  {
4276  this->mpSBMLDocument->getModel()->getListOfEvents()->remove(i);
4277  }
4278  }
4279 
4280  delete pSBMLEvent;
4281  this->mCOPASI2SBMLMap.erase(&event);
4282  pSBMLEvent = NULL;
4283 
4284  CCopasiMessage(CCopasiMessage::WARNING, "The event '%s' contained no event assignments, this is only supported for SBML Level 3 onwards. The event has not been exported.", event.getObjectName().c_str());
4285  }
4286 
4287  if (pSBMLEvent != NULL)
4288  {
4289  CSBMLExporter::setSBMLNotes(pSBMLEvent, &event);
4290  CSBMLExporter::updateMIRIAMAnnotation(&event, pSBMLEvent, this->mMetaIdMap);
4291  }
4292 }
4293 
4294 void CSBMLExporter::exportEventAssignments(const CEvent& event, Event* pSBMLEvent, CCopasiDataModel& dataModel)
4295 {
4296  // we export all event assignments
4297  // i) delete all event assignments and add them to a vector (see rules)
4298  std::map<std::string, EventAssignment*> assignmentMap;
4299 
4300  while (pSBMLEvent->getNumEventAssignments() != 0)
4301  {
4302  EventAssignment* pAssignment = pSBMLEvent->getEventAssignment(pSBMLEvent->getNumEventAssignments() - 1);
4303  assert(pAssignment != NULL);
4304  assignmentMap[pAssignment->getVariable()] = pAssignment;
4305  pSBMLEvent->getListOfEventAssignments()->remove(pSBMLEvent->getNumEventAssignments() - 1);
4306  }
4307 
4308  CCopasiVectorN< CEventAssignment >::const_iterator itAssignment = event.getAssignments().begin();
4309  CCopasiVectorN< CEventAssignment >::const_iterator endAssignment = event.getAssignments().end();
4310 
4311  for (; itAssignment != endAssignment; ++itAssignment)
4312  {
4313  // ii) for each assignment, check if there already is an old assignment
4314  // for this variable
4315  std::string key = (*itAssignment)->getTargetKey();
4316 
4317  // now we have to get the object for the key, check if the object is a
4318  // compartment, species or global parameter,
4319  const CCopasiObject* pObject = CCopasiRootContainer::getKeyFactory()->get(key);
4320  std::string objectType = pObject->getObjectType();
4321 
4322  if (objectType == "Reference")
4323  {
4324  pObject = pObject->getObjectParent();
4325  objectType = pObject->getObjectType();
4326  }
4327 
4328  if (objectType != "Compartment" && objectType != "Metabolite" && objectType != "ModelValue")
4329  {
4330  // create an error message that only compartment, species and
4331  // global parameters can be used in an event assignment
4332  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 32, std::string("event assignment for event called \"" + event.getObjectName() + "\"").c_str(), pObject->getObjectName().c_str());
4333  }
4334 
4335  std::map<const CCopasiObject*, SBase*>::const_iterator pos = this->mCOPASI2SBMLMap.find(pObject);
4336 
4337  if (pos == this->mCOPASI2SBMLMap.end())
4338  {
4339  // create a general error message that the export of an event
4340  // assignment failed and continue
4341  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 73, std::string("event assignment for object called \"" + pObject->getObjectName() + "\" in evend named \"" + event.getObjectName() + "\"").c_str(), "event assignment");
4342  continue;
4343  }
4344 
4345  // iii) for each assignment, check if the object that variable points to has been exported
4346  SBase* pSBase = pos->second;
4347  assert(pSBase != NULL);
4348 
4349  if (this->mHandledSBMLObjects.find(pSBase) == this->mHandledSBMLObjects.end())
4350  {
4351  // create an error message that the event assignment is for an object
4352  // that is no longer part of the model
4353  // ignore it and continue
4354  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 57, std::string("event assignment in event called \"" + event.getObjectName() + "\"").c_str(), pSBase->getId().c_str());
4355  continue;
4356  }
4357 
4358  // set the constant flag to false
4359  Compartment* pSBMLCompartment = NULL;
4360  Species* pSBMLSpecies = NULL;
4361  Parameter* pSBMLParameter = NULL;
4362 
4363  switch (pSBase->getTypeCode())
4364  {
4365  case SBML_COMPARTMENT:
4366  pSBMLCompartment = dynamic_cast<Compartment*>(pSBase);
4367  assert(pSBMLCompartment != NULL);
4368 
4369  if (pSBMLCompartment->getConstant() == true)
4370  {
4371  pSBMLCompartment->setConstant(false);
4372  }
4373 
4374  break;
4375 
4376  case SBML_SPECIES:
4377  pSBMLSpecies = dynamic_cast<Species*>(pSBase);
4378  assert(pSBMLSpecies != NULL);
4379 
4380  if (pSBMLSpecies->getConstant() == true)
4381  {
4382  pSBMLSpecies->setConstant(false);
4383  }
4384 
4385  break;
4386 
4387  case SBML_PARAMETER:
4388  pSBMLParameter = dynamic_cast<Parameter*>(pSBase);
4389  assert(pSBMLParameter != NULL);
4390 
4391  if (pSBMLParameter->getConstant() == true)
4392  {
4393  pSBMLParameter->setConstant(false);
4394  }
4395 
4396  break;
4397 
4398  default:
4399  fatalError();
4400  break;
4401  }
4402 
4403  std::string sbmlId = pSBase->getId();
4404  assert(sbmlId.find_first_not_of("\t\r\n ") != std::string::npos);
4405  std::map<std::string, EventAssignment*>::iterator pos2 = assignmentMap.find(sbmlId);
4406  EventAssignment* pAssignment = NULL;
4407 
4408  if (pos2 != assignmentMap.end())
4409  {
4410  pAssignment = dynamic_cast<EventAssignment*>(pos2->second);
4411  assert(pAssignment != NULL);
4412  pSBMLEvent->getListOfEventAssignments()->appendAndOwn(pAssignment);
4413  // remove the entry from the event assigment map
4414  assignmentMap.erase(pos2);
4415  }
4416 
4417  if (pAssignment == NULL)
4418  {
4419  pAssignment = pSBMLEvent->createEventAssignment();
4420  }
4421 
4422  pAssignment->setVariable(sbmlId);
4423  // iv) for each assignment, check the assignment expression if it is valid
4424  const CExpression * pExpression = (*itAssignment)->getExpressionPtr();
4425 
4426  if (pExpression != NULL)
4427  {
4428  std::vector<SBMLIncompatibility> result;
4430  *pExpression
4431  , dataModel
4432  , this->mSBMLLevel
4433  , this->mSBMLVersion
4434  , result
4435  , std::string("event assignment for variable \"" + sbmlId + "\" in event with id\"+" + event.getSBMLId() + "\"")
4436  , false
4437  , &mInitialValueMap);
4438 
4439  // collect directly used functions
4440  if (result.empty())
4441  {
4442  std::set<std::string> directlyUsedFunctionNames;
4443  CSBMLExporter::findDirectlyUsedFunctions(pExpression->getRoot(), directlyUsedFunctionNames);
4444  std::set<CFunction*> usedFunctions = CSBMLExporter::createFunctionSetFromFunctionNames(directlyUsedFunctionNames, CCopasiRootContainer::getFunctionList());
4445 
4446 # if defined _MSC_VER && _MSC_VER < 1201 // 1200 Identifies Visual C++ 6.0
4447  {
4448  it = usedFunctions.begin();
4449  usedFunctions.end();
4450 
4451  for (; it != end; ++it)
4452  this->mUsedFunctions.insert(*it);
4453  }
4454 #else
4455  this->mUsedFunctions.insert(usedFunctions.begin(), usedFunctions.end());
4456 #endif
4457  }
4458  else
4459  {
4460  this->mIncompatibilities.insert(this->mIncompatibilities.end(), result.begin(), result.end());
4461 
4462  if (!this->mIncompleteExport)
4463  {
4464  this->outputIncompatibilities();
4465  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 60, std::string("event assignment for variable with id \"" + sbmlId + "\"").c_str(), "event", event.getObjectName().c_str());
4466  }
4467  }
4468 
4469  //const CEvaluationNode* pOrigNode = pExpression->getRoot();
4470  const std::string& changedExpression = convertExpression(pExpression->getInfix(), mInitialValueMap);
4471  CEvaluationTree tree;
4472  tree.setInfix(changedExpression);
4473  const CEvaluationNode* pOrigNode = tree.getRoot();
4474 
4476  {
4477  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 70, std::string("event assignment for variable with id \"" + sbmlId + "\"").c_str(), "event", event.getObjectName().c_str());
4478  }
4479 
4480  // the next few lines replace references to species depending on whether
4481  // it is a reference to an amount or a reference to a concentration.
4482  // Other factors that influence this replacement are if the model
4483  // contains variable volumes or if the quantity units are set to CModel::number
4484  pOrigNode = this->replaceSpeciesReferences(pOrigNode, dataModel);
4485  assert(pOrigNode != NULL);
4486  // check if the assignment belongs to an amount species
4487  // if so, multiply the expression be the volume of the compartment that
4488  // the species belongs to.
4489  const CMetab* pMetab = dynamic_cast<const CMetab*>(pObject);
4490 
4491  if (pMetab != NULL)
4492  {
4493  std::map<const CCopasiObject*, SBase*>::const_iterator pos = this->mCOPASI2SBMLMap.find(pObject);
4494  assert(pos != this->mCOPASI2SBMLMap.end());
4495 
4496  if (dynamic_cast<const Species*>(pos->second)->getHasOnlySubstanceUnits() == true)
4497  {
4498  const CCompartment* pCompartment = pMetab->getCompartment();
4499 
4500  if (pCompartment->getDimensionality() != 0)
4501  {
4502  CEvaluationNode* pNode = CSBMLExporter::multiplyByObject(pOrigNode, pCompartment->getValueReference());
4503  assert(pNode != NULL);
4504 
4505  if (pNode != NULL)
4506  {
4507  delete pOrigNode;
4508  pOrigNode = pNode;
4509  }
4510  }
4511  }
4512  }
4513 
4514  ASTNode* pNode = this->convertToASTNode(pOrigNode, dataModel);
4515  delete pOrigNode;
4516  // convert local parameters that are referenced in an expression to global
4517  // parameters
4518  this->replace_local_parameters(pNode, dataModel);
4519 
4520  if (pNode != NULL)
4521  {
4522  pAssignment->setMath(pNode);
4523  delete pNode;
4524  }
4525  else
4526  {
4527  if (this->mIncompleteExport != true)
4528  {
4529  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 60, std::string("event assignment for variable with id \"" + sbmlId + "\"").c_str(), "event", event.getObjectName().c_str());
4530  }
4531  }
4532  }
4533  else
4534  {
4535  // create an error message that the event assignment could not be
4536  // created and continue
4537  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 73, std::string("event assignment for variable with id \"" + sbmlId + "\"").c_str(), "event assignment");
4538  continue;
4539  }
4540  }
4541 
4542  // v) delete the old unused assignments
4543  std::map<std::string, EventAssignment*>::iterator it = assignmentMap.begin(), endit = assignmentMap.end();
4544 
4545  while (it != endit)
4546  {
4547  delete it->second;
4548  ++it;
4549  }
4550 }
4551 
4552 void CSBMLExporter::checkForEvents(const CCopasiDataModel& dataModel, std::vector<SBMLIncompatibility>& result)
4553 {
4554  if (dataModel.getModel() != NULL && dataModel.getModel()->getEvents().size() > 0)
4555  {
4556  result.push_back(SBMLIncompatibility(7));
4557  }
4558 }
4559 
4561 {
4562  // make sure the idMap is already up to date
4563  // go through the existing map and create a new one with all SBML
4564  // objects updated with objects from the copied
4565  // model
4566  this->mCOPASI2SBMLMap.clear();
4567  std::map<CCopasiObject*, SBase*>& modelMap = const_cast<CCopasiDataModel&>(dataModel).getCopasi2SBMLMap();
4568  std::map<CCopasiObject*, SBase*>::const_iterator it = modelMap.begin();
4569  std::map<CCopasiObject*, SBase*>::const_iterator endit = modelMap.end();
4570 
4571  while (it != endit)
4572  {
4573  std::map<std::string, const SBase*>::iterator pos = it->second == NULL ? mIdMap.end() : this->mIdMap.find(it->second->getId());
4574 
4575  if (pos != this->mIdMap.end())
4576  {
4577  this->mCOPASI2SBMLMap.insert(std::pair<const CCopasiObject * const, SBase*>(it->first, const_cast<SBase*>(pos->second)));
4578  }
4579 
4580  ++it;
4581  }
4582 }
4583 
4584 KineticLaw* CSBMLExporter::createKineticLaw(CReaction& reaction, CCopasiDataModel& dataModel, unsigned int level, unsigned int version)
4585 {
4586  KineticLaw* pKLaw = NULL;
4587 
4588  if (!pKLaw)
4589  {
4590  /* create a new KineticLaw */
4591  pKLaw = new KineticLaw(level, version);
4592  }
4593 
4594  // create the local parameters
4595  size_t i, iMax = reaction.getFunctionParameters().size();
4596 
4597  for (i = 0; i < iMax; ++i)
4598  {
4599  const CFunctionParameter* pPara = reaction.getFunctionParameters()[i];
4600  // if the reaction calls a general function, all call parameters have a usage of VARIABLE
4601  // So local parameters will also have a usage of VARIABLE instead of PARAMETER
4602 
4603  if (pPara->getUsage() == CFunctionParameter::PARAMETER ||
4604  (reaction.getFunction() != NULL &&
4605  reaction.getFunction()->isReversible() == TriUnspecified &&
4607  {
4608  // only create a parameter if it is a local parameter,
4609  // and if it is not in the replacement map
4610  // otherwise the parameter already has been created
4611  if (reaction.isLocalParameter(i))
4612  {
4613  std::vector<std::string> v = reaction.getParameterMapping(pPara->getObjectName());
4614  assert(v.size() == 1);
4615  CCopasiObject* pTmpObject = CCopasiRootContainer::getKeyFactory()->get(v[0]);
4616  assert(pTmpObject != NULL);
4617  CCopasiParameter* pLocalParameter = dynamic_cast<CCopasiParameter*>(pTmpObject);
4618  assert(pLocalParameter != NULL);
4619 
4620  if (this->mParameterReplacementMap.find(pLocalParameter->getCN()) == this->mParameterReplacementMap.end())
4621  {
4622  // starting with SBML Level 3, the parameters of a kinetic law are expressed in
4623  // libsbml 4.1 does not handle this in a nice way,
4624  // so way have to distinguish between the levels we export
4625  // Actually here we could probably have gotten away with the old code, but
4626  // better safe than sorry
4627  Parameter* pSBMLPara = NULL;
4628 #if LIBSBML_VERSION >= 40100
4629 
4630  if (this->mSBMLLevel > 2)
4631  {
4632  pSBMLPara = pKLaw->createLocalParameter();
4633  }
4634  else
4635  {
4636 #endif // LIBSBML_VERSION
4637  pSBMLPara = pKLaw->createParameter();
4638 #if LIBSBML_VERSION >= 40100
4639  }
4640 
4641  pSBMLPara->setId(pPara->getObjectName().c_str());
4642 #endif // LIBSBML_VERSION
4643 
4644  // don't call setName on level 1 objects because this will also
4645  // change the id
4646  if (this->mpSBMLDocument->getLevel() > 1)
4647  {
4648  pSBMLPara->setName(pPara->getObjectName().c_str());
4649  }
4650 
4651  double value = reaction.getParameterValue(pPara->getObjectName());
4652 
4653  // if the value is NaN, leave the parameter value unset.
4654  if (!isnan(value))
4655  {
4656  pSBMLPara->setValue(value);
4657  }
4658  }
4659  }
4660  }
4661  }
4662 
4663  // the next few lines replace references to species depending on whether
4664  // it is a reference to an amount or a reference to a concentration.
4665  // Other factors that influence this replacement are if the model
4666  // contains variable volumes or if the quantity units are set to CModel::number
4667  CEvaluationNode* pExpression = CSBMLExporter::createKineticExpression(const_cast<CFunction*>(reaction.getFunction()), reaction.getParameterMappings());
4668 
4669  if (pExpression == NULL)
4670  {
4671  delete pKLaw;
4672  pKLaw = NULL;
4673  }
4674  else
4675  {
4676  /*
4677  ** If the reaction takes place in a single compartment, the rate law has
4678  ** to be converted from concentration/time to substance/time by
4679  ** multiplying the rate law with the volume of the compartment.
4680  */
4681  CEvaluationNode* pOrigNode = this->replaceSpeciesReferences(pExpression, dataModel);
4682  delete pExpression;
4683  assert(pOrigNode != NULL);
4684  ASTNode* pNode = this->convertToASTNode(pOrigNode, dataModel);
4685  // since toAST in CEvaluationNodeObject sets the CN of local parameters
4686  // as the name of the node, we have to go and replace all node names
4687  // that have a common name that is the CN of a local parameter
4688  this->restore_local_parameters(pNode, dataModel);
4689  delete pOrigNode;
4690  assert(pNode != NULL);
4691 
4692  if (reaction.getCompartmentNumber() == 1)
4693  {
4694  const CCompartment& compartment = (reaction.getChemEq().getSubstrates().size() != 0) ? (*reaction.getChemEq().getSubstrates()[0]->getMetabolite()->getCompartment()) : (*reaction.getChemEq().getProducts()[0]->getMetabolite()->getCompartment());
4695 
4696  if (compartment.getDimensionality() != 0)
4697  {
4698  // check if the importer has added a division by the volume
4699  // if so remove it instead of multiplying again
4700  ASTNode* pTNode = CSBMLExporter::isDividedByVolume(pNode, compartment.getSBMLId());
4701 
4702  if (pTNode)
4703  {
4704  if (pTNode->getNumChildren() == 0)
4705  {
4706  fatalError();
4707  }
4708 
4709  if (pTNode->getNumChildren() == 1)
4710  {
4711  ASTNode* pTmp = static_cast<ConverterASTNode*>(pTNode)->removeChild(0);
4712  delete pTNode;
4713  pTNode = pTmp;
4714  }
4715 
4716  delete pNode;
4717  pNode = pTNode;
4718  }
4719  else
4720  {
4721  pTNode = new ASTNode(AST_TIMES);
4722  ASTNode* pVNode = new ASTNode(AST_NAME);
4723  pVNode->setName(compartment.getSBMLId().c_str());
4724  pTNode->addChild(pVNode);
4725  pTNode->addChild(pNode);
4726  pNode = pTNode;
4727  }
4728  }
4729  }
4730 
4731  pKLaw->setMath(pNode);
4732  delete pNode;
4733  }
4734 
4735  return pKLaw;
4736 }
4737 
4738 ASTNode* CSBMLExporter::isDividedByVolume(const ASTNode* pRootNode, const std::string& compartmentId)
4739 {
4740  ASTNode* pResult = NULL;
4741 
4742  if (pRootNode->getType() == AST_DIVIDE || pRootNode->getType() == AST_TIMES)
4743  {
4744  ASTNode* pTmpResultNode = new ConverterASTNode(pRootNode->getType());
4745  unsigned int i, iMax = pRootNode->getNumChildren();
4746  bool found = false;
4747 
4748  for (i = 0; i < iMax; ++i)
4749  {
4750  const ASTNode* pChild = pRootNode->getChild(i);
4751 
4752  if (pRootNode->getType() == AST_DIVIDE && pChild->getType() == AST_NAME && pChild->getName() == compartmentId)
4753  {
4754  found = true;
4755  }
4756  else if ((!found) && (pChild->getType() == AST_DIVIDE || pChild->getType() == AST_TIMES))
4757  {
4758  ASTNode* pSubResult = CSBMLExporter::isDividedByVolume(pChild, compartmentId);
4759 
4760  if (pSubResult != NULL)
4761  {
4762  found = true;
4763 
4764  if (pSubResult->getNumChildren() > 1)
4765  {
4766  pTmpResultNode->addChild(pSubResult);
4767  }
4768  else if (pSubResult->getNumChildren() == 1)
4769  {
4770  pTmpResultNode->addChild(static_cast<ASTNode*>(static_cast<ConverterASTNode*>(pSubResult)->removeChild(0)));
4771  delete pSubResult;
4772  }
4773  else
4774  {
4775  delete pSubResult;
4776  }
4777  }
4778  else
4779  {
4780  pTmpResultNode->addChild(new ConverterASTNode(*pChild));
4781  }
4782  }
4783  else
4784  {
4785  pTmpResultNode->addChild(new ConverterASTNode(*pChild));
4786  }
4787  }
4788 
4789  if (found)
4790  {
4791  pResult = pTmpResultNode;
4792  }
4793  else
4794  {
4795  delete pTmpResultNode;
4796  }
4797  }
4798 
4799  return pResult;
4800 }
4801 
4802 CEvaluationNode* CSBMLExporter::createKineticExpression(CFunction* pFun, const std::vector<std::vector<std::string> >& arguments)
4803 {
4804  if (!pFun || pFun->getVariables().size() != arguments.size()) fatalError();
4805 
4807 
4808  if (pFun->getType() == CEvaluationTree::MassAction)
4809  {
4810  pResult = CSBMLExporter::createMassActionExpression(arguments, pFun->isReversible() == TriTrue);
4811  }
4812  else
4813  {
4815  this->mUsedFunctions.insert(pFun);
4816  size_t i, iMax = arguments.size();
4817  std::string cn;
4818 
4819  for (i = 0; i < iMax; ++i)
4820  {
4821  if (arguments[i].size() != 1) fatalError(); // we can't have arrays here.
4822 
4823  const CCopasiObject* pObject = CCopasiRootContainer::getKeyFactory()->get(arguments[i][0]);
4824 
4825  if (!pObject) fatalError();
4826 
4827  if (dynamic_cast<const CModel*>(pObject) != NULL)
4828  {
4829  cn = "<" + pObject->getCN() + ",Reference=Time>";
4830  }
4831  else if (dynamic_cast<const CCompartment*>(pObject) != NULL)
4832  {
4833  cn = "<" + pObject->getCN() + ",Reference=Volume>";
4834  }
4835  else if (dynamic_cast<const CMetab*>(pObject) != NULL)
4836  {
4837  cn = "<" + pObject->getCN() + ",Reference=Concentration>";
4838  }
4839  else if (dynamic_cast<const CModelValue*>(pObject) != NULL)
4840  {
4841  cn = "<" + pObject->getCN() + ",Reference=Value>";
4842  }
4843  else if (dynamic_cast<const CReaction*>(pObject) != NULL)
4844  {
4845  cn = "<" + pObject->getCN() + ",Reference=Flux>";
4846  }
4847  else if (dynamic_cast<const CCopasiParameter*>(pObject) != NULL)
4848  {
4849  // local parameter of a reaction
4850  // must set the node content to the SBML id if the parameter
4851  cn = "<" + pObject->getCN() + ">";
4852  }
4853  else
4854  {
4855  cn = "<" + pObject->getCN() + ">";
4856  }
4857 
4858  pFunctionCall->addChild(new CEvaluationNodeObject(CEvaluationNodeObject::CN, cn));
4859  }
4860 
4861  pResult = pFunctionCall;
4862  }
4863 
4864  return pResult;
4865 }
4866 
4867 /**
4868  * This method creates an ASTNode tree where all the species specified in
4869  * the given vector are multiplied. This is used to create the mass action
4870  * kinetic law.
4871  */
4873 {
4874  ASTNode* pNode = NULL;
4875  double multiplicity = vect[pos]->getMultiplicity();
4876 
4877  if (pos == vect.size() - 1)
4878  {
4879  pNode = new ASTNode(AST_NAME);
4880  const CMetab* pMetab = vect[pos]->getMetabolite();
4881  assert(pMetab);
4882 
4883  pNode->setName(pMetab->getSBMLId().c_str());
4884 
4885  /* if the stoichiometry is not 1.0, we have to add it to the exponent */
4886  if (multiplicity != 1.0)
4887  {
4888  ASTNode* pTmpNode1 = new ASTNode(AST_POWER);
4889  ASTNode* pTmpNode2 = new ASTNode(AST_REAL);
4890  pTmpNode2->setValue(multiplicity);
4891  pTmpNode1->addChild(pNode);
4892  pTmpNode1->addChild(pTmpNode2);
4893  pNode = pTmpNode1;
4894  }
4895  }
4896  else
4897  {
4898  pNode = new ASTNode(AST_TIMES);
4899  ASTNode* pChild = new ASTNode(AST_NAME);
4900  const CMetab* pMetab = vect[pos]->getMetabolite();
4901  assert(pMetab);
4902 
4903  pChild->setName(pMetab->getSBMLId().c_str());
4904 
4905  /* if the stoichiometry is not 1.0, we have to add it to the exponent */
4906  if (multiplicity != 1.0)
4907  {
4908  ASTNode* pTmpNode1 = new ASTNode(AST_POWER);
4909  ASTNode* pTmpNode2 = new ASTNode(AST_REAL);
4910  pTmpNode2->setValue(multiplicity);
4911  pTmpNode1->addChild(pChild);
4912  pTmpNode1->addChild(pTmpNode2);
4913  pChild = pTmpNode1;
4914  }
4915 
4916  pNode->addChild(pChild);
4917  pNode->addChild(CSBMLExporter::createTimesTree(vect, pos + 1));
4918  }
4919 
4920  return pNode;
4921 }
4922 
4923 /**
4924  * Go through a CEvaluationNode base tree and add the names
4925  * of all functions directly called in this tree to the set.
4926  */
4927 void CSBMLExporter::findDirectlyUsedFunctions(const CEvaluationNode* pRootNode, std::set<std::string>& result)
4928 {
4929  if (pRootNode == NULL) return;
4930 
4932  {
4933  result.insert(pRootNode->getData());
4934  }
4935 
4936  const CEvaluationNode* pChild = dynamic_cast<const CEvaluationNode*>(pRootNode->getChild());
4937 
4938  while (pChild != NULL)
4939  {
4941  pChild = dynamic_cast<const CEvaluationNode*>(pChild->getSibling());
4942  }
4943 }
4944 
4945 const std::vector<CFunction*> CSBMLExporter::findUsedFunctions(std::set<CFunction*>& functions, CFunctionDB* pFunctionDB)
4946 {
4947  std::map<CFunction*, std::set<CFunction*> > functionDependencyMap;
4948  std::set<CFunction*>::iterator setit = functions.begin(), setendit = functions.end();
4949 
4950  // create the dependency map
4951  while (setit != setendit)
4952  {
4953  std::set<std::string> usedFunctionNames;
4954  CSBMLExporter::findDirectlyUsedFunctions(const_cast<CFunction*>(*setit)->getRoot(), usedFunctionNames);
4955  std::set<CFunction*> usedFunctions = CSBMLExporter::createFunctionSetFromFunctionNames(usedFunctionNames, pFunctionDB);
4956  functionDependencyMap[*setit] = usedFunctions;
4957  ++setit;
4958  }
4959 
4960  bool found = true;
4961 
4962  while (found)
4963  {
4964  found = false;
4965  std::map<CFunction*, std::set<CFunction*> >::iterator mapIt = functionDependencyMap.begin(), mapEndit = functionDependencyMap.end();
4966 
4967  while (mapIt != mapEndit)
4968  {
4969  setit = mapIt->second.begin(), setendit = mapIt->second.end();
4970 
4971  while (setit != setendit)
4972  {
4973  // check if we already have the dependencies for the function
4974  if (functionDependencyMap.find(*setit) == functionDependencyMap.end())
4975  {
4976  found = true;
4977  std::set<std::string> usedFunctionNames;
4978  CSBMLExporter::findDirectlyUsedFunctions(const_cast<CFunction*>(*setit)->getRoot(), usedFunctionNames);
4979  std::set<CFunction*> usedFunctions = CSBMLExporter::createFunctionSetFromFunctionNames(usedFunctionNames, pFunctionDB);
4980  functionDependencyMap[*setit] = usedFunctions;
4981  }
4982 
4983  ++setit;
4984  }
4985 
4986  // if we changed the map, we have to stop the inner loop since the
4987  // map iterator is now invalid
4988  if (found)
4989  {
4990  break;
4991  }
4992 
4993  ++mapIt;
4994  }
4995  }
4996 
4997  bool removed = true;
4998  std::vector<CFunction*> functionVector;
4999 
5000  while (!functionDependencyMap.empty() && removed == true)
5001  {
5002  std::map<CFunction*, std::set<CFunction*> >::iterator mapIt = functionDependencyMap.begin(), mapEndit = functionDependencyMap.end();
5003  removed = false;
5004 
5005  // go through the map and find functions that have no dependencies
5006  while (mapIt != mapEndit)
5007  {
5008  if (mapIt->second.empty())
5009  {
5010  // set remove to true so that we know that at least one
5011  // function definition without dependencies has been found
5012  // otherwise we might end up in an endless loop
5013  removed = true;
5014  // add the function to the front
5015  functionVector.push_back(mapIt->first);
5016  // remove this function from the dependency list of all other
5017  // functions
5018  std::map<CFunction*, std::set<CFunction*> >::iterator mapIt2 = functionDependencyMap.begin(), mapEndit2 = functionDependencyMap.end();
5019 
5020  while (mapIt2 != mapEndit2)
5021  {
5022  std::set<CFunction*>::iterator pos = mapIt2->second.find(mapIt->first);
5023 
5024  if (pos != mapIt2->second.end())
5025  {
5026  mapIt2->second.erase(pos);
5027  }
5028 
5029  ++mapIt2;
5030  }
5031 
5032  // delete the entry of this function from the depenency map
5033  std::map<CFunction*, std::set<CFunction*> >::iterator pos = mapIt;
5034  ++mapIt;
5035  functionDependencyMap.erase(pos);
5036  }
5037  else
5038  {
5039  ++mapIt;
5040  }
5041  }
5042  }
5043 
5044  // check if the dependency map is empty. If not, there must be a dependency
5045  // loop
5046  if (!functionDependencyMap.empty())
5047  {
5048  CCopasiMessage(CCopasiMessage::EXCEPTION, "Dependency loop found in function definitions.");
5049  }
5050 
5051  return functionVector;
5052 }
5053 
5054 const std::set<CFunction*> CSBMLExporter::createFunctionSetFromFunctionNames(const std::set<std::string>& names, CFunctionDB* pFunDB)
5055 {
5056  std::set<CFunction*> result;
5057 
5058  if (pFunDB != NULL)
5059  {
5060  std::set<std::string>::const_iterator it = names.begin();
5061  std::set<std::string>::const_iterator endit = names.end();
5062  CFunction* pFun = NULL;
5063 
5064  while (it != endit)
5065  {
5066  pFun = dynamic_cast<CFunction*>(pFunDB->findFunction(*it));
5067 
5068  if (pFun == NULL)
5069  {
5070  CCopasiMessage(CCopasiMessage::ERROR, MCSBML + 15, (*it).c_str());
5071  }
5072  else
5073  {
5074  result.insert(pFun);
5075  }
5076 
5077  ++it;
5078  }
5079  }
5080 
5081  return result;
5082 }
5083 
5084 std::vector<CModelEntity*> CSBMLExporter::orderRules(const CCopasiDataModel& dataModel)
5085 {
5086  // make sure the vector of rules is filled
5087  std::map<CModelEntity*, std::set<const CModelEntity*> > dependencyMap;
5088  std::vector<CModelEntity*>::iterator it = this->mAssignmentVector.begin(), endit = this->mAssignmentVector.end();
5089 
5090  while (it != endit)
5091  {
5092  const CExpression* pExpr = (*it)->getExpressionPtr();
5093  assert(pExpr != NULL);
5094  std::set<const CModelEntity*> dependencies;
5095  CSBMLExporter::findModelEntityDependencies(pExpr->getRoot(), dataModel, dependencies);
5096  /*
5097  std::vector<CModelEntity*>::iterator it2 = this->mAssignmentVector.begin(), endit2 = this->mAssignmentVector.end();
5098  while (it2 != endit2)
5099  {
5100  std::set<const CModelEntity*>::iterator pos = dependencies.find(*it2);
5101  if (pos != dependencies.end())
5102  {
5103  dependencies.erase(pos);
5104  }
5105  ++it2;
5106  }
5107  */
5108  dependencyMap[*it] = dependencies;
5109  ++it;
5110  }
5111 
5112  // now sort the entities according to their dependencies
5113  bool error = false;
5114  std::vector<CModelEntity*> orderedRules;
5115  // remove all dependencies that don't have assignment rules themselves
5116  std::map<CModelEntity*, std::set<const CModelEntity*> > reducedDependencyMap;
5117  std::map<CModelEntity*, std::set<const CModelEntity*> >::iterator mapIt = dependencyMap.begin(), mapEndit = dependencyMap.end();
5118 
5119  while (mapIt != mapEndit)
5120  {
5121  std::set<const CModelEntity*> dependencies;
5122  std::set<const CModelEntity*>::const_iterator setIt = mapIt->second.begin(), setEndit = mapIt->second.end();
5123 
5124  while (setIt != setEndit)
5125  {
5126  if (dependencyMap.find(const_cast<CModelEntity*>(*setIt)) != dependencyMap.end())
5127  {
5128  dependencies.insert(*setIt);
5129  }
5130 
5131  ++setIt;
5132  }
5133 
5134  reducedDependencyMap[mapIt->first] = dependencies;
5135  ++mapIt;
5136  }
5137 
5138  while (error == false && !reducedDependencyMap.empty())
5139  {
5140  mapIt = reducedDependencyMap.begin(), mapEndit = reducedDependencyMap.end();
5141  // if we cant't remove anything in one round, we have a circular
5142  // dependency
5143  error = true;
5144 
5145  while (mapIt != mapEndit)
5146  {
5147  // if the dependency list of the item is not empty
5148  // we remove all entities that have been removed already
5149  if (!mapIt->second.empty())
5150  {
5151  std::vector<CModelEntity*>::iterator vIt = orderedRules.begin(), vEndit = orderedRules.end();
5152 
5153  while (vIt != vEndit)
5154  {
5155  std::set<const CModelEntity*>::iterator pos = mapIt->second.find(*vIt);
5156 
5157  if (pos != mapIt->second.end())
5158  {
5159  mapIt->second.erase(pos);
5160  }
5161 
5162  ++vIt;
5163  }
5164  }
5165 
5166  if (mapIt->second.empty())
5167  {
5168  // since the dependencies are empty, we can add the rule to be
5169  // exported
5170  orderedRules.push_back(mapIt->first);
5171  // now we know that we could remove something in this round
5172  error = false;
5173  break;
5174  }
5175 
5176  ++mapIt;
5177  }
5178 
5179  if (error == false)
5180  {
5181  reducedDependencyMap.erase(mapIt);
5182  }
5183  }
5184 
5185  return orderedRules;
5186 }
5187 
5188 void CSBMLExporter::findModelEntityDependencies(const CEvaluationNode* pNode, const CCopasiDataModel& dataModel, std::set<const CModelEntity*>& dependencies)
5189 {
5190  if (pNode == NULL) return;
5191 
5193  {
5194  const CEvaluationNodeObject* pObjectNode = dynamic_cast<const CEvaluationNodeObject*>(pNode);
5195  assert(pObjectNode != NULL);
5196 
5197  if (pObjectNode != NULL)
5198  {
5199  const CCopasiObject* pObject = dataModel.getDataObject(pObjectNode->getObjectCN());
5200 
5201  if (!pObject)
5202  {
5203  fatalError();
5204  }
5205 
5206  if (pObject->isReference())
5207  {
5208  pObject = pObject->getObjectParent();
5209  assert(pObject != NULL);
5210  }
5211 
5212  const CModelEntity* pME = dynamic_cast<const CModelEntity*>(pObject);
5213 
5214  if (pME != NULL)
5215  {
5216  dependencies.insert(pME);
5217  }
5218  }
5219  }
5220 
5221  const CEvaluationNode* pChild = dynamic_cast<const CEvaluationNode*>(pNode->getChild());
5222 
5223  while (pChild != NULL)
5224  {
5225  CSBMLExporter::findModelEntityDependencies(pChild, dataModel, dependencies);
5226  pChild = dynamic_cast<const CEvaluationNode*>(pChild->getSibling());
5227  }
5228 }
5229 
5231 {
5232  // expand all function calls in rules, events and
5233  // kinetic laws and delete the functions
5234  // initial assignments do not need to be considered since they can not be
5235  // exported to Level 1 anyway
5236  //
5237  // checking of the resulting expressions for piecewise functions
5238  // is done in the convertASTTreeToLevel1 method
5239  if (this->mpSBMLDocument == NULL) return;
5240 
5241  Model* pModel = this->mpSBMLDocument->getModel();
5242  Rule* pRule = NULL;
5243  Reaction* pReaction = NULL;
5244  KineticLaw* pLaw = NULL;
5245  unsigned int i, iMax = pModel->getNumRules();
5246 
5247  for (i = 0; i < iMax; ++i)
5248  {
5249  pRule = pModel->getRule(i);
5250  assert(pRule != NULL);
5251  const ASTNode* pMath = pRule->getMath();
5252  assert(pMath != NULL);
5253  std::string message = "rule for object with id \"";
5254  message += pRule->getVariable();
5255  message += "\"";
5256  ASTNode* pNewMath = CSBMLExporter::convertASTTreeToLevel1(pMath, this->mExportedFunctions, message);
5257  assert(pNewMath != NULL);
5258 
5259  if (pNewMath != NULL)
5260  {
5261  pRule->setMath(pNewMath);
5262  delete pNewMath;
5263  }
5264  }
5265 
5266  iMax = pModel->getNumReactions();
5267 
5268  for (i = 0; i < iMax; ++i)
5269  {
5270  pReaction = pModel->getReaction(i);
5271  assert(pReaction != NULL);
5272  pLaw = pReaction->getKineticLaw();
5273 
5274  // maybe we don't have a kinetic law for all reactions
5275  if (pLaw != NULL)
5276  {
5277  const ASTNode* pMath = pLaw->getMath();
5278  assert(pMath != NULL);
5279  std::string message = "kinetic law in reaction with id \"";
5280  message += pReaction->getId();
5281  message += "\"";
5282  ASTNode* pNewMath = CSBMLExporter::convertASTTreeToLevel1(pMath, this->mExportedFunctions, message);
5283  assert(pNewMath != NULL);
5284 
5285  if (pNewMath != NULL)
5286  {
5287  pLaw->setMath(pNewMath);
5288  delete pNewMath;
5289  }
5290  else
5291  {
5292  fatalError();
5293  }
5294  }
5295  }
5296 }
5297 
5298 ASTNode* CSBMLExporter::replaceL1IncompatibleNodes(const ASTNode* pNode)
5299 {
5300  // replace all inf, pi, nan and exponentiale nodes
5301  // replace all unsupported functions calls: SEC,CSC, COT, SINH, COSH,
5302  // TANH, SECH, CSCH, COTH, ARCSINH, ARCCOSH, ARCTANH, ARCSECH, ARCSCSH,
5303  // ARCCOTH
5304  if (pNode == NULL) return NULL;
5305 
5306  ASTNode* pResult = NULL;
5307  ASTNode* pChild = NULL;
5308  unsigned int i, iMax;
5309 
5310  switch (pNode->getType())
5311  {
5312  case AST_FUNCTION_PIECEWISE:
5313  break;
5314 
5315  case AST_CONSTANT_E:
5316  // replace by exp(1)
5317  pResult = new ASTNode(AST_FUNCTION_EXP);
5318  pChild = new ASTNode(AST_REAL);
5319  pChild->setValue(1.0);
5320  pResult->addChild(pChild);
5321  break;
5322 
5323  case AST_CONSTANT_PI:
5324  // replace by 2*ASIN(1.0)
5325  pResult = new ASTNode(AST_TIMES);
5326  pChild = new ASTNode(AST_REAL);
5327  pChild->setValue(2.0);
5328  pResult->addChild(pChild);
5329  pResult->addChild(new ASTNode(AST_FUNCTION_ARCSIN));
5330  pChild = new ASTNode(AST_REAL);
5331  pChild->setValue(1.0);
5332  pResult->getChild(1)->addChild(pChild);
5333  break;
5334 
5335  case AST_FUNCTION_SEC:
5336  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5337  assert(pChild);
5338  pResult = replace_SEC(pChild);
5339  // replace by 1/cos(X)
5340  delete pChild;
5341  break;
5342 
5343  case AST_FUNCTION_CSC:
5344  // replace by 1/sin(X)
5345  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5346  assert(pChild);
5347  pResult = replace_CSC(pChild);
5348  delete pChild;
5349  break;
5350 
5351  case AST_FUNCTION_COT:
5352  // replace by 1/tan(X)
5353  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5354  assert(pChild);
5355  pResult = replace_COT(pChild);
5356  delete pChild;
5357  break;
5358 
5359  case AST_FUNCTION_SINH:
5360  // replace by (e^X-e^(-X))/2
5361  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5362  assert(pChild);
5363  pResult = replace_SINH(pChild);
5364  delete pChild;
5365  break;
5366 
5367  case AST_FUNCTION_COSH:
5368  // replace by (e^X+e^(-X))/2
5369  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5370  assert(pChild);
5371  pResult = replace_COSH(pChild);
5372  delete pChild;
5373  break;
5374 
5375  case AST_FUNCTION_TANH:
5376  // replace by (e^X-e^(-X))/(e^X+e^(-X))
5377  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5378  assert(pChild);
5379  pResult = replace_TANH(pChild);
5380  delete pChild;
5381  break;
5382 
5383  case AST_FUNCTION_SECH:
5384  // replace by 2/(e^X+e^(-X))
5385  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5386  assert(pChild);
5387  pResult = replace_SECH(pChild);
5388  delete pChild;
5389  break;
5390 
5391  case AST_FUNCTION_CSCH:
5392  // replace by 2/(e^X-e^(-X))
5393  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5394  assert(pChild);
5395  pResult = replace_CSCH(pChild);
5396  delete pChild;
5397  break;
5398 
5399  case AST_FUNCTION_COTH:
5400  // replace by (e^X+e^(-X))/(e^X-e^(-X))
5401  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5402  assert(pChild);
5403  pResult = replace_COTH(pChild);
5404  delete pChild;
5405  break;
5406 
5407  case AST_FUNCTION_ARCSINH:
5408  // replace by log(X + sqrt(X^2 + 1))
5409  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5410  assert(pChild);
5411  pResult = replace_ARCSINH(pChild);
5412  delete pChild;
5413  break;
5414 
5415  case AST_FUNCTION_ARCCOSH:
5416  // replace by log(X + sqrt(X-1) * sqrt(X+1))
5417  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5418  assert(pChild);
5419  pResult = replace_ARCCOSH(pChild);
5420  delete pChild;
5421  break;
5422 
5423  case AST_FUNCTION_ARCTANH:
5424  // replace by 1/2 * (log(1+X) - log(1-X))
5425  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5426  assert(pChild);
5427  pResult = replace_ARCTANH(pChild);
5428  delete pChild;
5429  break;
5430 
5431  case AST_FUNCTION_ARCSECH:
5432  // replace by log(sqrt((1/X)-1) * sqrt(1+(1/X)) + 1/X)
5433  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5434  assert(pChild);
5435  pResult = replace_ARCSECH(pChild);
5436  delete pChild;
5437  break;
5438 
5439  case AST_FUNCTION_ARCCSCH:
5440  // replace by log(sqrt(1+ (1/ (X^2)))+(1/X))
5441  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(0));
5442  assert(pChild);
5443  pResult = replace_ARCCSCH(pChild);
5444  delete pChild;
5445  break;
5446 
5447  case AST_FUNCTION_ARCCOTH:
5448  // this one is difficult since it would need a piecewise definition which
5449  // is not available in SBML Level 1
5450  fatalError();
5451  break;
5452 
5453  case AST_REAL:
5454 
5455  // for nan and inf
5456  if (pNode->getReal() != pNode->getReal()) // NaN
5457  {
5458  // replace by 0/0
5459  pResult = new ASTNode(AST_RATIONAL);
5460  pResult->setValue(0L, 0L);
5461  }
5462  else if (pNode->getReal() > std::numeric_limits< C_FLOAT64 >::max())
5463  {
5464  // replace by 1/0
5465  pResult = new ASTNode(AST_RATIONAL);
5466  pResult->setValue(1L, 0L);
5467  }
5468  else if (pNode->getReal() < - std::numeric_limits< C_FLOAT64 >::max())
5469  {
5470  // replace by -1/0
5471  pResult = new ASTNode(AST_RATIONAL);
5472  pResult->setValue(-1L, 0L);
5473  }
5474  else
5475  {
5476  pResult = pNode->deepCopy();
5477  }
5478 
5479  break;
5480 
5481  default:
5482  // copy the node
5483  pResult = ConverterASTNode::shallowCopy(pNode);
5484  iMax = pNode->getNumChildren();
5485 
5486  for (i = 0; i < iMax; ++i)
5487  {
5488  pChild = CSBMLExporter::replaceL1IncompatibleNodes(pNode->getChild(i));
5489 
5490  if (pChild != NULL)
5491  {
5492  pResult->addChild(pChild);
5493  }
5494  else
5495  {
5496  delete pResult;
5497  pResult = NULL;
5498  }
5499  }
5500 
5501  break;
5502  }
5503 
5504  return pResult;
5505 }
5506 
5507 ASTNode* CSBMLExporter::convertASTTreeToLevel1(const ASTNode* pNode, const ListOfFunctionDefinitions& functions, std::string& message)
5508 {
5509  ASTNode* pExpanded = create_expression(pNode, &functions);
5510 
5511  if (pExpanded != NULL)
5512  {
5513  ASTNode* pReplaced = CSBMLExporter::replaceL1IncompatibleNodes(pExpanded);
5514  delete pExpanded;
5515 
5516  if (pReplaced == NULL)
5517  {
5518  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 62, message.c_str());
5519  }
5520 
5521  pExpanded = pReplaced;
5522  }
5523  else
5524  {
5525  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 61, message.c_str());
5526  }
5527 
5528  return pExpanded;
5529 }
5530 
5531 void CSBMLExporter::checkForPiecewiseFunctions(const CCopasiDataModel& dataModel, std::vector<SBMLIncompatibility>& result)
5532 {
5533  std::set<std::string> usedFunctionNames;
5534  const CModel* pModel = dataModel.getModel();
5535 
5536  if (pModel == NULL) return;
5537 
5538  // go through all model entities and check assignments
5539  // we don't have to check initial assignment since those can not be exported
5540  // anyway
5541  const CModelEntity* pME = NULL;
5542  const CCopasiVectorNS<CCompartment>& compartments = pModel->getCompartments();
5543  size_t i, iMax = compartments.size();
5544 
5545  for (i = 0; i < iMax; ++i)
5546  {
5547  pME = compartments[i];
5548 
5550  {
5551  const CEvaluationTree* pTree = pME->getExpressionPtr();
5552  CSBMLExporter::findDirectlyUsedFunctions(pTree->getRoot(), usedFunctionNames);
5553  CSBMLExporter::checkForPiecewiseFunctions(*pTree->getRoot(), result, pME->getObjectName(), "rule for compartment");
5554  }
5555  }
5556 
5557  const CCopasiVector<CMetab>& metabolites = pModel->getMetabolites();
5558 
5559  iMax = metabolites.size();
5560 
5561  for (i = 0; i < iMax; ++i)
5562  {
5563  pME = metabolites[i];
5564 
5566  {
5567  const CEvaluationTree* pTree = pME->getExpressionPtr();
5568  CSBMLExporter::findDirectlyUsedFunctions(pTree->getRoot(), usedFunctionNames);
5569  CSBMLExporter::checkForPiecewiseFunctions(*pTree->getRoot(), result, pME->getObjectName(), "rule for metabolite");
5570  }
5571  }
5572 
5573  const CCopasiVectorN<CModelValue>& modelvalues = pModel->getModelValues();
5574 
5575  iMax = modelvalues.size();
5576 
5577  for (i = 0; i < iMax; ++i)
5578  {
5579  pME = modelvalues[i];
5580 
5582  {
5583  const CEvaluationTree* pTree = pME->getExpressionPtr();
5584  CSBMLExporter::findDirectlyUsedFunctions(pTree->getRoot(), usedFunctionNames);
5585  CSBMLExporter::checkForPiecewiseFunctions(*pTree->getRoot(), result, pME->getObjectName(), "rule for global parameter");
5586  }
5587  }
5588 
5589  // go through all reactions and check the function called in the kinetic
5590  // laws.
5591  CReaction* pReaction = NULL;
5592  const CCopasiVectorNS<CReaction>& reactions = pModel->getReactions();
5593  iMax = reactions.size();
5594 
5595  for (i = 0; i < iMax; ++i)
5596  {
5597  pReaction = reactions[i];
5598 
5599  if (pReaction->getFunction() != NULL)
5600  {
5601  CSBMLExporter::findDirectlyUsedFunctions(pReaction->getFunction()->getRoot(), usedFunctionNames);
5602  }
5603  }
5604 
5605  // check indirectly called functions
5606  std::set<CFunction*> directlyUsedFunctions = CSBMLExporter::createFunctionSetFromFunctionNames(usedFunctionNames, CCopasiRootContainer::getFunctionList());
5607  std::vector<CFunction*> functions = CSBMLExporter::findUsedFunctions(directlyUsedFunctions, CCopasiRootContainer::getFunctionList());
5608  std::vector<CFunction*>::const_iterator it = functions.begin(), endit = functions.end();
5609 
5610  while (it != endit)
5611  {
5612  CSBMLExporter::checkForPiecewiseFunctions(*(*it)->getRoot(), result, (*it)->getObjectName(), "function");
5613  ++it;
5614  }
5615 }
5616 
5617 void CSBMLExporter::checkForPiecewiseFunctions(const CEvaluationNode& node, std::vector<SBMLIncompatibility>& result, const std::string& objectName, const std::string& objectType)
5618 {
5619 
5620  // no need to go through the children.
5621  // one incompatibility warning is enough
5623  {
5624  result.push_back(SBMLIncompatibility(8, objectType.c_str(), objectName.c_str()));
5625  }
5626  else
5627  {
5628  // store the old size
5629  // if one error has been found, we can stop looking at the other
5630  // children
5631  size_t size = result.size();
5632  const CEvaluationNode* pChild = dynamic_cast<const CEvaluationNode*>(node.getChild());
5633 
5634  while (pChild != NULL && result.size() == size)
5635  {
5636  CSBMLExporter::checkForPiecewiseFunctions(*pChild, result, objectName, objectType);
5637  pChild = dynamic_cast<const CEvaluationNode*>(pChild->getSibling());
5638  }
5639  }
5640 }
5641 
5642 /**
5643  * Remove all compartments, species, parameters and reactions and initial
5644  * assignments that did not end up in mHandledSBMLObjects during an export.
5645  */
5647 {
5648  if (this->mpSBMLDocument != NULL || this->mpSBMLDocument->getModel() == NULL)
5649  {
5650  // we need to reverse the COPASI2SBMLMap so that we can remove objects
5651  // from the map if we delete them
5652  std::map<const SBase*, const CCopasiObject*> reverseCOPASI2SBMLMap;
5653  std::map<const CCopasiObject*, SBase*>::const_iterator it = this->mCOPASI2SBMLMap.begin(), endit = this->mCOPASI2SBMLMap.end();
5654 
5655  while (it != endit)
5656  {
5657  reverseCOPASI2SBMLMap.insert(std::pair<const SBase*, const CCopasiObject*>(it->second, it->first));
5658  ++it;
5659  }
5660 
5661  Model* pModel = this->mpSBMLDocument->getModel();
5662  std::list<SBase*> removedObjects;
5663  SBase* pSBase = NULL;
5664  ListOf* pList = pModel->getListOfCompartments();
5665  unsigned int i;
5666  std::map<const SBase*, const CCopasiObject*>::const_iterator pos;
5667  std::map<const CCopasiObject*, SBase*>::iterator pos2;
5668 
5669  for (i = pList->size(); i != 0; --i)
5670  {
5671  pSBase = pList->get(i - 1);
5672 
5673  if (this->mHandledSBMLObjects.find(pSBase) == this->mHandledSBMLObjects.end())
5674  {
5675  removedObjects.push_back(pList->remove(i - 1));
5676  }
5677  }
5678 
5679  pList = pModel->getListOfSpecies();
5680 
5681  for (i = pList->size(); i != 0; --i)
5682  {
5683  pSBase = pList->get(i - 1);
5684 
5685  if (this->mHandledSBMLObjects.find(pSBase) == this->mHandledSBMLObjects.end())
5686  {
5687  removedObjects.push_back(pList->remove(i - 1));
5688  }
5689  }
5690 
5691  pList = pModel->getListOfParameters();
5692 
5693  for (i = pList->size(); i != 0; --i)
5694  {
5695  pSBase = pList->get(i - 1);
5696 
5697  // skip all specially marked items
5698  if (pSBase->getUserData() != NULL && strcmp((const char*)pSBase->getUserData(), "1") == 0)
5699  continue;
5700 
5701  if (this->mHandledSBMLObjects.find(pSBase) == this->mHandledSBMLObjects.end())
5702  {
5703  removedObjects.push_back(pList->remove(i - 1));
5704  }
5705  }
5706 
5707  pList = pModel->getListOfInitialAssignments();
5708 
5709  for (i = pList->size(); i != 0; --i)
5710  {
5711  pSBase = pList->get(i - 1);
5712 
5713  // skip all specially marked items
5714  if (pSBase->getUserData() != NULL && strcmp((const char*)pSBase->getUserData(), "1") == 0)
5715  continue;
5716 
5717  if (this->mHandledSBMLObjects.find(pSBase) == this->mHandledSBMLObjects.end())
5718  {
5719  // remove the item from the COPASI2SBMLmap
5720  pos = reverseCOPASI2SBMLMap.find(pSBase);
5721 
5722  if (pos != reverseCOPASI2SBMLMap.end() && pos->second != NULL)
5723  {
5724  pos2 = this->mCOPASI2SBMLMap.find(pos->second);
5725 
5726  if (pos2 != this->mCOPASI2SBMLMap.end())
5727  {
5728  this->mCOPASI2SBMLMap.erase(pos2);
5729  }
5730  }
5731 
5732  // we don't have to store the initial assignments since
5733  // there are no rules or events for initial assignments
5734  delete pList->remove(i - 1);
5735  }
5736  }
5737 
5738  pList = pModel->getListOfReactions();
5739 
5740  for (i = pList->size(); i != 0; --i)
5741  {
5742  pSBase = pList->get(i - 1);
5743 
5744  if (this->mHandledSBMLObjects.find(pSBase) == this->mHandledSBMLObjects.end())
5745  {
5746  // remove the item from the COPASI2SBMLmap
5747  pos = reverseCOPASI2SBMLMap.find(pSBase);
5748 
5749  if (pos != reverseCOPASI2SBMLMap.end() && pos->second != NULL)
5750  {
5751  pos2 = this->mCOPASI2SBMLMap.find(pos->second);
5752 
5753  if (pos2 != this->mCOPASI2SBMLMap.end())
5754  {
5755  this->mCOPASI2SBMLMap.erase(pos2);
5756  }
5757  }
5758 
5759  // we don't have to store the reactions since
5760  // there are no rules or events for reactions
5761  delete pList->remove(i - 1);
5762  }
5763  }
5764 
5765  // check each removed object, if there is a rule or an event in
5766  // the SBML model for this object
5767  std::list<SBase*>::iterator it2 = removedObjects.begin(), endit2 = removedObjects.end();
5768  pList = pModel->getListOfRules();
5769  ListOf* pEventList = pModel->getListOfEvents();
5770  ListOf* pInitialAssignmentList = pModel->getListOfInitialAssignments();
5771  Rule* pRule = NULL;
5772  Event* pEvent = NULL;
5773  InitialAssignment* pIA = NULL;
5774  unsigned int iMax;
5775 
5776  while (it2 != endit2)
5777  {
5778  iMax = pList->size();
5779 
5780  for (i = 0; i < iMax; ++i)
5781  {
5782  pRule = dynamic_cast<Rule*>(pList->get(i));
5783 
5784  if (pRule != NULL && pRule->isSetVariable() && pRule->getVariable() == (*it2)->getId())
5785  {
5786  // remove the item from the COPASI2SBMLmap
5787  pos = reverseCOPASI2SBMLMap.find(pRule);
5788 
5789  if (pos != reverseCOPASI2SBMLMap.end() && pos->second != NULL)
5790  {
5791  pos2 = this->mCOPASI2SBMLMap.find(pos->second);
5792 
5793  if (pos2 != this->mCOPASI2SBMLMap.end())
5794  {
5795  this->mCOPASI2SBMLMap.erase(pos2);
5796  }
5797  }
5798 
5799  delete pList->remove(i);
5800  break;
5801  }
5802  }
5803 
5804  iMax = pInitialAssignmentList->size();
5805 
5806  for (i = 0; i < iMax; ++i)
5807  {
5808  pIA = dynamic_cast<InitialAssignment*>(pInitialAssignmentList->get(i));
5809 
5810  if (pIA->isSetSymbol() && pIA->getSymbol() == (*it2)->getId())
5811  {
5812  // remove the item from the COPASI2SBMLmap
5813  pos = reverseCOPASI2SBMLMap.find(pIA);
5814 
5815  if (pos != reverseCOPASI2SBMLMap.end() && pos->second != NULL)
5816  {
5817  pos2 = this->mCOPASI2SBMLMap.find(pos->second);
5818 
5819  if (pos2 != this->mCOPASI2SBMLMap.end())
5820  {
5821  this->mCOPASI2SBMLMap.erase(pos2);
5822  }
5823  }
5824 
5825  delete pInitialAssignmentList->remove(i);
5826  break;
5827  }
5828  }
5829 
5830  iMax = pEventList->size();
5831 
5832  for (i = 0; i < iMax; ++i)
5833  {
5834  pEvent = dynamic_cast<Event*>(pEventList->get(i));
5835 
5836  if (pEvent != NULL)
5837  {
5838  ListOfEventAssignments* pEAL = pEvent->getListOfEventAssignments();
5839  EventAssignment* pEA = NULL;
5840  unsigned int j, jMax = pEAL->size();
5841 
5842  for (j = 0; j < jMax; ++j)
5843  {
5844  pEA = pEvent->getEventAssignment(j);
5845 
5846  if (pEA != NULL && pEA->getVariable() == (*it2)->getId())
5847  {
5848  // remove the item from the COPASI2SBMLmap
5849  pos = reverseCOPASI2SBMLMap.find(pEA);
5850 
5851  if (pos != reverseCOPASI2SBMLMap.end() && pos->second != NULL)
5852  {
5853  pos2 = this->mCOPASI2SBMLMap.find(pos->second);
5854 
5855  if (pos2 != this->mCOPASI2SBMLMap.end())
5856  {
5857  this->mCOPASI2SBMLMap.erase(pos2);
5858  }
5859  }
5860 
5861  delete pEAL->remove(j);
5862  break;
5863  }
5864  }
5865 
5866  // if there are no more event assignments, we can delete the
5867  // event
5868  if (pEAL->size() == 0)
5869  {
5870  // remove the item from the COPASI2SBMLmap
5871  pos = reverseCOPASI2SBMLMap.find(pEvent);
5872 
5873  if (pos != reverseCOPASI2SBMLMap.end() && pos->second != NULL)
5874  {
5875  pos2 = this->mCOPASI2SBMLMap.find(pos->second);
5876 
5877  if (pos2 != this->mCOPASI2SBMLMap.end())
5878  {
5879  this->mCOPASI2SBMLMap.erase(pos2);
5880  }
5881  }
5882 
5883  delete pEventList->remove(i);
5884  }
5885  }
5886  }
5887 
5888  // remove the item from the COPASI2SBMLmap
5889  pos = reverseCOPASI2SBMLMap.find(*it2);
5890 
5891  if (pos != reverseCOPASI2SBMLMap.end() && pos->second != NULL)
5892  {
5893  pos2 = this->mCOPASI2SBMLMap.find(pos->second);
5894 
5895  if (pos2 != this->mCOPASI2SBMLMap.end())
5896  {
5897  this->mCOPASI2SBMLMap.erase(pos2);
5898  }
5899  }
5900 
5901  delete *it2;
5902  ++it2;
5903  }
5904  }
5905 }
5906 
5907 CEvaluationNode* CSBMLExporter::createMassActionExpression(const std::vector<std::vector<std::string> >& arguments, bool isReversible)
5908 {
5909  assert((isReversible && arguments.size() == 4) || arguments.size() == 2);
5910  assert(arguments[0].size() == 1 && arguments[1].size() > 0);
5911  // create a multiplication of all items in arguments[1] and multiply that
5912  // with item arguments[0][0]
5913  std::set<std::string> finishedElements;
5914  std::vector<CEvaluationNode*> multiplicants;
5915  const CCopasiObject* pObject = CCopasiRootContainer::getKeyFactory()->get(arguments[0][0]);
5916  assert(pObject != NULL);
5917  multiplicants.push_back(new CEvaluationNodeObject(CEvaluationNodeObject::CN, "<" + pObject->getCN() + ",Reference=Value>"));
5918  std::vector<std::string>::const_iterator it = arguments[1].begin(), endit = arguments[1].end();
5919 
5920  while (it != endit)
5921  {
5922  if (finishedElements.find(*it) == finishedElements.end())
5923  {
5924 #ifdef __SUNPRO_CC
5925  // replaced count call with "old" count call
5926  // to compile on Sun Compiler
5927  unsigned int num = 0;
5928  std::count(arguments[1].begin(), arguments[1].end(), *it, num);
5929 #else
5930  size_t num = std::count(arguments[1].begin(), arguments[1].end(), *it);
5931 #endif // __SUNPRO_CC
5932  assert(num != 0);
5933  finishedElements.insert(*it);
5934  pObject = CCopasiRootContainer::getKeyFactory()->get(*it);
5935  assert(pObject != NULL);
5936 
5937  if (num == 1)
5938  {
5939  multiplicants.push_back(new CEvaluationNodeObject(CEvaluationNodeObject::CN, "<" + pObject->getCN() + ",Reference=Concentration>"));
5940  }
5941  else
5942  {
5943  std::ostringstream os;
5944  os << num;
5946  pOperator->addChild(new CEvaluationNodeObject(CEvaluationNodeObject::CN, "<" + pObject->getCN() + ",Reference=Concentration>"));
5948  multiplicants.push_back(pOperator);
5949  }
5950  }
5951 
5952  ++it;
5953  }
5954 
5955  std::vector<CEvaluationNode*>::reverse_iterator rIt = multiplicants.rbegin(), rEndit = multiplicants.rend();
5956  CEvaluationNode* pResult = *rIt;
5957  ++rIt;
5958  CEvaluationNode* pTmpNode = NULL;
5959 
5960  if (rIt != rEndit)
5961  {
5962  --rEndit;
5963 
5964  while (rIt != rEndit)
5965  {
5967  pTmpNode->addChild(*rIt);
5968  pTmpNode->addChild(pResult);
5969  pResult = pTmpNode;
5970  ++rIt;
5971  }
5972  }
5973 
5975  pTmpNode->addChild(*rIt);
5976  pTmpNode->addChild(pResult);
5977  pResult = pTmpNode;
5978 
5979  if (isReversible)
5980  {
5981  std::vector<std::vector<std::string> > tmpV;
5982  tmpV.push_back(arguments[2]);
5983  tmpV.push_back(arguments[3]);
5985  pTmpNode->addChild(pResult);
5986  pTmpNode->addChild(CSBMLExporter::createMassActionExpression(tmpV, false));
5987  pResult = pTmpNode;
5988  }
5989 
5990  return pResult;
5991 }
5992 
5993 bool CSBMLExporter::isValidSId(const std::string& id)
5994 {
5995  bool result = true;
5996 
5997  if (id.length() > 0)
5998  {
5999  char c = id[0];
6000  result = (c == '_' || (c > 64 && c < 91) || (c > 96 && c < 123));
6001  size_t i, iMax = id.length();
6002 
6003  for (i = 1; (i < iMax) && result; ++i)
6004  {
6005  c = id[i];
6006  result = (c == '_' || (c > 64 && c < 91) || (c > 96 && c < 123) || (c > 47 && c < 58));
6007  }
6008  }
6009  else
6010  {
6011  result = false;
6012  }
6013 
6014  return result;
6015 }
6016 
6017 /**
6018  * Remove the initial assignment for the entity with the given id
6019  * if there is any.
6020  */
6021 void CSBMLExporter::removeInitialAssignment(const std::string& sbmlId)
6022 {
6023  ListOfInitialAssignments* pList = this->mpSBMLDocument->getModel()->getListOfInitialAssignments();
6024  unsigned int i, iMax = pList->size();
6025 
6026  for (i = 0; i < iMax; ++i)
6027  {
6028  InitialAssignment* pIA = dynamic_cast<InitialAssignment*>(pList->get(i));
6029 
6030  if (pIA->getSymbol() == sbmlId)
6031  {
6032  pList->remove(i);
6033  break;
6034  }
6035  }
6036 }
6037 
6038 /**
6039  * Remove the rule for the entity with the given id
6040  * if there is any.
6041  */
6042 void CSBMLExporter::removeRule(const std::string& sbmlId)
6043 {
6044  ListOfRules* pList = this->mpSBMLDocument->getModel()->getListOfRules();
6045  unsigned int i, iMax = pList->size();
6046 
6047  for (i = 0; i < iMax; ++i)
6048  {
6049  Rule* pR = dynamic_cast<Rule*>(pList->get(i));
6050 
6051  if (pR->getVariable() == sbmlId)
6052  {
6053  pList->remove(i);
6054  break;
6055  }
6056  }
6057 }
6058 
6059 /**
6060  * Replaces references to species with reference to species divided by
6061  * volume if it is a reference to a concentration or by a reference to the
6062  * species times avogadros number if it is a reference to the amount.
6063  * The method also takes into consideration the substance units of the
6064  * model.
6065  */
6067 {
6068  CEvaluationNode* pResult = NULL;
6069  double factor = dataModel.getModel()->getQuantity2NumberFactor();
6070 
6072  {
6073  std::vector<CCopasiContainer*> containers;
6074  containers.push_back(const_cast<CModel*>(dataModel.getModel()));
6075  const CCopasiObject* pObject = dataModel.ObjectFromName(containers, dynamic_cast<const CEvaluationNodeObject*>(pOrigNode)->getObjectCN());
6076 
6077  if (pObject == NULL)
6078  {
6079  // this will be the object from the mInitialValueMap
6080  pResult = new CEvaluationNodeObject(CEvaluationNodeObject::CN, pOrigNode->getData());
6081  return pResult;
6082  }
6083 
6084  if (pObject->isReference())
6085  {
6086  const CCopasiObject* pParent = pObject->getObjectParent();
6087  // check if the parent is a metabolite
6088  const CMetab* pMetab = dynamic_cast<const CMetab*>(pParent);
6089 
6090  if (pMetab != NULL)
6091  {
6092  // check if the reference is to the concentration or to the
6093  // amount
6094  if (pObject->getObjectName() == "InitialConcentration" || pObject->getObjectName() == "Concentration")
6095  {
6096  // the concentration nodes only need to be replaced if we
6097  // are dealing with a model that has variable volumes
6098  if (this->mVariableVolumes == true)
6099  {
6100  // replace the node by the concentration node divided the the volume of
6101  // the (initial) volume of the node
6102  // although this is semantically incorrect, the result when
6103  // converted to SBML will be OK since the concentration ode
6104  // is converted to a reference to the species id which will
6105  // be interpreted as the amount due to the
6106  // hasOnlySubstanceUnits flag being set
6107  const CCompartment* pCompartment = pMetab->getCompartment();
6108  assert(pCompartment != NULL);
6110  // copy branch should be fine since the object node does
6111  // not have children
6112  pResult->addChild(pOrigNode->copyBranch());
6113 
6114  if (pObject->getObjectName() == "InitialConcentration")
6115  {
6116  pResult->addChild(new CEvaluationNodeObject(CEvaluationNodeObject::CN, "<" + pCompartment->getObject(CCopasiObjectName("Reference=InitialVolume"))->getCN() + ">"));
6117  }
6118  else
6119  {
6120  pResult->addChild(new CEvaluationNodeObject(CEvaluationNodeObject::CN, "<" + pCompartment->getObject(CCopasiObjectName("Reference=Volume"))->getCN() + ">"));
6121  }
6122  }
6123  else
6124  {
6125  // do nothing
6126  pResult = pOrigNode->copyBranch();
6127  }
6128  }
6129  else if (pObject->getObjectName() == "Rate")
6130  {
6131  std::string id = addRateOfIfItDoesNotExist(mpSBMLDocument, mIdMap, "rateOf");
6132  pResult = new CEvaluationNodeObject(CEvaluationNodeObject::INVALID, "<rateOf>");
6133  pResult->addChild(new CEvaluationNodeObject(CEvaluationNodeObject::CN, "<" + pObject->getObjectParent()->getObject(CCopasiObjectName("Reference=Concentration"))->getCN() + ">"));
6134  pResult->addChild(new CEvaluationNodeObject(CEvaluationNodeObject::CN, "<" + id + ">"));
6135  }
6136  else if (pObject->getObjectName() == "InitialParticleNumber" || pObject->getObjectName() == "ParticleNumber")
6137  {
6138  // if the units are not set to particle numbers anyway,
6139  // replace the node by the node times avogadros number
6140  if (dataModel.getModel()->getQuantityUnitEnum() != CModel::number)
6141  {
6142  if (this->mpAvogadro == NULL)
6143  {
6144  this->mpAvogadro = const_cast<CModel*>(dataModel.getModel())->createModelValue("quantity to number factor", dataModel.getModel()->getQuantity2NumberFactor());
6145  Parameter* pSBMLAvogadro = this->mpSBMLDocument->getModel()->createParameter();
6146  pSBMLAvogadro->setName("quantity to number factor");
6147  std::string sbmlId = CSBMLExporter::createUniqueId(this->mIdMap, mpAvogadro->getObjectName(), false);
6148  pSBMLAvogadro->setId(sbmlId);
6149  const_cast<CModelValue*>(this->mpAvogadro)->setSBMLId(sbmlId);
6150  this->mIdMap.insert(std::pair<const std::string, const SBase*>(sbmlId, pSBMLAvogadro));
6151 
6152  if (this->mSBMLLevel != 1)
6153  {
6154  pSBMLAvogadro->setConstant(true);
6155  }
6156  else
6157  {
6158  // Level 1 doesn't know the constant flag and
6159  // libSBML does not drop it automatically
6160  pSBMLAvogadro->setConstant(true);
6161  }
6162 
6163  pSBMLAvogadro->setValue(dataModel.getModel()->getQuantity2NumberFactor());
6164  this->mHandledSBMLObjects.insert(pSBMLAvogadro);
6165  this->mCOPASI2SBMLMap[this->mpAvogadro] = pSBMLAvogadro;
6166  this->mAvogadroCreated = true;
6167  }
6168 
6170  // copyBranch should be save here since object nodes can't
6171  // have children
6172  pResult->addChild(pOrigNode->copyBranch());
6173  pResult->addChild(new CEvaluationNodeObject(CEvaluationNodeObject::CN, "<" + this->mpAvogadro->getCN() + ",Reference=InitialValue>"));
6174  }
6175  else
6176  {
6177  // the result is the same as the original node
6178  pResult = pOrigNode->copyBranch();
6179  }
6180 
6181  if (this->mVariableVolumes == false)
6182  {
6183  // multiply by the volume as well
6184  const CCompartment* pCompartment = pMetab->getCompartment();
6185 
6186  if (pCompartment->getDimensionality() != 0)
6187  {
6189  pTmpNode->addChild(pResult);
6190 
6191  if (pObject->getObjectName() == "InitialParticleNumber")
6192  {
6193  pTmpNode->addChild(new CEvaluationNodeObject(CEvaluationNodeObject::CN, "<" + pCompartment->getObject(CCopasiObjectName("Reference=InitialVolume"))->getCN() + ">"));
6194  }
6195  else
6196  {
6197  pTmpNode->addChild(new CEvaluationNodeObject(CEvaluationNodeObject::CN, "<" + pCompartment->getObject(CCopasiObjectName("Reference=Volume"))->getCN() + ">"));
6198  }
6199 
6200  pResult = pTmpNode;
6201  }
6202  }
6203  }
6204  else
6205  {
6206  fatalError();
6207  }
6208  }
6209  }
6210  }
6211  // check if there is a division by avogadros number and if so, just
6212  // drop the division instead of introducing a new multiplication
6214  {
6215  // check if one of the child nodes is a reference to a species
6216 
6217  const CEvaluationNode* pLeft = dynamic_cast<const CEvaluationNode*>(pOrigNode->getChild());
6218  const CEvaluationNode* pRight = dynamic_cast<const CEvaluationNode*>(pLeft->getSibling());
6219  assert(pLeft != NULL);
6220  assert(pRight != NULL);
6221 
6223  {
6224  std::vector<CCopasiContainer*> containers;
6225  containers.push_back(const_cast<CModel*>(dataModel.getModel()));
6226  const CCopasiObject* pObject = dataModel.ObjectFromName(containers, dynamic_cast<const CEvaluationNodeObject*>(pLeft)->getObjectCN());
6227 
6228  if (pObject != NULL && pObject->isReference())
6229  {
6230  const CCopasiObject* pParent = pObject->getObjectParent();
6231  // check if the parent is a metabolite
6232  const CMetab* pMetab = dynamic_cast<const CMetab*>(pParent);
6233 
6234  if (pMetab != NULL)
6235  {
6236  // check if pRight is a number or a parameter that
6237  // corresponds to avogadros number
6241  {
6242  double value = dynamic_cast<const CEvaluationNodeNumber*>(pRight)->getValue();
6243 
6244  if (fabs((factor - value) / factor) <= 1e-3)
6245  {
6246  // copyBranch should be OK since the node has no
6247  // children anyway
6248  pResult = pLeft->copyBranch();
6249  }
6250  }
6251  else if (CEvaluationNode::type(pRight->getType()) == CEvaluationNode::OBJECT)
6252  {
6253  const CCopasiObject* pObject2 = dataModel.ObjectFromName(containers, dynamic_cast<const CEvaluationNodeObject*>(pRight)->getObjectCN());
6254 
6255  if (pObject2 != NULL && pObject2->isReference())
6256  {
6257  const CCopasiObject* pObjectParent2 = pObject2->getObjectParent();
6258  const CModelValue* pMV = dynamic_cast<const CModelValue*>(pObjectParent2);
6259 
6260  if (pMV != NULL && pMV->getStatus() == CM