SBMLImporter.cpp

Go to the documentation of this file.

// Copyright (C) 2010 - 2015 by Pedro Mendes, Virginia Tech Intellectual
// Properties, Inc., University of Heidelberg, and The University
// of Manchester.
// All rights reserved.

// Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
// Properties, Inc., EML Research, gGmbH, University of Heidelberg,
// and The University of Manchester.
// All rights reserved.

// Copyright (C) 2004 - 2007 by Pedro Mendes, Virginia Tech Intellectual
// Properties, Inc. and EML Research, gGmbH.
// All rights reserved.

#ifdef WIN32
# pragma warning (disable: 4786)
# pragma warning (disable: 4243)
// warning C4355: 'this' : used in base member initializer list
# pragma warning (disable: 4355)
#endif  // WIN32

#define USE_LAYOUT 1

#include <iostream>
#include <vector>
#include <sstream>
#include <map>
#include <limits>
#include <cmath>

#include <sbml/SBMLReader.h>
#include <sbml/SBMLDocument.h>
#include <sbml/Compartment.h>
#include <sbml/Species.h>
#include <sbml/SpeciesReference.h>
#include <sbml/Reaction.h>
#include <sbml/LocalParameter.h>

#if LIBSBML_VERSION >= 50400
#include <sbml/SBMLTransforms.h>
#include "IdList.h" // this file is missing from the libSBML distribution!
#include <sbml/conversion/ConversionProperties.h>
#include <sbml/packages/layout/extension/LayoutModelPlugin.h>
#define INIT_DEFAULTS(element) \
  {\
    element.initDefaults();\
  }
#else
#define INIT_DEFAULTS(element) \
  {\
  }
#endif // LIBSBML_VERSION

#include <sbml/KineticLaw.h>
#include <sbml/math/FormulaFormatter.h>
#include <sbml/Model.h>
#include <sbml/UnitKind.h>
#include <sbml/Unit.h>
#include <sbml/Parameter.h>
#include <sbml/InitialAssignment.h>
#include <sbml/Rule.h>
#include <sbml/FunctionDefinition.h>
#include <sbml/UnitDefinition.h>

#include "copasi.h"

#include "report/CKeyFactory.h"
#include "model/CModel.h"
#include "model/CCompartment.h"
#include "model/CMetab.h"
#include "model/CReaction.h"
#include "model/CModelValue.h"
#include "model/CEvent.h"
#include "function/CNodeK.h"
#include "function/CFunctionDB.h"
#include "function/CEvaluationTree.h"
#include "function/CExpression.h"
#include "function/CFunctionParameters.h"
#include "report/CCopasiObjectReference.h"
#include "utilities/CCopasiTree.h"
#include "utilities/CNodeIterator.h"
#include "CopasiDataModel/CCopasiDataModel.h"
#include "report/CCopasiRootContainer.h"
#include "MIRIAM/CRDFGraphConverter.h"
#include "compareExpressions/CEvaluationNodeNormalizer.h"
#include "commandline/CLocaleString.h"
#include "commandline/COptions.h"

#include "SBMLImporter.h"
#include "SBMLUtils.h"
#include "ConverterASTNode.h"
#include "utilities/CProcessReport.h"
#include "copasi/commandline/CConfigurationFile.h"

#include "layout/SBMLDocumentLoader.h"
#include "layout/CListOfLayouts.h"

#include "utilities/CCopasiMessage.h"

#if LIBSBML_VERSION >= 50903 && LIBSBML_HAS_PACKAGE_COMP

#include <sbml/util/PrefixTransformer.h>
#include <sbml/packages/comp/extension/CompModelPlugin.h>

class CPrefixNameTransformer : public PrefixTransformer
{
public:
  CPrefixNameTransformer() {}

  static void replaceStringInPlace(std::string& subject, const std::string& search,
                                   const std::string& replace)
  {
    size_t pos = 0;

    while ((pos = subject.find(search, pos)) != std::string::npos)
      {
        subject.replace(pos, search.length(), replace);
        pos += replace.length();
      }
  }

  static inline std::string &rtrim(std::string &str)
  {
    size_t endpos = str.find_last_not_of(" \t");

    if (std::string::npos != endpos)
      {
        str = str.substr(0, endpos + 1);
      }

    return str;
  }

  const std::string& cleanName(std::string& prefix)
  {
    std::replace(prefix.begin(), prefix.end(), '_', ' ');
    replaceStringInPlace(prefix, "  ", " ");
    rtrim(prefix);
    return prefix;
  }

  virtual int transform(SBase* element)
  {
    if (element == NULL || getPrefix().empty())
      return LIBSBML_OPERATION_SUCCESS;

    // set up ids
    PrefixTransformer::transform(element);

    // skip local parameters, as they are not renamed
    if (element->getTypeCode() == SBML_LOCAL_PARAMETER)
      return LIBSBML_OPERATION_SUCCESS;

    // setup names
    if (element->isSetName())
      {
        std::stringstream newName;
        std::string prefix = getPrefix();
        newName << element->getName() << " (" << cleanName(prefix) << ")";
        element->setName(newName.str());
      }

    return LIBSBML_OPERATION_SUCCESS;
  }
};

#endif

// static
C_FLOAT64 SBMLImporter::round(const C_FLOAT64 & x)
{
  return
    x < 0.0 ? -floor(-x + 0.5) : floor(x + 0.5);
}

// static
bool SBMLImporter::areApproximatelyEqual(const double & x, const double & y, const double & t)
{
  double Scale =
    (fabs(x) + fabs(y)) * t;

  // Avoid underflow
  if (Scale < 100.0 * std::numeric_limits< C_FLOAT64 >::min())
    return true;

  return 2 * fabs(x - y) < Scale;
}

std::string getOriginalSBMLId(Parameter* parameter)
{
  if (parameter == NULL) return "";

  if (!parameter->isSetAnnotation()) return "";

  XMLNode* node = parameter->getAnnotation();

  if (node->getNumChildren() < 1) return "";

  for (unsigned int i = 0; i < node->getNumChildren(); ++i)
    {
      const XMLNode& current = node->getChild(i);

      if (current.getNamespaces().containsUri("http://copasi.org/initialValue"))
        {
          return current.getAttrValue("parent");
        }
    }

  return "";
}

std::string getInitialCNForSBase(SBase* sbase, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
{
  std::map<CCopasiObject*, SBase*>::const_iterator it;

  for (it = copasi2sbmlmap.begin(); it != copasi2sbmlmap.end(); ++it)
    {
      if (it->second != sbase)
        continue;

      CMetab *metab = dynamic_cast<CMetab*>(it->first);

      if (metab != NULL)
        return metab->getInitialConcentrationReference()->getCN();

      CCompartment *comp = dynamic_cast<CCompartment*>(it->first);

      if (comp != NULL)
        return comp->getInitialValueReference()->getCN();

      CModelValue *param = dynamic_cast<CModelValue*>(it->first);

      if (param != NULL)
        return param->getInitialValueReference()->getCN();
    }

  return "";
}

/**
 * Creates and returns a Copasi CModel from the SBMLDocument given as argument.
 */
CModel* SBMLImporter::createCModelFromSBMLDocument(SBMLDocument* sbmlDocument, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
{
  Model* sbmlModel = sbmlDocument->getModel();

  /* Create an empty model and set the title. */
  this->mpCopasiModel = new CModel(mpDataModel);
  copasi2sbmlmap[this->mpCopasiModel] = sbmlModel;
  this->mpCopasiModel->setLengthUnit(CModel::m);
  this->mpCopasiModel->setAreaUnit(CModel::m2);
  this->mpCopasiModel->setVolumeUnit(CModel::l);
  this->mpCopasiModel->setTimeUnit(CModel::s);
  this->mpCopasiModel->setQuantityUnit(CModel::Mol);
  this->mpCopasiModel->setSBMLId(sbmlModel->getId());

  unsigned C_INT32 step = 0, totalSteps = 0;
  size_t hStep = C_INVALID_INDEX;

  mImportStep = 1;

  if (mpImportHandler && !mpImportHandler->progressItem(mhImportStep)) return NULL;

  SBMLImporter::importMIRIAM(sbmlModel, this->mpCopasiModel);
  UnitDefinition *pSubstanceUnits = NULL;
  UnitDefinition *pTimeUnits = NULL;
  UnitDefinition *pVolumeUnits = NULL;
  UnitDefinition *pAreaUnits = NULL;
  UnitDefinition *pLengthUnits = NULL;
  /* Set standard units to match the standard units of SBML files. */

  // for SBML L3 files the default units are defined on the model
  if (this->mLevel > 2)
    {
      this->mpCopasiModel->setAvogadro(6.02214179e23);
      this->mAvogadroSet = true;

      // we make copies of the unit definitions so that we do not have to remember
      // if we created them or not
      std::string units;
      UnitDefinition* pUDef = NULL;
      Unit unit(sbmlModel->getLevel(), sbmlModel->getVersion());
      INIT_DEFAULTS(unit);

      if (sbmlModel->isSetSubstanceUnits())
        {
          units = sbmlModel->getSubstanceUnits();
          assert(units != "");
          pUDef =  sbmlModel->getUnitDefinition(units);

          if (pUDef != NULL)
            {
              pSubstanceUnits = new UnitDefinition(*pUDef);
            }
          else
            {
              if (units == "mole")
                {
                  unit.setKind(UNIT_KIND_MOLE);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pSubstanceUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pSubstanceUnits->addUnit(&unit);
                }
              else if (units == "item")
                {
                  unit.setKind(UNIT_KIND_ITEM);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pSubstanceUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pSubstanceUnits->addUnit(&unit);
                }
              else if (units == "dimensionless")
                {
                  unit.setKind(UNIT_KIND_DIMENSIONLESS);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pSubstanceUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pSubstanceUnits->addUnit(&unit);
                }
              else
                {
                  std::string message = "COPASI can't handle substance unit \"" + units + "\". Setting unit for substances to dimensionless.";
                  CCopasiMessage(CCopasiMessage::WARNING, message.c_str());
                  unit.setKind(UNIT_KIND_DIMENSIONLESS);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pSubstanceUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pSubstanceUnits->addUnit(&unit);
                }
            }
        }

      if (sbmlModel->isSetTimeUnits())
        {
          units = sbmlModel->getTimeUnits();
          assert(units != "");
          pUDef =  sbmlModel->getUnitDefinition(units);

          if (pUDef != NULL)
            {
              pTimeUnits = new UnitDefinition(*pUDef);
            }
          else
            {
              if (units == "second")
                {
                  unit.setKind(UNIT_KIND_SECOND);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pTimeUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pTimeUnits->addUnit(&unit);
                }
              else if (units == "dimensionless")
                {
                  unit.setKind(UNIT_KIND_DIMENSIONLESS);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pTimeUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pTimeUnits->addUnit(&unit);
                }
              else
                {
                  std::string message = "COPASI can't handle time unit \"" + units + "\". Setting unit for time to dimensionless.";
                  CCopasiMessage(CCopasiMessage::WARNING, message.c_str());
                  unit.setKind(UNIT_KIND_DIMENSIONLESS);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pTimeUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pTimeUnits->addUnit(&unit);
                }
            }
        }

      if (sbmlModel->isSetVolumeUnits())
        {
          units = sbmlModel->getVolumeUnits();
          assert(units != "");
          pUDef =  sbmlModel->getUnitDefinition(units);

          if (pUDef != NULL)
            {
              pVolumeUnits = new UnitDefinition(*pUDef);
            }
          else
            {
              if (units == "litre")
                {
                  unit.setKind(UNIT_KIND_LITRE);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pVolumeUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pVolumeUnits->addUnit(&unit);
                }
              else if (units == "dimensionless")
                {
                  unit.setKind(UNIT_KIND_DIMENSIONLESS);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pVolumeUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pVolumeUnits->addUnit(&unit);
                }
              else
                {
                  std::string message = "COPASI can't handle volume unit \"" + units + "\". Setting unit for volume to dimensionless.";
                  CCopasiMessage(CCopasiMessage::WARNING, message.c_str());
                  unit.setKind(UNIT_KIND_DIMENSIONLESS);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pVolumeUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pVolumeUnits->addUnit(&unit);
                }
            }
        }

      if (sbmlModel->isSetAreaUnits())
        {
          units = sbmlModel->getAreaUnits();
          assert(units != "");
          pUDef =  sbmlModel->getUnitDefinition(units);

          if (pUDef != NULL)
            {
              pAreaUnits = new UnitDefinition(*pUDef);
            }
          else
            {
              if (units == "dimensionless")
                {
                  unit.setKind(UNIT_KIND_DIMENSIONLESS);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pAreaUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pAreaUnits->addUnit(&unit);
                }
              else
                {
                  std::string message = "COPASI can't handle area unit \"" + units + "\". Setting unit for area to dimensionless.";
                  CCopasiMessage(CCopasiMessage::WARNING, message.c_str());
                  unit.setKind(UNIT_KIND_DIMENSIONLESS);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pAreaUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pAreaUnits->addUnit(&unit);
                }
            }
        }

      if (sbmlModel->isSetLengthUnits())
        {
          units = sbmlModel->getLengthUnits();
          assert(units != "");
          pUDef =  sbmlModel->getUnitDefinition(units);

          if (pUDef != NULL)
            {
              pLengthUnits = new UnitDefinition(*pUDef);
            }
          else
            {
              if (units == "litre")
                {
                  unit.setKind(UNIT_KIND_LITRE);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pLengthUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pLengthUnits->addUnit(&unit);
                }
              else if (units == "metre")
                {
                  unit.setKind(UNIT_KIND_METRE);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pLengthUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pLengthUnits->addUnit(&unit);
                }
              else if (units == "dimensionless")
                {
                  unit.setKind(UNIT_KIND_DIMENSIONLESS);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pLengthUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pLengthUnits->addUnit(&unit);
                }
              else
                {
                  std::string message = "COPASI can't handle length unit \"" + units + "\". Setting unit for length to dimensionless.";
                  CCopasiMessage(CCopasiMessage::WARNING, message.c_str());
                  unit.setKind(UNIT_KIND_DIMENSIONLESS);
                  unit.setExponent(1);
                  unit.setMultiplier(1.0);
                  unit.setScale(0);
                  pLengthUnits = new UnitDefinition(sbmlModel->getLevel(), sbmlModel->getVersion());
                  pLengthUnits->addUnit(&unit);
                }
            }
        }
    }
  else
    {
      if (sbmlModel->getNumUnitDefinitions() != 0)
        {
          unsigned int counter;

          // we make copies o the unit definitions so that we do not have to remember
          // if we created them or not
          for (counter = 0; counter < sbmlModel->getNumUnitDefinitions(); counter++)
            {
              UnitDefinition* uDef = sbmlModel->getUnitDefinition(counter);
              std::string unitId = uDef->getId();

              if (unitId == "substance")
                {
                  pSubstanceUnits = new UnitDefinition(*uDef);
                }
              else if (unitId == "time")
                {
                  pTimeUnits = new UnitDefinition(*uDef);
                }
              else if (unitId == "volume")
                {
                  pVolumeUnits = new UnitDefinition(*uDef);
                }
              else if ((unitId == "area"))
                {
                  pAreaUnits = new UnitDefinition(*uDef);
                }
              else if ((unitId == "length"))
                {
                  pLengthUnits = new UnitDefinition(*uDef);
                }
            }
        }
    }

  // create the default units if some unit has not been specified
  if (pSubstanceUnits == NULL)
    {
      if (this->mLevel > 2)
        {
          // issue a warning
          CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 91, "substance", "mole" , "substance");
        }

      // create the default units
      pSubstanceUnits = new UnitDefinition(this->mLevel, this->mVersion);
      pSubstanceUnits->setId("dummy_substance");
      pSubstanceUnits->setName("dummy_substance");
      Unit* pUnit = pSubstanceUnits->createUnit();
      pUnit->setKind(UNIT_KIND_MOLE);
      pUnit->setExponent(1);
      pUnit->setMultiplier(1.0);
      pUnit->setScale(0);
    }

  if (pTimeUnits == NULL)
    {
      if (this->mLevel > 2)
        {
          // issue a warning
          CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 91, "time", "second" , "time");
        }

      // create the default units
      pTimeUnits = new UnitDefinition(this->mLevel, this->mVersion);
      pTimeUnits->setId("dummy_time");
      pTimeUnits->setName("dummy_time");
      Unit* pUnit = pTimeUnits->createUnit();
      pUnit->setKind(UNIT_KIND_SECOND);
      pUnit->setExponent(1);
      pUnit->setMultiplier(1.0);
      pUnit->setScale(0);
    }

  if (pVolumeUnits == NULL)
    {
      if (this->mLevel > 2)
        {
          // issue a warning
          CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 91, "volume", "litre" , "volume");
        }

      // create the default units
      pVolumeUnits = new UnitDefinition(this->mLevel, this->mVersion);
      pVolumeUnits->setId("dummy_volume");
      pVolumeUnits->setName("dummy_volume");
      Unit* pUnit = pVolumeUnits->createUnit();
      pUnit->setKind(UNIT_KIND_LITRE);
      pUnit->setExponent(1);
      pUnit->setMultiplier(1.0);
      pUnit->setScale(0);
    }

  if (pAreaUnits == NULL)
    {
      if (this->mLevel > 2)
        {
          // issue a warning
          CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 91, "area", "m^2" , "area");
        }

      // create the default units
      pAreaUnits = new UnitDefinition(this->mLevel, this->mVersion);
      pAreaUnits->setId("dummy_area");
      pAreaUnits->setName("dummy_area");
      Unit* pUnit = pAreaUnits->createUnit();
      pUnit->setKind(UNIT_KIND_METRE);
      pUnit->setExponent(2);
      pUnit->setMultiplier(1.0);
      pUnit->setScale(0);
    }

  if (pLengthUnits == NULL)
    {
      if (this->mLevel > 2)
        {
          // issue a warning
          CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 91, "length", "m" , "length");
        }

      // create the default units
      pLengthUnits = new UnitDefinition(this->mLevel, this->mVersion);
      pLengthUnits->setId("dummy_length");
      pLengthUnits->setName("dummy_length");
      Unit* pUnit = pLengthUnits->createUnit();
      pUnit->setKind(UNIT_KIND_METRE);
      pUnit->setExponent(1);
      pUnit->setMultiplier(1.0);
      pUnit->setScale(0);
    }

  // now we have some common code to actually import the units

  // handle the substance units
  assert(pSubstanceUnits != NULL);

  if (pSubstanceUnits != NULL)
    {
      std::pair<CModel::QuantityUnit, bool> qUnit;

      try
        {
          qUnit = this->handleSubstanceUnit(pSubstanceUnits);
        }
      catch (...)
        {
          std::ostringstream os;
          os << "Error while importing substance units.";

          // check if the last message on the stack is an exception
          // and if so, add the message text to the current exception
          if (CCopasiMessage::peekLastMessage().getType() == CCopasiMessage::EXCEPTION)
            {
              // we only want the message, not the timestamp line
              std::string text = CCopasiMessage::peekLastMessage().getText();
              os << text.substr(text.find("\n"));
            }

          CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
        }

      if (qUnit.second == false)
        {
          // the unit could not be handled, give an error message and
          // set the units to mole
          CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 66, "substance", "Mole");
          this->mpCopasiModel->setQuantityUnit(CModel::Mol);
        }
      else
        {
          this->mpCopasiModel->setQuantityUnit(qUnit.first);
        }

      // check if the extends units are set and if they are equal to the substance units
      // otherwise issue a warning
      if (this->mLevel > 2)
        {
          if (!sbmlModel->isSetExtentUnits())
            {
              CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 92);
            }
          else
            {
              const UnitDefinition* pExtendsUnits = sbmlModel->getUnitDefinition(sbmlModel->getExtentUnits());

              if (pExtendsUnits != NULL)
                {
                  if (!areSBMLUnitDefinitionsIdentical(pSubstanceUnits, pExtendsUnits))
                    {
                      CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 92);
                    }
                }
              else
                {
                  CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 66, "extends", "the same units as the substances");
                }
            }
        }

      delete pSubstanceUnits;
      pSubstanceUnits = NULL;
    }

  // handle the time units
  assert(pTimeUnits != NULL);

  if (pTimeUnits != NULL)
    {
      std::pair<CModel::TimeUnit, bool> tUnit;

      try
        {
          tUnit = this->handleTimeUnit(pTimeUnits);
        }
      catch (...)
        {
          std::ostringstream os;
          os << "Error while importing time units.";

          // check if the last message on the stack is an exception
          // and if so, add the message text to the current exception
          if (CCopasiMessage::peekLastMessage().getType() == CCopasiMessage::EXCEPTION)
            {
              // we only want the message, not the timestamp line
              std::string text = CCopasiMessage::peekLastMessage().getText();
              os << text.substr(text.find("\n"));
            }

          CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
        }

      if (tUnit.second == false)
        {
          // the unit could not be handled, give an error message and
          // set the units to second
          CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 66, "time", "second");
          this->mpCopasiModel->setTimeUnit(CModel::s);
        }
      else
        {
          this->mpCopasiModel->setTimeUnit(tUnit.first);
        }

      delete pTimeUnits;
      pTimeUnits = NULL;
    }

  // handle the volume units
  assert(pVolumeUnits != NULL);

  if (pVolumeUnits != NULL)
    {
      std::pair<CModel::VolumeUnit, bool> vUnit;

      try
        {
          vUnit = this->handleVolumeUnit(pVolumeUnits);
        }
      catch (...)
        {
          std::ostringstream os;
          os << "Error while importing volume units.";

          // check if the last message on the stack is an exception
          // and if so, add the message text to the current exception
          if (CCopasiMessage::peekLastMessage().getType() == CCopasiMessage::EXCEPTION)
            {
              // we only want the message, not the timestamp line
              std::string text = CCopasiMessage::peekLastMessage().getText();
              os << text.substr(text.find("\n"));
            }

          CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
        }

      if (vUnit.second == false)
        {
          // the unit could not be handled, give an error message and
          // set the units to litre
          CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 66, "volume", "litre");
          this->mpCopasiModel->setVolumeUnit(CModel::l);
        }
      else
        {
          this->mpCopasiModel->setVolumeUnit(vUnit.first);
        }

      delete pVolumeUnits;
      pVolumeUnits = NULL;
    }

  // handle the area units
  assert(pAreaUnits != NULL);

  if (pAreaUnits != NULL)
    {
      std::pair<CModel::AreaUnit, bool> vUnit;

      try
        {
          vUnit = this->handleAreaUnit(pAreaUnits);
        }
      catch (...)
        {
          std::ostringstream os;
          os << "Error while importing area units.";

          // check if the last message on the stack is an exception
          // and if so, add the message text to the current exception
          if (CCopasiMessage::peekLastMessage().getType() == CCopasiMessage::EXCEPTION)
            {
              // we only want the message, not the timestamp line
              std::string text = CCopasiMessage::peekLastMessage().getText();
              os << text.substr(text.find("\n"));
            }

          CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
        }

      if (vUnit.second == false)
        {
          // the unit could not be handled, give an error message and
          // set the units to litre
          CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 66, "area", "square meter");
          this->mpCopasiModel->setAreaUnit(CModel::m2);
        }
      else
        {
          this->mpCopasiModel->setAreaUnit(vUnit.first);
        }

      delete pAreaUnits;
      pAreaUnits = NULL;
    }

  // handle the length units
  assert(pLengthUnits != NULL);

  if (pLengthUnits != NULL)
    {
      std::pair<CModel::LengthUnit, bool> vUnit;

      try
        {
          vUnit = this->handleLengthUnit(pLengthUnits);
        }
      catch (...)
        {
          std::ostringstream os;
          os << "Error while importing length units.";

          // check if the last message on the stack is an exception
          // and if so, add the message text to the current exception
          if (CCopasiMessage::peekLastMessage().getType() == CCopasiMessage::EXCEPTION)
            {
              // we only want the message, not the timestamp line
              std::string text = CCopasiMessage::peekLastMessage().getText();
              os << text.substr(text.find("\n"));
            }

          CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
        }

      if (vUnit.second == false)
        {
          // the unit could not be handled, give an error message and
          // set the units to litre
          CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 66, "length", "meter");
          this->mpCopasiModel->setLengthUnit(CModel::m);
        }
      else
        {
          this->mpCopasiModel->setLengthUnit(vUnit.first);
        }

      delete pLengthUnits;
      pLengthUnits = NULL;
    }

  // go through all compartments and species and check if the units are
  // consistent
  checkElementUnits(sbmlModel, this->mpCopasiModel, this->mLevel, this->mVersion);
  std::string title;

  if (this->isStochasticModel(sbmlModel))
    {
      this->mpCopasiModel->setModelType(CModel::stochastic);
    }
  else
    {
      this->mpCopasiModel->setModelType(CModel::deterministic);
    }

  SBMLImporter::importNotes(this->mpCopasiModel, sbmlModel);

  title = sbmlModel->getName();

  if (title == "")
    {
      title = "NoName";
    }

  size_t idCount = 0;

  while (mpDataModel->getObject("Model=" + title) != NULL)
    {
      std::stringstream str; str << sbmlModel->getName() << "_" << ++idCount;
      title = str.str();
    }

  this->mpCopasiModel->setObjectName(title.c_str());
  // fill the set of SBML species reference ids because
  // we need this to check for references to species references in all expressions
  // as long as we do not support these references
  SBMLImporter::updateSBMLSpeciesReferenceIds(sbmlModel, this->mSBMLSpeciesReferenceIds);

  /* import the functions */
  unsigned int counter;
  CCopasiVectorN< CFunction >* functions = &(this->functionDB->loadedFunctions());

  size_t num = (*functions).size();

  if (mpImportHandler)
    {
      mImportStep = 2;

      if (!mpImportHandler->progressItem(mhImportStep)) return NULL;

      totalSteps = (unsigned C_INT32) num;
      hStep = mpImportHandler->addItem("Importing function definitions",
                                       step,
                                       &totalSteps);
    }

  this->sbmlIdMap.clear();

  for (counter = 0; counter < num; ++counter)
    {
      CFunction* tree = (*functions)[counter];

      if (!tree->getSBMLId().empty())
        {
          this->sbmlIdMap[tree] = tree->getSBMLId();
          tree->setSBMLId("");
        }
    }

  CFunctionDB* pTmpFunctionDB = this->importFunctionDefinitions(sbmlModel, copasi2sbmlmap);
  // try to find global parameters that represent avogadros number
  this->findAvogadroConstant(sbmlModel, this->mpCopasiModel->getQuantity2NumberFactor());

  std::map<std::string, CCompartment*> compartmentMap;

  /* Create the compartments */
  num = sbmlModel->getNumCompartments();

  if (mpImportHandler)
    {
      mpImportHandler->finishItem(hStep);
      mImportStep = 3;

      if (!mpImportHandler->progressItem(mhImportStep)) return NULL;

      step = 0;
      totalSteps = (unsigned C_INT32) num;
      hStep = mpImportHandler->addItem("Importing compartments...",
                                       step,
                                       &totalSteps);
    }

  for (counter = 0; counter < num; counter++)
    {
      Compartment* sbmlCompartment = sbmlModel->getCompartment(counter);

      if (sbmlCompartment == NULL)
        {
          fatalError();
        }

      CCompartment* pCopasiCompartment = NULL;

      try
        {
          pCopasiCompartment = this->createCCompartmentFromCompartment(sbmlCompartment, this->mpCopasiModel, copasi2sbmlmap, sbmlModel);
        }
      catch (...)
        {
          std::ostringstream os;
          os << "Error while importing compartment \"";
          os << sbmlCompartment->getId() << "\".";

          // check if the last message on the stack is an exception
          // and if so, add the message text to the current exception
          if (CCopasiMessage::peekLastMessage().getType() == CCopasiMessage::EXCEPTION)
            {
              // we only want the message, not the timestamp line
              std::string text = CCopasiMessage::peekLastMessage().getText();
              os << text.substr(text.find("\n"));
            }

          CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
        }

      std::string key = sbmlCompartment->getId();

      if (mLevel == 1)
        {
          key = sbmlCompartment->getName();
        }

      compartmentMap[key] = pCopasiCompartment;
      ++step;

      if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return NULL;
    }

  /* Create all species */
  num = sbmlModel->getNumSpecies();

  if (mpImportHandler)
    {
      mpImportHandler->finishItem(hStep);
      mImportStep = 4;

      if (!mpImportHandler->progressItem(mhImportStep)) return NULL;

      step = 0;
      totalSteps = (unsigned C_INT32) num;
      hStep = mpImportHandler->addItem("Importing species...",
                                       step,
                                       &totalSteps);
    }

  for (counter = 0; counter < num; ++counter)
    {
      Species* sbmlSpecies = sbmlModel->getSpecies(counter);

      if (sbmlSpecies == NULL)
        {
          fatalError();
        }

      CCompartment* pCopasiCompartment = compartmentMap[sbmlSpecies->getCompartment()];

      if (pCopasiCompartment != NULL)
        {
          CMetab* pCopasiMetabolite = NULL;

          try
            {
              pCopasiMetabolite = this->createCMetabFromSpecies(sbmlSpecies, this->mpCopasiModel, pCopasiCompartment, copasi2sbmlmap, sbmlModel);
            }
          catch (...)
            {
              std::ostringstream os;
              os << "Error while importing species \"";
              os << sbmlSpecies->getId() << "\".";

              // check if the last message on the stack is an exception
              // and if so, add the message text to the current exception
              if (CCopasiMessage::peekLastMessage().getType() == CCopasiMessage::EXCEPTION)
                {
                  // we only want the message, not the timestamp line
                  std::string text = CCopasiMessage::peekLastMessage().getText();
                  os << text.substr(text.find("\n"));
                }

              CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
            }

          std::string key;

          if (this->mLevel == 1)
            {
              key = sbmlSpecies->getName();
            }
          else
            {
              key = sbmlSpecies->getId();
            }

          this->speciesMap[key] = pCopasiMetabolite;
        }
      else
        {
          CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 5 , sbmlSpecies->getCompartment().c_str(), sbmlSpecies->getId().c_str());
        }

      ++step;

      if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return NULL;
    }

  /* Create the global Parameters */
  num = sbmlModel->getNumParameters();

  if (mpImportHandler)
    {
      mpImportHandler->finishItem(hStep);
      mImportStep = 5;

      if (!mpImportHandler->progressItem(mhImportStep)) return NULL;

      step = 0;
      totalSteps = (unsigned C_INT32) num;
      hStep = mpImportHandler->addItem("Importing global parameters...",
                                       step,
                                       &totalSteps);
    }

  for (counter = 0; counter < num; ++counter)
    {
      Parameter* sbmlParameter = sbmlModel->getParameter(counter);

      if (sbmlParameter == NULL)
        {
          fatalError();
        }

      try
        {
          std::string sbmlId = getOriginalSBMLId(sbmlParameter);

          // don't import parameter if it is one of our special ones instead get the cn to the target
          if (!sbmlId.empty())
            {
              Species* species = sbmlModel->getSpecies(sbmlId);

              if (species != NULL)
                {
                  mKnownInitalValues[sbmlParameter->getId()] = getInitialCNForSBase(species, copasi2sbmlmap);
                  continue;
                }

              Compartment *comp = sbmlModel->getCompartment(sbmlId);

              if (comp != NULL)
                {
                  mKnownInitalValues[sbmlParameter->getId()] = getInitialCNForSBase(comp, copasi2sbmlmap);
                  continue;
                }

              Parameter* param = sbmlModel->getParameter(sbmlId);

              if (param != NULL)
                {
                  mKnownInitalValues[sbmlParameter->getId()] = getInitialCNForSBase(param, copasi2sbmlmap);
                  continue;
                }
            }

          this->createCModelValueFromParameter(sbmlParameter, this->mpCopasiModel, copasi2sbmlmap);
        }
      catch (...)
        {
          std::ostringstream os;
          os << "Error while importing parameter \"";
          os << sbmlParameter->getId() << "\".";

          // check if the last message on the stack is an exception
          // and if so, add the message text to the current exception
          if (CCopasiMessage::peekLastMessage().getType() == CCopasiMessage::EXCEPTION)
            {
              // we only want the message, not the timestamp line
              std::string text = CCopasiMessage::peekLastMessage().getText();
              os << text.substr(text.find("\n"));
            }

          CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
        }

      ++step;

      if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return NULL;
    }

  // the information that is collected here is used in the creation of the reactions to
  // adjust the stoichiometry of substrates and products with the conversion factos from
  // the SBML model

  // now that all parameters have been imported, we can check if the model
  // defines a global conversion factor
  if (this->mLevel > 2 && sbmlModel->isSetConversionFactor())
    {
      this->mConversionFactorFound = true;
      std::string id = sbmlModel->getConversionFactor();
      assert(id != "");
      std::map<std::string, const CModelValue*>::const_iterator pos = this->mSBMLIdModelValueMap.find(id);
      assert(pos != this->mSBMLIdModelValueMap.end());

      if (pos != this->mSBMLIdModelValueMap.end())
        {
          this->mpModelConversionFactor = pos->second;
        }
      else
        {
          fatalError();
        }
    }

  // we need to go through all species again and find conversion factors
  // right now I don't want to change the order of importing model values
  // and species since I am not sure what implications that would have
  if (this->mLevel > 2)
    {
      num = sbmlModel->getNumSpecies();

      for (counter = 0; counter < num; ++counter)
        {
          Species* sbmlSpecies = sbmlModel->getSpecies(counter);

          if (sbmlSpecies->isSetConversionFactor())
            {
              this->mConversionFactorFound = true;
              std::map<std::string, const CModelValue*>::const_iterator pos = this->mSBMLIdModelValueMap.find(sbmlSpecies->getConversionFactor());

              if (pos != this->mSBMLIdModelValueMap.end())
                {
                  this->mSpeciesConversionParameterMap[sbmlSpecies->getId()] = pos->second;
                }
            }
        }
    }

  if (!this->mIgnoredParameterUnits.empty())
    {
      std::ostringstream os;
      std::vector<std::string>::iterator errorIt = this->mIgnoredParameterUnits.begin();

      while (errorIt != this->mIgnoredParameterUnits.end())
        {
          os << *errorIt << ", ";
          ++errorIt;
        }

      std::string s = os.str();
      CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 26, s.substr(0, s.size() - 2).c_str());
    }

  /* Create all reactions */
  num = sbmlModel->getNumReactions();

  if (mpImportHandler)
    {
      mpImportHandler->finishItem(hStep);
      mImportStep = 6;

      if (!mpImportHandler->progressItem(mhImportStep)) return NULL;

      step = 0;
      totalSteps = (unsigned C_INT32) num;
      hStep = mpImportHandler->addItem("Importing reactions...",
                                       step,
                                       &totalSteps);
    }

  this->mDivisionByCompartmentReactions.clear();
  this->mFastReactions.clear();
  this->mReactionsWithReplacedLocalParameters.clear();

  for (counter = 0; counter < num; counter++)
    {
      try
        {
          this->createCReactionFromReaction(sbmlModel->getReaction(counter), sbmlModel, this->mpCopasiModel, copasi2sbmlmap, pTmpFunctionDB);
        }
      catch (...)
        {
          std::ostringstream os;
          os << "Error while importing reaction \"";
          os << sbmlModel->getReaction(counter)->getId() << "\".";

          // check if the last message on the stack is an exception
          // and if so, add the message text to the current exception
          if (CCopasiMessage::peekLastMessage().getType() == CCopasiMessage::EXCEPTION)
            {
              // we only want the message, not the timestamp line
              std::string text = CCopasiMessage::peekLastMessage().getText();
              os << text.substr(text.find("\n"));
            }

          CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
        }

      ++step;

      if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return NULL;
    }

  if (!this->mDivisionByCompartmentReactions.empty())
    {
      // create the error message
      std::string idList;
      std::set<std::string>::const_iterator it = this->mDivisionByCompartmentReactions.begin();
      std::set<std::string>::const_iterator endit = this->mDivisionByCompartmentReactions.end();

      while (it != endit)
        {
          idList += (*it);
          idList += ", ";
          ++it;
        }

      idList = idList.substr(0, idList.length() - 2);
      CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 17, idList.c_str());
    }

  if (!this->mFastReactions.empty())
    {
      // create the error message
      std::string idList;
      std::set<std::string>::const_iterator it = this->mFastReactions.begin();
      std::set<std::string>::const_iterator endit = this->mFastReactions.end();

      while (it != endit)
        {
          idList += (*it);
          idList += ", ";
          ++it;
        }

      idList = idList.substr(0, idList.length() - 2);
      CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 29, idList.c_str());
    }

  if (!this->mReactionsWithReplacedLocalParameters.empty())
    {
      // create the error message
      std::string idList;
      std::set<std::string>::const_iterator it = this->mReactionsWithReplacedLocalParameters.begin();
      std::set<std::string>::const_iterator endit = this->mReactionsWithReplacedLocalParameters.end();

      while (it != endit)
        {
          idList += (*it);
          idList += ", ";
          ++it;
        }

      idList = idList.substr(0, idList.length() - 2);
      CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 87, idList.c_str());
    }

  // import the initial assignments
  // we do this after the reactions since intial assignments can reference reaction ids.
  if (mpImportHandler)
    {
      mpImportHandler->finishItem(hStep);
      mImportStep = 7;

      if (!mpImportHandler->progressItem(mhImportStep)) return NULL;
    }

  importInitialAssignments(sbmlModel, copasi2sbmlmap, this->mpCopasiModel);

  // import all rules
  // we have to import them after the reactions since they can reference a reaction
  // id which has to be known when importing the rule

  /* Create the rules */
  // rules should be imported after reactions because we use the mSBMLSpeciesReferenceIds to determine if a
  // rule changes a species reference (stoichiometry
  // Since COPASI does not support this, we need to ignore the rule
  this->areRulesUnique(sbmlModel);
  num = sbmlModel->getNumRules();

  if (mpImportHandler)
    {
      mpImportHandler->finishItem(hStep);
      mImportStep = 8;

      if (!mpImportHandler->progressItem(mhImportStep)) return NULL;

      step = 0;
      totalSteps = (unsigned C_INT32) num;
      hStep = mpImportHandler->addItem("Importing global parameters...",
                                       step,
                                       &totalSteps);
    }

  for (counter = 0; counter < num; ++counter)
    {
      Rule* sbmlRule = sbmlModel->getRule(counter);

      if (sbmlRule == NULL)
        {
          fatalError();
        }

      // improve the error message a bit.
      try
        {
          this->importSBMLRule(sbmlRule, copasi2sbmlmap, sbmlModel);
        }
      catch (...)
        {
          std::ostringstream os;
          os << "Error while importing rule for variable \"";
          os << sbmlRule->getVariable() << "\".";

          // check if the last message on the stack is an exception
          // and if so, add the message text to the current exception
          if (CCopasiMessage::peekLastMessage().getType() == CCopasiMessage::EXCEPTION)
            {
              // we only want the message, not the time stamp line
              std::string text = CCopasiMessage::peekLastMessage().getText();
              os << text.substr(text.find("\n"));
            }

          CCopasiMessage(CCopasiMessage::EXCEPTION, os.str().c_str());
        }

      ++step;

      if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return NULL;
    }

  if (sbmlModel->getNumConstraints() > 0)
    {
      CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 49);
    }

  // TODO Create all constraints
  // TODO Since we don't have constraints yet, there is no code here.
  // TODO When implementing import of constraints, don't forget to replace calls to
  // TODO explicitly time dependent functions in the constraints math expression.

  // import all events
  // events should be imported after reactions because we use the mSBMLSpeciesReferenceIds to determine if an
  // event assignment changes a species reference (stoichiometry
  // Since COPASI does not support this, we need to ignore the event assignment
  if (mpImportHandler)
    {
      mpImportHandler->finishItem(hStep);
      mImportStep = 9;

      if (!mpImportHandler->progressItem(mhImportStep)) return NULL;
    }

  this->importEvents(sbmlModel, this->mpCopasiModel, copasi2sbmlmap);

  this->mpCopasiModel->setCompileFlag();

  if (this->mUnsupportedRuleFound)
    {
      CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 3);
    }

  if (this->mRateRuleForSpeciesReferenceIgnored == true)
    {
      CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 94 , "Rate rule" , "Rate rule");
    }

  if (this->mEventAssignmentForSpeciesReferenceIgnored == true)
    {
      CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 94 , "Event assignment", "event assignment");
    }

  if (this->mConversionFactorFound == true)
    {
      CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 100);
    }

  if (this->mEventPrioritiesIgnored == true)
    {
      CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 98);
    }

  if (this->mInitialTriggerValues == true)
    {
      CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 97);
    }

  if (this->mNonPersistentTriggerFound == true)
    {
      CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 99);
    }

  if (mpImportHandler)
    {
      mpImportHandler->finishItem(hStep);
      mImportStep = 10;

      if (!mpImportHandler->progressItem(mhImportStep)) return NULL;
    }

  // unset the hasOnlySubstanceUnits flag on all such species
  std::map<Species*, Compartment*>::iterator it = this->mSubstanceOnlySpecies.begin();
  std::map<Species*, Compartment*>::iterator endIt = this->mSubstanceOnlySpecies.end();

  while (it != endIt)
    {
      it->first->setHasOnlySubstanceUnits(false);
      ++it;
    }

  setInitialValues(this->mpCopasiModel, copasi2sbmlmap);
  // evaluate and apply the initial expressions
  this->applyStoichiometricExpressions(copasi2sbmlmap, sbmlModel);

  // now we apply the conversion factors
  if (this->mLevel > 2)
    {
      this->applyConversionFactors();
    }

  this->removeUnusedFunctions(pTmpFunctionDB, copasi2sbmlmap);

  // remove the temporary avogadro parameter if one was created
  if (this->mAvogadroCreated == true)
    {
      const Parameter* pParameter = *this->mPotentialAvogadroNumbers.begin();
      ListOf* pList = sbmlModel->getListOfParameters();
      unsigned i, iMax = pList->size();

      for (i = 0; i < iMax; ++i)
        {
          if (pList->get(i)->getId() == pParameter->getId())
            {
              pList->remove(i);
              break;
            }
        }
    }

  delete pTmpFunctionDB;

  // create a warning if the delay function is used in the model
  if (this->mDelayFound)
    {
      CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 36);
    }

  // give a warning that units on pure number as allowed in SBML L3 and above
  // are ignored by COPASI
  if (this->mLevel > 2 && this->mUnitOnNumberFound)
    {
      CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 93);
    }

  this->mpCopasiModel->forceCompile(this->mpImportHandler);
  return this->mpCopasiModel;
}

std::string isKnownCustomFunctionDefinition(const FunctionDefinition* sbmlFunction,
    const std::string& sNamespace,
    const std::string& elementName,
    const std::string& definition)
{
  FunctionDefinition* current = const_cast<FunctionDefinition*>(sbmlFunction);

  if (current == NULL) return "";

  if (!current->isSetAnnotation()) return "";

  const XMLNode* element = current->getAnnotation();

  if (element == NULL) return "";

  for (unsigned int i = 0 ; i < element->getNumChildren(); ++i)
    {
      const XMLNode& annot = element->getChild(i);

      if (annot.getURI() == sNamespace &&
          annot.getName() == elementName &&
          annot.getAttrValue("definition") == definition)
        {
          return current->getId();
        }
    }

  return "";
}

bool addToKnownFunctionToMap(std::map<std::string, std::string>& map, const FunctionDefinition* sbmlFunction)
{
  if (!sbmlFunction->isSetAnnotation())
    return false;

  std::string id = isKnownCustomFunctionDefinition(sbmlFunction,
                   "http://sbml.org/annotations/symbols",
                   "symbols",
                   "http://en.wikipedia.org/wiki/Derivative");

  if (!id.empty())
    {
      map[id] = "RATE";
      return true;
    }

  id = isKnownCustomFunctionDefinition(sbmlFunction,
                                       "http://sbml.org/annotations/distribution",
                                       "distribution",
                                       "http://www.uncertml.org/distributions/normal");

  if (!id.empty())
    {
      map[id] = "RNORMAL";
      return true;
    }

  id = isKnownCustomFunctionDefinition(sbmlFunction,
                                       "http://sbml.org/annotations/distribution",
                                       "distribution",
                                       "http://www.uncertml.org/distributions/uniform");

  if (!id.empty())
    {
      map[id] = "RUNIFORM";
      return true;
    }

  id = isKnownCustomFunctionDefinition(sbmlFunction,
                                       "http://sbml.org/annotations/distribution",
                                       "distribution",
                                       "http://www.uncertml.org/distributions/gamma");

  if (!id.empty())
    {
      map[id] = "RGAMMA";
      return true;
    }

  id = isKnownCustomFunctionDefinition(sbmlFunction,
                                       "http://sbml.org/annotations/distribution",
                                       "distribution",
                                       "http://www.uncertml.org/distributions/poisson");

  if (!id.empty())
    {
      map[id] = "RPOISSON";
      return true;
    }

  return false;
}

int
AstStrCmp(const void *s1, const void *s2)
{
  const char* a = static_cast<const ASTNode *>(s1)->getName();
  const char* b = static_cast<const ASTNode *>(s2)->getName();

  if (a == NULL && b == NULL) return  0;

  if (a == NULL && b != NULL) return -1;

  if (a != NULL && b == NULL) return  1;

  return strcmp(
           a,
           b);
}

/**
 * This function checks the function definition for unused arguments (that would not
 * be properly displayed in the CopasiUI). If found, the function body will be replaced
 * with one including all arguments.
 */
void ensureAllArgsAreBeingUsedInFunctionDefinition(const FunctionDefinition* sbmlFunction)
{
  if (sbmlFunction == NULL || sbmlFunction->getNumArguments() == 0 || sbmlFunction->getBody() == NULL) return;

  // get all variables
  List *variables = sbmlFunction->getBody()->getListOfNodes(ASTNode_isName);

  // find unused ones
  std::vector<std::string> unused;

  for (unsigned int i = 0; i < sbmlFunction->getNumArguments(); ++i)
    {
      const ASTNode *arg = sbmlFunction->getArgument(i);

      if (variables->find(arg, AstStrCmp) == NULL)
        {
          if (arg->getName() != NULL)
            unused.push_back(arg->getName());
        }
    }

  // get rid of the list
  delete variables;

  // let us hope this is empty
  if (unused.size() == 0)
    return;

  // it is not, so modify the function definition to include all of them
  std::stringstream str;
  str << "lambda(";

  for (unsigned int i = 0; i < sbmlFunction->getNumArguments(); ++i)
    str << sbmlFunction->getArgument(i)->getName() << ", ";

  char* formula = SBML_formulaToString(sbmlFunction->getBody());
  str << formula;

  std::vector<std::string>::iterator it;

  for (it = unused.begin(); it != unused.end(); ++it)
    str << " + 0*" << *it;

  str << ")";

  // update the function definition
  const_cast<FunctionDefinition*>(sbmlFunction)->setMath(SBML_parseFormula(str.str().c_str()));

  // free the formula
  free(formula);
}

CFunction* SBMLImporter::createCFunctionFromFunctionDefinition(const FunctionDefinition* sbmlFunction, CFunctionDB* pTmpFunctionDB, Model* pSBMLModel, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
{

  ensureAllArgsAreBeingUsedInFunctionDefinition(sbmlFunction);

  addToKnownFunctionToMap(mKnownCustomUserDefinedFunctions, sbmlFunction);

  CFunction* pTmpFunction = this->createCFunctionFromFunctionTree(sbmlFunction, pSBMLModel, copasi2sbmlmap);

  if (pTmpFunction == NULL)
    {
      CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 14, sbmlFunction->getId().c_str());
      return NULL;
    }

  std::string sbmlId = sbmlFunction->getId();
  pTmpFunction->setSBMLId(sbmlId);
  // check if the id is already taken by another function definition, maybe
  // from an earlier import, if this is the case, delete the id on the old
  // function definition
  // if we don't do this, two functions might have the same SBML id during
  // export which makes the exporter code so much more difficult
  size_t i, iMax = this->functionDB->loadedFunctions().size();

  for (i = 0; i < iMax; ++i)
    {
      CFunction* pFun = this->functionDB->loadedFunctions()[i];

      if (pFun->getSBMLId() == sbmlId)
        {
          pFun->setSBMLId("");
        }
    }

  std::string functionName = sbmlFunction->getName();

  if (functionName == "")
    {
      functionName = sbmlFunction->getId();
    }

  unsigned int counter = 1;
  std::ostringstream numberStream;
  std::string appendix = "";
  CFunction * pExistingFunction = NULL;

  while ((pExistingFunction = functionDB->findFunction(functionName + appendix)))
    {
      if (areEqualFunctions(pExistingFunction, pTmpFunction))
        {
          pdelete(pTmpFunction);
          pTmpFunction = pExistingFunction;

          break;
        }

      // We need to check whether the functions are identical.
      numberStream.str("");
      numberStream << "_" << counter;
      counter++;
      appendix = numberStream.str();
    }

  if (pTmpFunction != pExistingFunction)
    {
      pTmpFunction->setObjectName(functionName + appendix);
      functionDB->add(pTmpFunction, true);
      pTmpFunctionDB->add(pTmpFunction, false);
    }

  if (pTmpFunction->getType() == CEvaluationTree::UserDefined)
    {
      SBMLImporter::importMIRIAM(sbmlFunction, pTmpFunction);
      SBMLImporter::importNotes(pTmpFunction, sbmlFunction);
    }

  return pTmpFunction;
}

CFunction* SBMLImporter::createCFunctionFromFunctionTree(const FunctionDefinition* pSBMLFunction, Model* pSBMLModel, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
{
  CFunction* pFun = NULL;

  if (!pSBMLFunction->isSetMath())
    return NULL;

  ConverterASTNode root(*pSBMLFunction->getMath());

  if (SBMLImporter::isDelayFunctionUsed(&root))
    {
      CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 85, pSBMLFunction->getId().c_str());
    }

  this->preprocessNode(&root, pSBMLModel, copasi2sbmlmap);

  if (root.getType() != AST_LAMBDA)
    {
      CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 11, pSBMLFunction->getId().c_str());
      return NULL;
    }

  // get the number of children.
  // the first n-1 children are the parameters for the function
  // the last child is the actual function
  pFun = new CKinFunction();
  unsigned int i, iMax = root.getNumChildren() - 1;
  std::set<std::string> variableNames;

  for (i = 0; i < iMax; ++i)
    {
      ASTNode* pVarNode = root.getChild(i);

      if (pVarNode->getType() != AST_NAME)
        {
          delete pFun;
          CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 12, pSBMLFunction->getId().c_str());
        }

      pFun->addVariable(pVarNode->getName());
      variableNames.insert(pVarNode->getName());
    }

  // now we check if the AST tree has a node that represents the time
  // object
  // find a unique name for the time variable
  std::ostringstream sstream;
  std::string timeVariableName = "time";
  unsigned int postfix = 1;

  while (variableNames.find(timeVariableName) != variableNames.end())
    {
      sstream.str("");
      sstream << "time_" << postfix;
      timeVariableName = sstream.str();
      ++postfix;
    }

  if (this->replaceTimeNodesInFunctionDefinition(root.getChild(iMax), timeVariableName))
    {
      // add another variable to the function
      ASTNode* pVarNode = new ASTNode(AST_NAME);
      pVarNode->setName(timeVariableName.c_str());
      ASTNode* pTmpNode = root.removeChild(iMax);
      root.addChild(pVarNode);
      root.addChild(pTmpNode);
      // increase iMax since we now have one more child
      ++iMax;
      pFun->addVariable(timeVariableName);
      this->mExplicitelyTimeDependentFunctionDefinitions.insert(pSBMLFunction->getId());
    }

  pFun->setTree(*root.getChild(iMax));
  CCopasiTree<CEvaluationNode>::iterator treeIt = pFun->getRoot();

  // if the root node already is an object node, this has to be dealt with separately
  if (dynamic_cast<CEvaluationNodeObject*>(&(*treeIt)))
    {
      CEvaluationNodeVariable* pVariableNode = new CEvaluationNodeVariable(CEvaluationNodeVariable::ANY, (*treeIt).getData().substr(1, (*treeIt).getData().length() - 2));
      pFun->setRoot(pVariableNode);
    }
  else
    {
      while (treeIt != NULL)
        {
          if (dynamic_cast<CEvaluationNodeObject*>(&(*treeIt)))
            {
              CEvaluationNodeVariable* pVariableNode = new CEvaluationNodeVariable(CEvaluationNodeVariable::ANY, (*treeIt).getData().substr(1, (*treeIt).getData().length() - 2));

              if ((*treeIt).getParent())
                {
                  (*treeIt).getParent()->addChild(pVariableNode, &(*treeIt));
                  (*treeIt).getParent()->removeChild(&(*treeIt));
                }

              delete &(*treeIt);
              treeIt = pVariableNode;
            }

          ++treeIt;
        }
    }

  pFun->updateTree();

  if (!pFun->compile())
    {
      delete pFun;
      CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 28, pSBMLFunction->getId().c_str());
    }

  if (pFun->getRoot() == NULL)
    {
      delete pFun;
      CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 13, pSBMLFunction->getId().c_str());
    }

  return pFun;
}

CFunction* getFunctionForKey(CCopasiVectorN<CFunction> &functionDb, const std::string& key)
{
  CFunction* pFunc = NULL;
  CCopasiVectorN<CFunction>::iterator it = functionDb.begin(), endit = functionDb.end();

  while (it != endit)
    {
      if ((*it)->getKey() == key)
        {
          pFunc = dynamic_cast<CFunction*>(*it);
          break;
        }

      ++it;
    }

  return pFunc;
}

/**
 * Creates and returns a COPASI CCompartment from the SBML Compartment
 * given as argument.
 */
CCompartment*
SBMLImporter::createCCompartmentFromCompartment(const Compartment* sbmlCompartment, CModel* copasiModel, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap, const Model* /*pSBMLModel*/)
{
  if (sbmlCompartment->isSetUnits())
    {
      std::string cU = sbmlCompartment->getUnits();
    }

  unsigned int dimensionality = 3;

  if (sbmlCompartment->isSetSpatialDimensions())
    {
      dimensionality = sbmlCompartment->getSpatialDimensions();

      // starting with SBML Level 3, the spatial dimensions of a compartment can be
      // any rational number
      double dDimensionality = sbmlCompartment->getSpatialDimensions();

      if (this->mLevel > 2)
        {
          dDimensionality = sbmlCompartment->getSpatialDimensionsAsDouble();
        }

      // check if the unsigned int dimensionality corresponds to the double
      // dimensionality, otherwise give an error message
      // Actually if we wanted to be correct, we would have to check if the
      // part before the komma also matches and maybe have a separate error message if not
      dDimensionality -= dimensionality;

      if (fabs(dDimensionality) > 1e-9)
        {
          CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 89, sbmlCompartment->getId().c_str());
          dimensionality = 3;
        }
    }
  else
    {
      CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 90, sbmlCompartment->getId().c_str());
      dimensionality = 3;
    }

  if (dimensionality > 3)
    {
      CCopasiMessage Message(CCopasiMessage::WARNING,
                             "Reaction with id \"%s\" has dimensions of %d, this is not supported by COPASI. COPASI will assume that the compartment is three dimensional."
                             , sbmlCompartment->getId().c_str(), dimensionality);
      dimensionality = 3;
      //fatalError();
    }

  std::string name = sbmlCompartment->getName();

  if (name == "")
    {
      name = sbmlCompartment->getId();
    }

  std::string appendix = "";
  unsigned int counter = 2;
  std::ostringstream numberStream;

  while (copasiModel->getCompartments().getIndex(name + appendix) != C_INVALID_INDEX)
    {
      numberStream.str("");
      numberStream << "_" << counter;
      counter++;
      appendix = numberStream.str();
    }

  double value;
  CCompartment* copasiCompartment = copasiModel->createCompartment(name + appendix, value);

  if (this->mLevel == 1)
    {
      copasiCompartment->setSBMLId(sbmlCompartment->getName());
    }
  else
    {
      copasiCompartment->setSBMLId(sbmlCompartment->getId());
    }

  // set the dimension of the compartment
  copasiCompartment->setDimensionality(dimensionality);

  //DebugFile << "Created Compartment: " << copasiCompartment->getObjectName() << std::endl;
  SBMLImporter::importMIRIAM(sbmlCompartment, copasiCompartment);
  SBMLImporter::importNotes(copasiCompartment, sbmlCompartment);
  copasi2sbmlmap[copasiCompartment] = const_cast<Compartment*>(sbmlCompartment);
  return copasiCompartment;
}

/**
 * Creates and returns a Copasi CMetab from the given SBML Species object.
 */
CMetab*
SBMLImporter::createCMetabFromSpecies(const Species* sbmlSpecies, CModel* copasiModel, CCompartment* copasiCompartment, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap, const Model* /*pSBMLModel*/)
{
  if (sbmlSpecies->isSetSubstanceUnits())
    {
      std::string cU = sbmlSpecies->getSubstanceUnits();
    }

  std::map<CCopasiObject*, SBase*>::iterator it = copasi2sbmlmap.find(copasiCompartment);

  if (it == copasi2sbmlmap.end())
    {
      fatalError();
    }

  Compartment* pSBMLCompartment = (Compartment*)it->second;

  if (sbmlSpecies->getHasOnlySubstanceUnits() == true)
    {
      this->mSubstanceOnlySpecies.insert(std::make_pair(const_cast<Species*>(sbmlSpecies), pSBMLCompartment));
    }

  std::string name = sbmlSpecies->getName();

  if (name == "")
    {
      name = sbmlSpecies->getId();
    }

  std::string appendix = "";
  unsigned int counter = 2;
  std::ostringstream numberStream;

  while (copasiCompartment->getMetabolites().getIndex(name + appendix) != C_INVALID_INDEX)
    {
      numberStream.str("");
      numberStream << "_" << counter;
      counter++;
      appendix = numberStream.str();
    }

  CMetab* copasiMetabolite = copasiModel->createMetabolite(name + appendix, copasiCompartment->getObjectName());

  if (copasiMetabolite == NULL)
    {
      //DebugFile << "Could not create Copasi species." << std::endl;
      fatalError();
    }

  if (sbmlSpecies->getConstant() || sbmlSpecies->getBoundaryCondition())
    {
      copasiMetabolite->setStatus(CModelEntity::FIXED);
    }
  else
    {
      copasiMetabolite->setStatus(CModelEntity::REACTIONS);
    }

  // also check if the compartment has a spatialSize of 0 because this also implies hasOnlySubstanceUnits for the species in this compartment
  // In Level 3 files this no longer seems to be the case however, so we only
  // do this for L1 and L2 files
  if (pSBMLCompartment->isSetSpatialDimensions() && pSBMLCompartment->getSpatialDimensions() == 0 && this->mLevel < 3)
    {
      this->mSubstanceOnlySpecies.insert(std::make_pair(const_cast<Species*>(sbmlSpecies), pSBMLCompartment));
    }

  //DebugFile << "Created species: " << copasiMetabolite->getObjectName() << std::endl;
  copasi2sbmlmap[copasiMetabolite] = const_cast<Species*>(sbmlSpecies);

  if (this->mLevel == 1)
    {
      copasiMetabolite->setSBMLId(sbmlSpecies->getName());
    }
  else
    {
      copasiMetabolite->setSBMLId(sbmlSpecies->getId());
    }

  SBMLImporter::importMIRIAM(sbmlSpecies, copasiMetabolite);
  SBMLImporter::importNotes(copasiMetabolite, sbmlSpecies);
  return copasiMetabolite;
}

CFunction* findFunction(CCopasiVectorN < CFunction > &  db, const CFunction* func)
{
  size_t i, iMax = db.size();

  for (i = 0; i < iMax; ++i)
    {
      CFunction* pFun = (db[i]);

      if (*pFun == *func)
        return pFun;
    }

  return NULL;
}

/**
 * Creates and returns a COPASI CReaction object from the given SBML
 * Reaction object.
 */
CReaction*
SBMLImporter::createCReactionFromReaction(Reaction* sbmlReaction, Model* pSBMLModel, CModel* copasiModel, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap, CFunctionDB* pTmpFunctionDB)
{
  if (sbmlReaction == NULL)
    {
      fatalError();
    }

  std::string name = sbmlReaction->getName();

  if (name == "")
    {
      name = sbmlReaction->getId();
    }

  /* Check if the name of the reaction is unique. */
  std::string appendix = "";
  unsigned int counter = 2;
  std::ostringstream numberStream;

  while (copasiModel->getReactions().getIndex(name + appendix) != C_INVALID_INDEX)
    {
      numberStream.str("");
      numberStream << "_" << counter;
      counter++;
      appendix = numberStream.str();
    }

  /* create a new reaction with the unique name */
  CReaction* copasiReaction = copasiModel->createReaction(name + appendix);

  if (copasiReaction == NULL)
    {
      fatalError();
    }

  copasiReaction->setReversible(sbmlReaction->getReversible());

  if (this->mLevel == 1)
    {
      copasiReaction->setSBMLId(sbmlReaction->getName());
    }
  else
    {
      copasiReaction->setSBMLId(sbmlReaction->getId());
    }

  copasi2sbmlmap[copasiReaction] = const_cast<Reaction*>(sbmlReaction);

  if (sbmlReaction->isSetFast() && sbmlReaction->getFast() == true)
    {
      const_cast<Reaction*>(sbmlReaction)->setFast(false);
      this->mFastReactions.insert(sbmlReaction->getId());
    }

  // store if the reaction involves species that need a conversion factor
  bool mConversionFactorNeeded = false;

  /* Add all substrates to the reaction */
  unsigned int num = sbmlReaction->getNumReactants();
  bool singleCompartment = true;
  const CCompartment* compartment = NULL;
  bool hasOnlySubstanceUnitPresent = false;
  bool ignoreMassAction = false;

  for (counter = 0; counter < num; counter++)
    {
      const SpeciesReference* sr = sbmlReaction->getReactant(counter);

      if (sr == NULL)
        {
          delete copasiReaction;
          fatalError();
        }

      C_FLOAT64 stoi = 1.0;

      if (this->mLevel < 3)
        {
          // We may not use Mass Action kinetics if we do not know the actual stoichiometry.
          if (sr->isSetStoichiometryMath())
            {
              ignoreMassAction = true;
            }
          else
            {
              stoi = sr->getStoichiometry() / sr->getDenominator();
            }
        }
      else
        {
          // We may not use Mass Action kinetics if we do not know the actual stoichiometry.
          if (sr->isSetId() &&
              (pSBMLModel->getInitialAssignment(sr->getId()) ||
               pSBMLModel->getRule(sr->getId())))
            {
              ignoreMassAction = true;
            }

          if (sr->isSetStoichiometry())
            {
              stoi = sr->getStoichiometry();
            }
        }

      std::map<std::string, CMetab*>::iterator pos;
      pos = this->speciesMap.find(sr->getSpecies());

      // check if we have a conversion factor
      if (pos == this->speciesMap.end())
        {
          delete copasiReaction;
          fatalError();
        }
      else
        {
          // check if there is a conversion factor on the species
          if (this->mpModelConversionFactor != NULL || this->mSpeciesConversionParameterMap.find(pos->second->getSBMLId()) != this->mSpeciesConversionParameterMap.end())
            {
              mConversionFactorNeeded = true;
            }
        }

      std::map<CCopasiObject*, SBase*>::const_iterator spos = copasi2sbmlmap.find(pos->second);
      assert(spos != copasi2sbmlmap.end());
      Species* pSBMLSpecies = dynamic_cast<Species*>(spos->second);
      assert(pSBMLSpecies != NULL);
      hasOnlySubstanceUnitPresent = (hasOnlySubstanceUnitPresent | (pSBMLSpecies->getHasOnlySubstanceUnits() == true));

      if (compartment == NULL)
        {
          compartment = pos->second->getCompartment();
        }
      else
        {
          if (singleCompartment && compartment != pos->second->getCompartment())
            {
              singleCompartment = false;
            }
        }

      copasiReaction->addSubstrate(pos->second->getKey(), stoi);

      // we need to store the id of the species reference if it is set because SBML Level 3 allows
      // references to species references and if we want to support his, we need the id to import
      // expressions that reference a species reference
      if (this->mLevel > 2)
        {
          CCopasiVector<CChemEqElement>::const_iterator it = copasiReaction->getChemEq().getSubstrates().begin(), endit = copasiReaction->getChemEq().getSubstrates().end();
          CChemEqElement* pElement = NULL;

          while (it != endit)
            {
              if ((*it)->getMetaboliteKey() == pos->second->getKey())
                {
                  pElement = const_cast<CChemEqElement*>(*it);
                  break;
                }

              ++it;
            }

          assert(pElement != NULL);

          if (pElement != NULL)
            {
              copasi2sbmlmap[pElement] = const_cast<SpeciesReference*>(sr);

              this->mChemEqElementSpeciesIdMap[pElement] = std::pair<std::string, CChemEq::MetaboliteRole>(sr->getSpecies(), CChemEq::SUBSTRATE);
            }
        }

      // find the CChemEqElement that belongs to the added substrate
      if (this->mLevel < 3 && sr->isSetStoichiometryMath())
        {
          CChemEq& chemEq = copasiReaction->getChemEq();
          CCopasiVector < CChemEqElement >::const_iterator it = chemEq.getSubstrates().begin();
          CCopasiVector < CChemEqElement >::const_iterator end = chemEq.getSubstrates().end();
          CChemEqElement* pChemEqElement = NULL;

          while (it != end)
            {
              if ((*it)->getMetabolite() == pos->second)
                {
                  pChemEqElement = (*it);
                  break;
                }

              ++it;
            }

          assert(pChemEqElement != NULL);
          mStoichiometricExpressionMap.insert(std::make_pair(sr->getStoichiometryMath()->getMath(), pChemEqElement));
        }
    }

  /* Add all products to the reaction */
  num = sbmlReaction->getNumProducts();

  for (counter = 0; counter < num; counter++)
    {
      const SpeciesReference* sr = sbmlReaction->getProduct(counter);

      if (sr == NULL)
        {
          delete copasiReaction;
          fatalError();
        }

      C_FLOAT64 stoi = 1.0;

      if (this->mLevel < 3)
        {
          if (sr->isSetStoichiometryMath())
            {
              ignoreMassAction = true;
            }
          else
            {
              stoi = sr->getStoichiometry() / sr->getDenominator();
            }
        }
      else
        {
          // We may not use Mass Action kinetics if we do not know the actual stoichiometry.
          if (sr->isSetId() &&
              (pSBMLModel->getInitialAssignment(sr->getId()) ||
               pSBMLModel->getRule(sr->getId())))
            {
              ignoreMassAction = true;
            }

          if (sr->isSetStoichiometry())
            {
              stoi = sr->getStoichiometry();
            }
        }

      std::map<std::string, CMetab*>::iterator pos;
      pos = this->speciesMap.find(sr->getSpecies());

      if (pos == this->speciesMap.end())
        {
          delete copasiReaction;
          fatalError();
        }
      else
        {
          // check if there is a conversion factor on the species
          if (this->mpModelConversionFactor != NULL || this->mSpeciesConversionParameterMap.find(pos->second->getSBMLId()) != this->mSpeciesConversionParameterMap.end())
            {
              mConversionFactorNeeded = true;
            }
        }

      std::map<CCopasiObject*, SBase*>::const_iterator spos = copasi2sbmlmap.find(pos->second);
      assert(spos != copasi2sbmlmap.end());
      Species* pSBMLSpecies = dynamic_cast<Species*>(spos->second);
      assert(pSBMLSpecies != NULL);
      hasOnlySubstanceUnitPresent = (hasOnlySubstanceUnitPresent | (pSBMLSpecies->getHasOnlySubstanceUnits() == true));

      if (compartment == NULL)
        {
          compartment = pos->second->getCompartment();
        }
      else
        {
          if (singleCompartment && compartment != pos->second->getCompartment())
            {
              singleCompartment = false;
            }
        }

      copasiReaction->addProduct(pos->second->getKey(), stoi);

      // we need to store the id of the species reference if it is set because SBML Level 3 allows
      // references to species references and if we want to support his, we need the id to import
      // expressions that reference a species reference
      if (this->mLevel > 2)
        {
          CCopasiVector<CChemEqElement>::const_iterator it = copasiReaction->getChemEq().getProducts().begin(), endit = copasiReaction->getChemEq().getProducts().end();
          CChemEqElement* pElement = NULL;

          while (it != endit)
            {
              if ((*it)->getMetaboliteKey() == pos->second->getKey())
                {
                  pElement = const_cast<CChemEqElement*>(*it);
                  break;
                }

              ++it;
            }

          assert(pElement != NULL);

          if (pElement != NULL)
            {
              copasi2sbmlmap[pElement] = const_cast<SpeciesReference*>(sr);
              this->mChemEqElementSpeciesIdMap[pElement] = std::pair<std::string, CChemEq::MetaboliteRole>(sr->getSpecies(), CChemEq::PRODUCT);
            }
        }

      if (sr->isSetStoichiometryMath())
        {
          CChemEq& chemEq = copasiReaction->getChemEq();
          CCopasiVector < CChemEqElement >::const_iterator it = chemEq.getProducts().begin();
          CCopasiVector < CChemEqElement >::const_iterator end = chemEq.getProducts().end();
          CChemEqElement* pChemEqElement = NULL;

          while (it != end)
            {
              if ((*it)->getMetabolite() == pos->second)
                {
                  pChemEqElement = (*it);
                  break;
                }

              ++it;
            }

          assert(pChemEqElement != NULL);
          mStoichiometricExpressionMap.insert(std::make_pair(sr->getStoichiometryMath()->getMath(), pChemEqElement));
        }
    }

  /* Add all modifiers to the reaction */
  num = sbmlReaction->getNumModifiers();

  for (counter = 0; counter < num; counter++)
    {
      const ModifierSpeciesReference* sr = sbmlReaction->getModifier(counter);

      if (sr == NULL)
        {
          delete copasiReaction;
          fatalError();
        }

      std::map<std::string, CMetab*>::iterator pos;
      pos = this->speciesMap.find(sr->getSpecies());

      if (pos == this->speciesMap.end())
        {
          delete copasiReaction;
          fatalError();
        }

      std::map<CCopasiObject*, SBase*>::const_iterator spos = copasi2sbmlmap.find(pos->second);
      assert(spos != copasi2sbmlmap.end());
      Species* pSBMLSpecies = dynamic_cast<Species*>(spos->second);
      assert(pSBMLSpecies != NULL);
      hasOnlySubstanceUnitPresent = (hasOnlySubstanceUnitPresent | (pSBMLSpecies->getHasOnlySubstanceUnits() == true));
      copasiReaction->addModifier(pos->second->getKey());

      // we need to store the id of the species reference if it is set because SBML Level 3 allows
      // references to species references and if we want to support his, we need the id to import
      // expressions that reference a species reference
      if (this->mLevel > 2)
        {

          CCopasiVector<CChemEqElement>::const_iterator it = copasiReaction->getChemEq().getModifiers().begin(), endit = copasiReaction->getChemEq().getModifiers().end();
          CChemEqElement* pElement = NULL;

          while (it != endit)
            {
              if ((*it)->getMetaboliteKey() == pos->second->getKey())
                {
                  pElement = const_cast<CChemEqElement*>(*it);
                  break;
                }

              ++it;
            }

          assert(pElement != NULL);

          if (pElement != NULL)
            {
              copasi2sbmlmap[pElement] = const_cast<ModifierSpeciesReference*>(sr);
            }
        }
    }

  /* in the newly created CFunction set the types for all parameters and
   * either a mapping or a value
   */
  const KineticLaw* kLaw = sbmlReaction->getKineticLaw();

  if (kLaw != NULL)
    {
      const ListOfParameters* pParamList = NULL;

      if (this->mLevel > 2)
        {
          pParamList = kLaw->getListOfLocalParameters();
        }
      else
        {
          pParamList = kLaw->getListOfParameters();
        }

      for (counter = 0; counter < pParamList->size(); ++counter)
        {
          const Parameter* pSBMLParameter = pParamList->get(counter);
          std::string id;

          if (this->mLevel == 1)
            {
              id = pSBMLParameter->getName();
            }
          else
            {
              id = pSBMLParameter->getId();
            }

          double value;

          if (pSBMLParameter->isSetValue() && pSBMLParameter->getValue() == pSBMLParameter->getValue()) // make sure it is not set to NaN
            {
              value = pSBMLParameter->getValue();
            }
          else
            {
              // Set value to NaN and create a warning if it is the first time
              // this happend
              value = std::numeric_limits<C_FLOAT64>::quiet_NaN();

              if (!this->mIncompleteModel)
                {
                  this->mIncompleteModel = true;
                  CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 42, pSBMLParameter->getId().c_str());
                }
            }

          copasiReaction->getParameters().addParameter(id, CCopasiParameter::DOUBLE, value);
        }

      const ASTNode* kLawMath = kLaw->getMath();

      if (kLawMath == NULL || kLawMath->getType() == AST_UNKNOWN)
        {
          copasiReaction->setFunction(NULL);
          CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 56, sbmlReaction->getId().c_str());
        }
      else
        {
          // check for references to species references in the expression because we don't support them yet
          if (!SBMLImporter::findIdInASTTree(kLawMath, this->mSBMLSpeciesReferenceIds).empty())
            {
              CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 95);
            }

          ConverterASTNode* node = new ConverterASTNode(*kLawMath);
          this->preprocessNode(node, pSBMLModel, copasi2sbmlmap, sbmlReaction);

          if (node == NULL)
            {
              delete copasiReaction;
              fatalError();
            }

          /* if it is a single compartment reaction, we have to divide the whole kinetic
          ** equation by the compartment because COPASI expects
          ** kinetic laws that specify concentration/time for single compartment
          ** reactions.
          */
          if (singleCompartment)
            {

              if (compartment != NULL)
                {
                  // only divide if it is not a 0-dimensional compartment
                  if (compartment->getDimensionality() != 0)
                    {
                      // check if the root node is a multiplication node and it's first child
                      // is a compartment node, if so, drop those two and make the second child
                      // new new root node
                      ConverterASTNode* tmpNode1 = this->isMultipliedByVolume(node, compartment->getSBMLId());

                      if (tmpNode1)
                        {
                          delete node;
                          node = tmpNode1;

                          if (node->getType() == AST_DIVIDE && node->getNumChildren() != 2)
                            {
                              delete tmpNode1;
                              fatalError();
                            }
                        }
                      else
                        {
                          tmpNode1 = new ConverterASTNode();
                          tmpNode1->setType(AST_DIVIDE);
                          tmpNode1->addChild(node);
                          ConverterASTNode* tmpNode2 = new ConverterASTNode();
                          tmpNode2->setType(AST_NAME);
                          tmpNode2->setName(compartment->getSBMLId().c_str());
                          tmpNode1->addChild(tmpNode2);
                          node = tmpNode1;
                          std::map<CCopasiObject*, SBase*>::const_iterator pos = copasi2sbmlmap.find(const_cast<CCompartment*>(compartment));
                          assert(pos != copasi2sbmlmap.end());
                          Compartment* pSBMLCompartment = dynamic_cast<Compartment*>(pos->second);
                          assert(pSBMLCompartment != NULL);

                          if (!hasOnlySubstanceUnitPresent && ((this->mLevel == 1 && pSBMLCompartment->isSetVolume()) || (this->mLevel >= 2 && pSBMLCompartment->isSetSize())) && pSBMLCompartment->getSize() == 1.0)
                            {
                              // we have to check if all species used in the reaction
                              // have the hasOnlySubstance flag set

                              if (node->getChild(0)->getType() == AST_FUNCTION && (!this->containsVolume(node->getChild(0), compartment->getSBMLId())))
                                {
                                  // add the id of the reaction to the set so that we can create an error message later.
                                  this->mDivisionByCompartmentReactions.insert(sbmlReaction->getId());
                                }
                            }
                        }
                    }
                }
              else
                {
                  delete node;
                  delete copasiReaction;
                  fatalError();
                }
            }

          /* Create a new user defined CKinFunction */
          if (!sbmlId2CopasiCN(node, copasi2sbmlmap, copasiReaction->getParameters()))
            {
              CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 27, copasiReaction->getObjectName().c_str());
            }

          CEvaluationNode* pExpressionTreeRoot = CEvaluationTree::fromAST(node);
          delete node;
          node = NULL;

          if (pExpressionTreeRoot)
            {
              CExpression KineticLawExpression;
              KineticLawExpression.setRoot(pExpressionTreeRoot);

              // check if the expression is constant flux
              CCopasiObject* pParamObject = SBMLImporter::isConstantFlux(pExpressionTreeRoot, copasiModel, pTmpFunctionDB);

              if (pParamObject != NULL)
                {
                  // assert that the object really is a local or global
                  // parameter
                  assert(dynamic_cast<const CCopasiParameter*>(pParamObject) || dynamic_cast<const CModelValue*>(pParamObject));
                  std::string functionName;

                  if (copasiReaction->isReversible())
                    {
                      // set the function to Constant flux (reversible)
                      functionName = "Constant flux (reversible)";
                    }
                  else
                    {
                      // set the function to Constant flux (irreversible)
                      functionName = "Constant flux (irreversible)";
                    }

                  CFunction* pCFFun = dynamic_cast<CFunction*>(this->functionDB->findFunction(functionName));
                  assert(pCFFun != NULL);
                  CEvaluationNodeCall* pCallNode = NULL;

                  if (CEvaluationNode::type(pExpressionTreeRoot->getType()) == CEvaluationNode::OBJECT)
                    {
                      pCallNode = new CEvaluationNodeCall(CEvaluationNodeCall::EXPRESSION, "dummy_call");
                      // add the parameter
                      pCallNode->addChild(pExpressionTreeRoot->copyBranch());
                    }
                  else
                    {
                      pCallNode = dynamic_cast<CEvaluationNodeCall*>(pExpressionTreeRoot->copyBranch());
                      assert(pCallNode != NULL);
                    }

                  if (pParamObject->getObjectType() == "Parameter")
                    {
                      pParamObject->setObjectName("v");
                      dynamic_cast<CEvaluationNode*>(pCallNode->getChild())->setData("<" + pParamObject->getCN() + ">");
                    }

                  copasiReaction->setFunction(pCFFun);
                  // map the parameter
                  this->doMapping(copasiReaction, pCallNode);
                  delete pCallNode;
                }
              else
                {
                  // check if the root node is a simple function call
                  if (this->isSimpleFunctionCall(pExpressionTreeRoot))
                    {
                      // if yes, we check if it corresponds to an already existing function
                      std::string functionName = pExpressionTreeRoot->getData();
                      CFunction* pImportedFunction = dynamic_cast<CFunction*>(functionDB->findFunction(functionName));
                      assert(pImportedFunction);
                      std::vector<CEvaluationNodeObject*>* v = NULL;

                      // only check for mass action if there is no conversion factor involved
                      // for any of the species involved in the reaction (substrates and products)
                      if (!mConversionFactorNeeded && !ignoreMassAction)
                        {
                          v = this->isMassAction(pImportedFunction, copasiReaction->getChemEq(), static_cast<const CEvaluationNodeCall*>(pExpressionTreeRoot));

                          if (!v)
                            {
                              CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 27, copasiReaction->getObjectName().c_str());
                            }
                        }

                      if (v && !v->empty())
                        {
                          CFunction* pFun = NULL;

                          if (copasiReaction->isReversible())
                            {
                              pFun = static_cast<CFunction*>(this->functionDB->findFunction("Mass action (reversible)"));
                            }
                          else
                            {
                              pFun = static_cast<CFunction*>(this->functionDB->findFunction("Mass action (irreversible)"));
                            }

                          if (!pFun)
                            {
                              fatalError();
                            }

                          // do the mapping
                          CEvaluationNodeCall* pCallNode = new CEvaluationNodeCall(CEvaluationNodeCall::EXPRESSION, "dummy_call");
                          size_t i, iMax = v->size();

                          for (i = 0; i < iMax; ++i)
                            {
                              pCallNode->addChild((*v)[i]);
                            }

                          this->renameMassActionParameters(pCallNode);
                          copasiReaction->setFunction(pFun);
                          this->doMapping(copasiReaction, pCallNode);
                          delete pCallNode;
                        }
                      else
                        {
                          // find corresponding reaction does *not* return the function if it is already in the function db
                          // we better change this to return this instance.
                          CFunction* pExistingFunction = findCorrespondingFunction(&KineticLawExpression, copasiReaction);

                          // if it does, we set the existing function for this reaction
                          if (pExistingFunction)
                            {
                              copasiReaction->setFunction(pExistingFunction);
                              // do the mapping
                              doMapping(copasiReaction, dynamic_cast<const CEvaluationNodeCall*>(pExpressionTreeRoot));
                            }
                          // else we take the function from the pTmpFunctionDB, copy it and set the usage correctly
                          else
                            {
                              // replace the variable nodes in pImportedFunction with  nodes from
                              std::map<std::string, std::string > arguments;
                              const CFunctionParameters& funParams = pImportedFunction->getVariables();
                              const CEvaluationNode* pTmpNode = static_cast<const CEvaluationNode*>(pExpressionTreeRoot->getChild());
                              size_t i, iMax = funParams.size();

                              for (i = 0; (i < iMax) && pTmpNode; ++i)
                                {
                                  if (!(pTmpNode->getType() == CEvaluationNode::OBJECT)) fatalError();

                                  arguments[funParams[i]->getObjectName()] = pTmpNode->getData().substr(1, pTmpNode->getData().length() - 2);
                                  pTmpNode = static_cast<const CEvaluationNode*>(pTmpNode->getSibling());
                                }

                              assert((i == iMax) && pTmpNode == NULL);

                              CEvaluationNode* pTmpExpression = variables2objects(pImportedFunction->getRoot()->copyBranch(), arguments);

                              CExpression TmpTree2;
                              TmpTree2.setRoot(pTmpExpression);

                              // code to fix bug 1874
                              // since this is a user defined function, we want to retain the original name
                              // next setFunctionFromExpressionTree will take this name if it is not 'Expression'
                              // (the default)
                              TmpTree2.setObjectName(functionName);

                              CFunction * pNewFunction = copasiReaction->setFunctionFromExpressionTree(TmpTree2, copasi2sbmlmap, this->functionDB);

                              if (pNewFunction != NULL &&
                                  pNewFunction->getType() == CEvaluationTree::UserDefined)
                                {
                                  // code to fix Bug 1015
                                  if (!pNewFunction->isSuitable(copasiReaction->getChemEq().getSubstrates().size(),
                                                                copasiReaction->getChemEq().getProducts().size(),
                                                                copasiReaction->isReversible() ? TriTrue : TriFalse))
                                    {
                                      pNewFunction->setReversible(TriUnspecified);
                                    }

                                  pTmpFunctionDB->add(pNewFunction, false);
                                  mUsedFunctions.insert(pNewFunction->getObjectName());
                                }
                            }
                        }

                      pdelete(v);
                    }
                  else
                    {
                      std::vector<CEvaluationNodeObject*>* v = NULL;

                      // only check for mass action if there is no conversion factor involved
                      // for any of the species involved in the reaction (substrates and products)
                      if (!mConversionFactorNeeded && !ignoreMassAction)
                        {
                          v = this->isMassAction(&KineticLawExpression, copasiReaction->getChemEq());

                          if (!v)
                            {
                              CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 27, copasiReaction->getObjectName().c_str());
                            }
                        }

                      if (v && !v->empty())
                        {
                          CFunction* pFun = NULL;

                          if (copasiReaction->isReversible())
                            {
                              pFun = static_cast<CFunction*>(this->functionDB->findFunction("Mass action (reversible)"));
                            }
                          else
                            {
                              pFun = static_cast<CFunction*>(this->functionDB->findFunction("Mass action (irreversible)"));
                            }

                          if (!pFun)
                            {
                              fatalError();
                            }

                          // do the mapping
                          CEvaluationNodeCall* pCallNode = new CEvaluationNodeCall(CEvaluationNodeCall::EXPRESSION, "dummy_call");
                          size_t i, iMax = v->size();

                          for (i = 0; i < iMax; ++i)
                            {
                              pCallNode->addChild((*v)[i]);
                            }

                          // rename the function parameters to k1 and k2
                          this->renameMassActionParameters(pCallNode);
                          copasiReaction->setFunction(pFun);
                          this->doMapping(copasiReaction, pCallNode);
                          delete pCallNode;
                        }
                      else
                        {
                          CFunction* pNonconstFun = copasiReaction->setFunctionFromExpressionTree(KineticLawExpression, copasi2sbmlmap, this->functionDB);

                          if (pNonconstFun != NULL &&
                              pNonconstFun->getType() == CEvaluationTree::UserDefined)
                            {
                              // code to fix Bug 1015
                              if (!pNonconstFun->isSuitable(copasiReaction->getChemEq().getSubstrates().size(),
                                                            copasiReaction->getChemEq().getProducts().size(),
                                                            copasiReaction->isReversible() ? TriTrue : TriFalse))
                                {
                                  pNonconstFun->setReversible(TriUnspecified);
                                }

                              pTmpFunctionDB->add(pNonconstFun, false);
                              mUsedFunctions.insert(pNonconstFun->getObjectName());
                            }
                        }

                      pdelete(v);
                    }
                }
            }
          else
            {
              // error message
              CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 8, copasiReaction->getObjectName().c_str());
            }
        }
    }
  else
    {
      /* if no KineticLaw was defined for the reaction. */
      copasiReaction->setFunction(NULL);
    }

  //DebugFile << "Created reaction: " << copasiReaction->getObjectName() << std::endl;
  SBMLImporter::importMIRIAM(sbmlReaction, copasiReaction);
  SBMLImporter::importNotes(copasiReaction, sbmlReaction);

  return copasiReaction;
}

/**
 * Creates a map of each parameter of the function definition and its
 * corresponding parameter in the function call.
 */
std::map<std::string, ASTNode*>
SBMLImporter::createBVarMap(const ASTNode* uDefFunction, const ASTNode* function)
{
  /* the first n-1 children, where n is the number of children, of a function definition ASTnode are the
   * arguments to the function. These correspond to the m=n-1 children of the
   * function call.
   */
  if (uDefFunction->getNumChildren() != function->getNumChildren() + 1)
    {
      std::string functionName = uDefFunction->getName();
      fatalError();
    }

  std::map<std::string, ASTNode*> varMap;
  unsigned int counter;

  for (counter = 0; counter < uDefFunction->getNumChildren() - 1; counter++)
    {
      varMap[uDefFunction->getChild(counter)->getName()] = function->getChild(counter);
    }

  return varMap;
}

/**
 * Returns the user defined SBML function definition that belongs to the given
 * name, or NULL if none can be found.
 */
const FunctionDefinition*
SBMLImporter::getFunctionDefinitionForName(const std::string name, const Model* sbmlModel)
{
  const FunctionDefinition* fDef = NULL;
  unsigned int counter;

  for (counter = 0; counter < sbmlModel->getNumFunctionDefinitions(); counter++)
    {
      std::string functionName = sbmlModel->getFunctionDefinition(counter)->getName();

      if (sbmlModel->getFunctionDefinition(counter)->isSetId())
        {
          functionName = sbmlModel->getFunctionDefinition(counter)->getId();
        }

      if (functionName == name)
        {
          fDef = sbmlModel->getFunctionDefinition(counter);
          break;
        }
    }

  return fDef;
}

#ifdef XXXX
/**
 * Replaces the variables in a function definition with the actual function
 * parameters that were used when the function was called. The function returns
 * a pointer to the ConverterAST node with the replaced variables.
 */
ConverterASTNode*
SBMLImporter::replaceBvars(const ASTNode* node, std::map<std::string, ASTNode*> bvarMap)
{
  ConverterASTNode* newNode = NULL;

  if (node->isName())
    {
      /* check if name matches any in bvarMap */
      /* if yes, replace node with node in bvarMap */
      /* node needs to be set to be a deep copy of the replacement */
      if (bvarMap.find(node->getName()) != bvarMap.end())
        {
          newNode = new ConverterASTNode(*bvarMap[node->getName()]);
        }
    }
  else
    {
      newNode = new ConverterASTNode(*node);
      newNode->setChildren(new List());
      unsigned int counter;

      for (counter = 0; counter < node->getNumChildren(); counter++)
        {
          newNode->addChild(this->replaceBvars(node->getChild(counter), bvarMap));
        }
    }

  return newNode;
}
#endif // XXXX

/**
 * Constructor that initializes speciesMap and the FunctionDB object
 */
SBMLImporter::SBMLImporter():
  mIgnoredSBMLMessages(),
  speciesMap(),
  functionDB(NULL),
  mIncompleteModel(false),
  mUnsupportedRuleFound(false),
  mUnsupportedRateRuleFound(false),
  mUnsupportedAssignmentRuleFound(false),
  mUnitOnNumberFound(false),
  mAssignmentToSpeciesReferenceFound(false),
  mLevel(0),
  mOriginalLevel(0),
  mVersion(0),
  sbmlIdMap(),
  mUsedFunctions(),
  mpDataModel(NULL),
  mpCopasiModel(NULL),
  mFunctionNameMapping(),
  mDivisionByCompartmentReactions(),
  mpImportHandler(NULL),
  mImportStep(0),
  mhImportStep(C_INVALID_INDEX),
  mTotalSteps(0),
  mSubstanceOnlySpecies(),
  mFastReactions(),
  mReactionsWithReplacedLocalParameters(),
  mExplicitelyTimeDependentFunctionDefinitions(),
  mIgnoredParameterUnits(),
  mStoichiometricExpressionMap(),
  mDelayFound(false),
  mPotentialAvogadroNumbers(),
  mAvogadroCreated(false),
  mImportCOPASIMIRIAM(true),
  mDelayNodeMap(),
  mUsedSBMLIds(),
  mUsedSBMLIdsPopulated(false),
  mAvogadroSet(false),
  mKnownCustomUserDefinedFunctions(),
#if LIBSBML_VERSION >= 40100
  mpModelConversionFactor(NULL),
  mChemEqElementSpeciesIdMap(),
  mSpeciesConversionParameterMap(),
  mSBMLIdModelValueMap(),
  mSBMLSpeciesReferenceIds(),
  mRateRuleForSpeciesReferenceIgnored(false),
  mEventAssignmentForSpeciesReferenceIgnored(false),
  mConversionFactorFound(false),
# if LIBSBML_VERSION >= 40200
  mEventPrioritiesIgnored(false),
  mInitialTriggerValues(false),
  mNonPersistentTriggerFound(false)
# endif // LIBSBML_VERSION >= 40200
#endif // LIBSBML_VERSION >= 40100
{
  this->speciesMap = std::map<std::string, CMetab*>();
  this->functionDB = NULL;
  this->mIncompleteModel = false;
  this->mUnsupportedRuleFound = false;
  this->mUnsupportedRateRuleFound = false;
  this->mUnsupportedAssignmentRuleFound = false;
  this->mUnitOnNumberFound = false;
  this->mAssignmentToSpeciesReferenceFound = false;
  this->mpImportHandler = NULL;
  this->mDelayFound = false;
  this->mAvogadroCreated = false;
  // these data structures are used to handle the new conversion factores in
  // SBML L3 models and references to species references
  this->mpModelConversionFactor = NULL;
  this->mChemEqElementSpeciesIdMap.clear();
  this->mSpeciesConversionParameterMap.clear();
  this->mSBMLIdModelValueMap.clear();
  this->mRateRuleForSpeciesReferenceIgnored = false;
  this->mEventAssignmentForSpeciesReferenceIgnored = false;
  this->mConversionFactorFound = false;
  this->mEventPrioritiesIgnored = false;
  this->mInitialTriggerValues = false;
  this->mNonPersistentTriggerFound = false;

  this->mIgnoredSBMLMessages.insert(10501);
  this->mIgnoredSBMLMessages.insert(10512);
  this->mIgnoredSBMLMessages.insert(10513);
  this->mIgnoredSBMLMessages.insert(10522);
  this->mIgnoredSBMLMessages.insert(10533);
  this->mIgnoredSBMLMessages.insert(10541);
  this->mIgnoredSBMLMessages.insert(10551);
  this->mIgnoredSBMLMessages.insert(10562);
  this->mIgnoredSBMLMessages.insert(80701);
  this->mIgnoredSBMLMessages.insert(99505);
}

/**
 * Destructor that does nothing.
 */
SBMLImporter::~SBMLImporter()
{}

/**
 * This functions replaces all species nodes for species that are in the substanceOnlySpeciesVector.
 * With the node multiplied by the volume of the species compartment.
void SBMLImporter::replaceSubstanceOnlySpeciesNodes(ConverterASTNode* node, const std::map<Species*, Compartment*>& substanceOnlySpecies)
{
  if (node != NULL)
    {
      if (node->getType() == AST_NAME)
        {
          std::map<Species*, Compartment*>::const_iterator it = substanceOnlySpecies.begin();
          std::map<Species*, Compartment*>::const_iterator endIt = substanceOnlySpecies.end();
          while (it != endIt)
            {
              if (it->first->getId() == node->getName())
                {
                  // replace node
                  List* l = new List();
                  ConverterASTNode* child1 = new ConverterASTNode(AST_NAME);
                  child1->setName(node->getName());
                  ConverterASTNode* child2 = new ConverterASTNode(AST_NAME);
                  child2->setName(it->second->getId().c_str());
                  l->add(child1);
                  l->add(child2);
                  node->setChildren(l);
                  node->setType(AST_TIMES);
                  break;
                }
              ++it;
            }
        }
      else
        {
          unsigned int counter;
          for (counter = 0;counter < node->getNumChildren();counter++)
            {
              this->replaceSubstanceOnlySpeciesNodes((ConverterASTNode*)node->getChild(counter), substanceOnlySpecies);
            }
        }
    }
}
 */

/**
 * Function reads an SBML file with libsbml and converts it to a Copasi CModel
 */
CModel* SBMLImporter::readSBML(std::string filename,
                               CFunctionDB* funDB,
                               SBMLDocument*& pSBMLDocument,
                               std::map<CCopasiObject*, SBase*>& copasi2sbmlmap,
                               CListOfLayouts *& prLol,
                               CCopasiDataModel* pDataModel)
{
  // convert filename to the locale encoding
  std::ifstream file(CLocaleString::fromUtf8(filename).c_str());

  if (!file)
    {
      CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 50, filename.c_str());
    }

  std::ostringstream stringStream;
  char c;

  while (file.get(c))
    {
      stringStream << c;
    }

  file.clear();
  file.close();
  return this->parseSBML(stringStream.str(), funDB,
                         pSBMLDocument, copasi2sbmlmap, prLol, pDataModel);
}

/**
 * Function parses an SBML document with libsbml and converts it to a COPASI CModel
 * object which is returned. Deletion of the returned pointer is up to the
 * caller.
 */
CModel*
SBMLImporter::parseSBML(const std::string& sbmlDocumentText,
                        CFunctionDB* funDB,
                        SBMLDocument *& pSBMLDocument,
                        std::map<CCopasiObject*, SBase*>& copasi2sbmlmap,
                        CListOfLayouts *& prLol,
                        CCopasiDataModel* pDataModel)
{
  this->mUsedSBMLIdsPopulated = false;
  this->mAvogadroSet = false;
  mpDataModel = pDataModel;
  assert(mpDataModel != NULL);

  this->mpCopasiModel = NULL;

  if (funDB != NULL)
    {
      this->functionDB = funDB;
      SBMLReader* reader = new SBMLReader();

      mImportStep = 0;

      if (mpImportHandler)
        {
          mpImportHandler->setName("Importing SBML file...");
          mTotalSteps = 11;
          mhImportStep = mpImportHandler->addItem("Step",
                                                  mImportStep,
                                                  &mTotalSteps);
        }

      unsigned C_INT32 step = 0, totalSteps = 0;
      size_t hStep = C_INVALID_INDEX;

      if (this->mpImportHandler != 0)
        {
          step = 0;
          totalSteps = 1;
          hStep = mpImportHandler->addItem("Reading SBML file...",
                                           step,
                                           &totalSteps);
        }

      SBMLDocument* sbmlDoc = reader->readSBMLFromString(sbmlDocumentText);

      if (mpImportHandler) mpImportHandler->finishItem(hStep);

      if (this->mpImportHandler != 0)
        {
          step = 0;
          totalSteps = 1;
          hStep = mpImportHandler->addItem("Checking consistency...",
                                           step,
                                           &totalSteps);
        }

      if (CCopasiRootContainer::getConfiguration()->validateUnits())
        sbmlDoc->setApplicableValidators(AllChecksON);
      else
        sbmlDoc->setApplicableValidators(AllChecksON & UnitsCheckOFF);

#if LIBSBML_VERSION > 50800
      // libSBML is validating comp models after 5.8.0 this would throw an
      // error in case external references can't be resolved.
      sbmlDoc->setLocationURI(pDataModel->getReferenceDirectory());
#endif

      bool checkResult = sbmlDoc->checkConsistency();

#if LIBSBML_VERSION > 50800

      sbmlDoc->setLocationURI(pDataModel->getReferenceDirectory());

      // the new libsbml includes complete layout validation, and flags many
      // things as errors, that would cause COPASI to reject the model (even
      // though the layout code has been written to anticipate all these possible
      // error sources). Thus downgrade these error messages
      sbmlDoc->getErrorLog()->changeErrorSeverity(LIBSBML_SEV_ERROR, LIBSBML_SEV_WARNING, "layout");
#endif

      if (mpImportHandler) mpImportHandler->finishItem(hStep);

      if (checkResult != 0)
        {
          int fatal = -1;
          unsigned int i, iMax = sbmlDoc->getNumErrors();

          for (i = 0; (i < iMax) && (fatal == -1); ++i)
            {
              const XMLError* pSBMLError = sbmlDoc->getError(i);

              // we check if the model contained a required package
              if (sbmlDoc->getLevel() > 2 && pSBMLError->getErrorId() == 99107)
                {
                  CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 96);
                }

              CCopasiMessage::Type messageType = CCopasiMessage::RAW;

              switch (pSBMLError->getSeverity())
                {
                  case LIBSBML_SEV_INFO:

                    if (mIgnoredSBMLMessages.find(pSBMLError->getErrorId()) != mIgnoredSBMLMessages.end())
                      {
                        messageType = CCopasiMessage::WARNING_FILTERED;
                      }
                    else
                      {
                        messageType = CCopasiMessage::WARNING;
                      }

                    CCopasiMessage(messageType, MCSBML + 40, "INFO", pSBMLError->getErrorId(), pSBMLError->getLine(), pSBMLError->getColumn(), pSBMLError->getMessage().c_str());
                    break;

                  case LIBSBML_SEV_WARNING:

#if LIBSBML_VERSION > 50800

                    // filter layout warnings always
                    if (pSBMLError->getPackage() == "layout")
                      messageType = CCopasiMessage::WARNING_FILTERED;
                    else
#endif
                      if (mIgnoredSBMLMessages.find(pSBMLError->getErrorId()) != mIgnoredSBMLMessages.end())
                        {
                          messageType = CCopasiMessage::WARNING_FILTERED;
                        }
                      else
                        {
                          messageType = CCopasiMessage::WARNING;
                        }

                    CCopasiMessage(messageType, MCSBML + 40, "WARNING", pSBMLError->getErrorId(), pSBMLError->getLine(), pSBMLError->getColumn(), pSBMLError->getMessage().c_str());
                    break;

                  case LIBSBML_SEV_ERROR:

                    if (mIgnoredSBMLMessages.find(pSBMLError->getErrorId()) != mIgnoredSBMLMessages.end())
                      {
                        messageType = CCopasiMessage::ERROR_FILTERED;
                      }

                    CCopasiMessage(messageType, MCSBML + 40, "ERROR", pSBMLError->getErrorId(), pSBMLError->getLine(), pSBMLError->getColumn(), pSBMLError->getMessage().c_str());
                    break;

                  case LIBSBML_SEV_FATAL:

                    // treat unknown as fatal
                  default:

                    //CCopasiMessage(CCopasiMessage::TRACE, MCSBML + 40,"FATAL",pSBMLError->getLine(),pSBMLError->getColumn(),pSBMLError->getMessage().c_str());
                    if (pSBMLError->getErrorId() == 10804)
                      {
                        // this error indicates a problem with a notes element
                        // although libsbml flags this as fatal, we would still
                        // like to read the model
                        CCopasiMessage(messageType, MCSBML + 40, "ERROR", pSBMLError->getErrorId(), pSBMLError->getLine(), pSBMLError->getColumn(), pSBMLError->getMessage().c_str());
                      }
                    else
                      {
                        fatal = i;
                      }

                    break;
                }

              //std::cerr << pSBMLError->getMessage() << std::endl;
            }

          if (fatal != -1)
            {
              const XMLError* pSBMLError = sbmlDoc->getError(fatal);
              CCopasiMessage Message(CCopasiMessage::EXCEPTION, MCXML + 2,
                                     pSBMLError->getLine(),
                                     pSBMLError->getColumn(),
                                     pSBMLError->getMessage().c_str());

              if (mpImportHandler) mpImportHandler->finishItem(mhImportStep);

              return NULL;
            }
        }
      else
        {
          if (sbmlDoc->getLevel() > 2)
            {
              // we check if the model contained a required package
              unsigned int i, iMax = sbmlDoc->getNumErrors();

              for (i = 0; i < iMax ; ++i)
                {
                  const XMLError* pSBMLError = sbmlDoc->getError(i);

                  if (pSBMLError->getErrorId() == 99107)
                    {
                      if (pSBMLError->getMessage().find("'spatial'") != std::string::npos)
                        {
                          CCopasiMessage(CCopasiMessage::ERROR,
                                         "The model you tried to open requires the SBML spatial package. "
                                         "This version of COPASI does not support spatial models.");
                        }
                      else
                        {
                          CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 96);
                        }
                    }
                }
            }
        }

#if LIBSBML_VERSION >= 50700

      if (sbmlDoc->getPlugin("comp") != NULL && sbmlDoc->isSetPackageRequired("comp"))
        {

          if (sbmlDoc->getNumErrors(LIBSBML_SEV_ERROR) > 0)
            {
              std::string message =
                "The SBML model you are trying to import uses the Hierarchical Modeling extension. "
                "In order to import this model in COPASI, it has to be flattened, however flattening is "
                "not possible because of the following errors:\n" + sbmlDoc->getErrorLog()->toString();

              // escape potential percent signs that would lead to an error when passed to CCopasiMessage
              CPrefixNameTransformer::replaceStringInPlace(message, "%", "%%");

              CCopasiMessage(CCopasiMessage::EXCEPTION, message.c_str());
            }

#if LIBSBML_VERSION >= 50903 && LIBSBML_HAS_PACKAGE_COMP
          // apply the name transformer
          CompModelPlugin* mPlug = dynamic_cast<CompModelPlugin*>(sbmlDoc->getModel()->getPlugin("comp"));
          CPrefixNameTransformer trans;

          if (mPlug != NULL)
            {
              mPlug->setTransformer(&trans);
            }

#endif //LIBSBML_VERSION >= 50903 && LIBSBML_HAS_PACKAGE_COMP

          // the sbml comp package is used, and the required flag is set, so it stands to reason
          // that we need to flatten the document
          sbmlDoc->getErrorLog()->clearLog();

          ConversionProperties props;
          props.addOption("flatten comp");
          props.addOption("leavePorts", false);
          props.addOption("basePath", pDataModel->getReferenceDirectory());

          if (sbmlDoc->convert(props) != LIBSBML_OPERATION_SUCCESS)
            {
              std::string message =
                "The SBML model you are trying to import uses the Hierarchical Modeling extension. "
                "In order to import this model in COPASI, it has to be flattened, however flattening failed";

              if (sbmlDoc->getNumErrors() == 0)
                message += ".";
              else
                message += " with the following errors:\n" + sbmlDoc->getErrorLog()->toString();

              CCopasiMessage(CCopasiMessage::EXCEPTION, message.c_str());
            }
        }

#endif

      if (sbmlDoc->getModel() == NULL)
        {
          CCopasiMessage Message(CCopasiMessage::ERROR, MCSBML + 2);

          if (mpImportHandler) mpImportHandler->finishItem(mhImportStep);

          return NULL;
        }

      delete reader;
      pSBMLDocument = sbmlDoc;
      this->mLevel = pSBMLDocument->getLevel();
      // remember the original level of the document because we convert Level 1 documents to Level 2
      // For the import of rules, we need to remember that is was actually a level 1 document
      // because otherwise we throw error messages on rules on parameters since the parameters
      //  have been set to constant by the conversation to Level 2
      this->mOriginalLevel = this->mLevel;
      this->mVersion = pSBMLDocument->getVersion();

      if (mLevel == 1)
        {
          unsigned int i, iMax = pSBMLDocument->getModel()->getNumCompartments();

          for (i = 0; i < iMax; ++i)
            {
              Compartment* pCompartment = pSBMLDocument->getModel()->getCompartment(i);
              pCompartment->setSize(pCompartment->getVolume());
            }

          pSBMLDocument->setLevelAndVersion(2, 1);
          mLevel = pSBMLDocument->getLevel();
        }

      this->mpCopasiModel = this->createCModelFromSBMLDocument(sbmlDoc, copasi2sbmlmap);

      prLol = new CListOfLayouts("ListOfLayouts", mpDataModel);
      Model* sbmlmodel = pSBMLDocument->getModel();

      if (sbmlmodel && prLol)
        {
          LayoutModelPlugin *lmPlugin = (LayoutModelPlugin*)sbmlmodel->getPlugin("layout");

          if (lmPlugin != NULL)
            SBMLDocumentLoader::readListOfLayouts(*prLol,
                                                  *lmPlugin->getListOfLayouts(),
                                                  copasi2sbmlmap);
        }
    }
  else
    {
      if (mpImportHandler) mpImportHandler->finishItem(mhImportStep);

      fatalError();
    }

  if (mpImportHandler) mpImportHandler->finishItem(mhImportStep);

  return this->mpCopasiModel;
}

/**
 * Returns the copasi QuantityUnit corresponding to the given SBML
 *  Substance UnitDefinition.
 */
std::pair<CModel::QuantityUnit, bool>
SBMLImporter::handleSubstanceUnit(const UnitDefinition* uDef)
{
  bool result = false;
  CModel::QuantityUnit qUnit = CModel::Mol;

  if (uDef == NULL)
    {
      //DebugFile << "Argument to handleSubstanceUnit is NULL pointer." << std::endl;
      fatalError();
    }

  if (uDef->getNumUnits() == 1)
    {
      const Unit* u = uDef->getUnit(0);

      if (u == NULL)
        {
          //DebugFile << "Expected Unit, got NULL pointer." << std::endl;
          fatalError();
        }

      if ((u->getKind() == UNIT_KIND_MOLE)
          || u->getKind() == UNIT_KIND_AVOGADRO
         )
        {
          double multiplier = u->getMultiplier();
          int scale = u->getScale();

          if (multiplier != 1)
            {
              // check if the multiplier is a multiple of 10
              // so that we might be able to convert it to a scale that makes
              // sense
              double tmp = log10(multiplier);

              if (areApproximatelyEqual(tmp, round(tmp)))
                {
                  scale += (int) round(tmp);
                  multiplier = 1;
                }
            }

          if ((u->getExponent() == 1) &&
              areApproximatelyEqual(multiplier, 1.0) &&
              ((scale % 3) == 0) &&
              (scale < 1) &&
              (scale > -16))
            {
              switch (scale)
                {
                  case 0:
                    qUnit = CModel::Mol;
                    result = true;
                    break;

                  case - 3:
                    qUnit = CModel::mMol;
                    result = true;
                    break;

                  case - 6:
                    qUnit = CModel::microMol;
                    result = true;
                    break;

                  case - 9:
                    qUnit = CModel::nMol;
                    result = true;
                    break;

                  case - 12:
                    qUnit = CModel::pMol;
                    result = true;
                    break;

                  case - 15:
                    qUnit = CModel::fMol;
                    result = true;
                    break;

                  default:
                    //DebugFile << "Error. This value should never have been reached for the scale of the liter unit." << std::endl;
                    result = false;
                    break;
                }
            }
          else
            {
              result = false;
            }
        }
      else if ((u->getKind() == UNIT_KIND_ITEM))
        {
          double multiplier = u->getMultiplier();
          int scale = u->getScale();

          if (multiplier != 1)
            {
              // check if the multiplier is a multiple of 10
              // so that we might be able to convert it to a scale that makes
              // sense
              double tmp = log10(multiplier);

              if (areApproximatelyEqual(tmp, round(tmp)))
                {
                  scale += (int)round(tmp);
                  multiplier = 1;
                }
            }

          if ((u->getExponent() == 1) &&
              areApproximatelyEqual(multiplier, 1.0) &&
              (scale == 0 || scale == 1))
            {
              if (u->getScale() == 1)
                {
                  CCopasiMessage Message(CCopasiMessage::ERROR, MCSBML + 30);
                }
              else
                {
                  result = true;
                  qUnit = CModel::number;
                }
            }
          else
            {
              result = false;
            }
        }
      else if ((u->getKind() == UNIT_KIND_DIMENSIONLESS))
        {
          double multiplier = u->getMultiplier();
          int scale = u->getScale();

          if (multiplier != 1)
            {
              // check if the multiplier is a multiple of 10
              // so that we might be able to convert it to a scale that makes
              // sense
              double tmp = log10(multiplier);

              if (areApproximatelyEqual(tmp, round(tmp)))
                {
                  scale += (int)round(tmp);
                  multiplier = 1;
                }
            }

          if ((u->getExponent() == 1) &&
              areApproximatelyEqual(multiplier, 1.0) &&
              scale == 0)
            {
              result = true;
              qUnit = CModel::dimensionlessQuantity;
            }
          else
            {
              result = false;
            }
        }
      else
        {
          result = false;
        }
    }
  else
    {
      result = false;
    }

  return std::make_pair(qUnit, result);
}

/**
 * Returns the copasi TimeUnit corresponding to the given SBML Time
 *  UnitDefinition.
 */
std::pair<CModel::TimeUnit, bool>
SBMLImporter::handleTimeUnit(const UnitDefinition* uDef)
{
  bool result = false;
  CModel::TimeUnit tUnit = CModel::s;

  if (uDef == NULL)
    {
      //DebugFile << "Argument to handleTimeUnit is NULL pointer." << std::endl;
      fatalError();
    }

  if (uDef->getNumUnits() == 1)
    {
      const Unit* u = uDef->getUnit(0);

      if (u == NULL)
        {
          //DebugFile << "Expected Unit, got NULL pointer." << std::endl;
          fatalError();
        }

      if ((u->getKind() == UNIT_KIND_SECOND))
        {
          double multiplier = u->getMultiplier();
          int scale = u->getScale();

          // this is more difficult, we have to check several possible
          // combinations of scales and multipliers.
          // Valid multipliers are 60, 3600 and 86400 which correspond to a
          // minute an hour and a day each of those has to have a scale of 0
          if (scale == 0)
            {
              // check if the multiplier is 1, 60, 3600 or 86400
              // if not, try to make the multiplier 1
              if (!areApproximatelyEqual(multiplier, 1.0) &&
                  !areApproximatelyEqual(multiplier, 60.0) &&
                  !areApproximatelyEqual(multiplier, 3600.0) &&
                  !areApproximatelyEqual(multiplier, 86400.0))
                {
                  double tmp = log10(multiplier);

                  if (!areApproximatelyEqual(tmp, round(tmp)))
                    {
                      result = false;
                    }
                  else
                    {
                      multiplier = 1;
                      scale += (int)round(tmp);
                    }
                }
            }
          else
            {
              // the multiplier must be 1
              if (!areApproximatelyEqual(multiplier, 1.0))
                {
                  // make the multiplier 1 and check if the scale is an integer,
                  // if not, try to make the 0 and check if the multiplier becomes one
                  // of the valid multipliers 1,60, 3600 or 86400
                  double tmp = log10(multiplier);

                  if (!areApproximatelyEqual(tmp, round(tmp)))
                    {
                      // try to make the scale 0
                      multiplier *= pow(10.0, (double)scale);
                      scale = 0;
                    }
                  else
                    {
                      // make the multiplier 1
                      multiplier = 1;
                      scale += (int) round(tmp);
                    }
                }
            }

          if ((u->getExponent() == 1) && ((scale % 3) == 0) && (scale < 1) && (scale > -16))
            {
              if (areApproximatelyEqual(multiplier, 1.0))
                {
                  switch (scale)
                    {
                      case 0:
                        tUnit = CModel::s;
                        result = true;
                        break;

                      case - 3:
                        tUnit = CModel::ms;
                        result = true;
                        break;

                      case - 6:
                        tUnit = CModel::micros;
                        result = true;
                        break;

                      case - 9:
                        tUnit = CModel::ns;
                        result = true;
                        break;

                      case - 12:
                        tUnit = CModel::ps;
                        result = true;
                        break;

                      case - 15:
                        tUnit = CModel::fs;
                        result = true;
                        break;

                      default:
                        //DebugFile << "Error. This value should never have been reached for the scale of the time unit." << std::endl;
                        result = false;
                        break;
                    }
                }
              else if ((scale == 0) &&
                       areApproximatelyEqual(multiplier, 60.0))
                {
                  tUnit = CModel::min;
                  result = true;
                }
              else if ((scale == 0) &&
                       areApproximatelyEqual(multiplier, 3600.0))
                {
                  tUnit = CModel::h;
                  result = true;
                }
              else if ((scale == 0) &&
                       areApproximatelyEqual(multiplier, 86400.0))
                {
                  tUnit = CModel::d;
                  result = true;
                }
              else
                {
                  result = false;
                }
            }
          else
            {
              result = false;
            }
        }
      else if ((u->getKind() == UNIT_KIND_DIMENSIONLESS))
        {
          double multiplier = u->getMultiplier();
          int scale = u->getScale();

          if (multiplier != 1)
            {
              // check if the multiplier is a multiple of 10
              // so that we might be able to convert it to a scale that makes
              // sense
              double tmp = log10(multiplier);

              if (areApproximatelyEqual(tmp, round(tmp)))
                {
                  scale += (int)round(tmp);
                  multiplier = 1;
                }
            }

          if ((u->getExponent() == 1) &&
              areApproximatelyEqual(multiplier, 1.0) &&
              scale == 0)
            {
              result = true;
              tUnit = CModel::dimensionlessTime;
            }
          else
            {
              result = false;
            }
        }
      else
        {
          result = false;
        }
    }
  else
    {
      result = false;
    }

  return std::make_pair(tUnit, result);
}

/**
 * Returns the copasi LengthUnit corresponding to the given SBML length
 *  UnitDefinition.
 */
std::pair<CModel::LengthUnit, bool>
SBMLImporter::handleLengthUnit(const UnitDefinition* uDef)
{
  bool result = false;
  CModel::LengthUnit lUnit = CModel::m;

  if (uDef == NULL)
    {
      //DebugFile << "Argument to handleLengthUnit is NULL pointer." << std::endl;
      fatalError();
    }

  if (uDef->getNumUnits() == 1)
    {
      const Unit* u = uDef->getUnit(0);

      if (u == NULL)
        {
          //DebugFile << "Expected Unit, got NULL pointer." << std::endl;
          CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 54, "length", uDef->getId().c_str());
        }

      if ((u->getKind() == UNIT_KIND_METRE || u->getKind() == UNIT_KIND_METER) && u->getExponent() == 1)
        {
          double multiplier = u->getMultiplier();
          int scale = u->getScale();

          if (multiplier != 1)
            {
              // check if the multiplier is a multiple of 10
              // so that we might be able to convert it to a scale that makes
              // sense
              double tmp = log10(multiplier);

              if (areApproximatelyEqual(tmp, round(tmp)))
                {
                  scale += (int) round(tmp);
                  multiplier = 1;
                }
            }

          if (areApproximatelyEqual(multiplier, 1.0) &&
              ((scale % 3) == 0 || scale == -1 || scale == -2) &&
              (scale < 1) &&
              (scale > -16))
            {
              switch (scale)
                {
                  case 0:
                    lUnit = CModel::m;
                    result = true;
                    break;

                  case - 1:
                    lUnit = CModel::dm;
                    result = true;
                    break;

                  case - 2:
                    lUnit = CModel::cm;
                    result = true;
                    break;

                  case - 3:
                    lUnit = CModel::mm;
                    result = true;
                    break;

                  case - 6:
                    lUnit = CModel::microm;
                    result = true;
                    break;

                  case - 9:
                    lUnit = CModel::nm;
                    result = true;
                    break;

                  case - 12:
                    lUnit = CModel::pm;
                    result = true;
                    break;

                  case - 15:
                    lUnit = CModel::fm;
                    result = true;
                    break;

                  default:
                    //DebugFile << "Error. This value should never have been reached for the scale of the liter unit." << std::endl;
                    result = false;
                    break;
                }
            }
          else
            {
              result = false;
            }
        }
      else if ((u->getKind() == UNIT_KIND_DIMENSIONLESS))
        {
          double multiplier = u->getMultiplier();
          int scale = u->getScale();

          if (multiplier != 1)
            {
              // check if the multiplier is a multiple of 10
              // so that we might be able to convert it to a scale that makes
              // sense
              double tmp = log10(multiplier);

              if (areApproximatelyEqual(tmp, round(tmp)))
                {
                  scale += (int)round(tmp);
                  multiplier = 1;
                }
            }

          if ((u->getExponent() == 1) &&
              areApproximatelyEqual(multiplier, 1.0) &&
              scale == 0)
            {
              result = true;
              lUnit = CModel::dimensionlessLength;
            }
          else
            {
              result = false;
            }
        }
      else
        {
          result = false;
        }
    }
  else
    {
      result = false;
    }

  return std::make_pair(lUnit, result);
}

/**
 * Returns the copasi AreaUnit corresponding to the given SBML area
 *  UnitDefinition.
 */
std::pair<CModel::AreaUnit, bool>
SBMLImporter::handleAreaUnit(const UnitDefinition* uDef)
{
  bool result = false;
  CModel::AreaUnit aUnit = CModel::m2;

  if (uDef == NULL)
    {
      //DebugFile << "Argument to handleLengthUnit is NULL pointer." << std::endl;
      fatalError();
    }

  if (uDef->getNumUnits() == 1)
    {
      const Unit* u = uDef->getUnit(0);

      if (u == NULL)
        {
          //DebugFile << "Expected Unit, got NULL pointer." << std::endl;
          CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 54, "area", uDef->getId().c_str());
        }

      if ((u->getKind() == UNIT_KIND_METRE || u->getKind() == UNIT_KIND_METER) && u->getExponent() == 2)
        {
          double multiplier = u->getMultiplier();
          int scale = u->getScale();

          if (multiplier != 1)
            {
              // check if the multiplier is a multiple of 10
              // so that we might be able to convert it to a scale that makes
              // sense
              double tmp = log10(multiplier);

              if (areApproximatelyEqual(tmp, round(tmp)))
                {
                  scale += (int) round(tmp);
                  multiplier = 1;
                }
            }

          if (areApproximatelyEqual(multiplier, 1.0) &&
              ((scale % 3) == 0 || scale == -1 || scale == -2) &&
              (scale < 1) &&
              (scale > -16))
            {
              switch (scale)
                {
                  case 0:
                    aUnit = CModel::m2;
                    result = true;
                    break;

                  case - 1:
                    aUnit = CModel::dm2;
                    result = true;
                    break;

                  case - 2:
                    aUnit = CModel::cm2;
                    result = true;
                    break;

                  case - 3:
                    aUnit = CModel::mm2;
                    result = true;
                    break;

                  case - 6:
                    aUnit = CModel::microm2;
                    result = true;
                    break;

                  case - 9:
                    aUnit = CModel::nm2;
                    result = true;
                    break;

                  case - 12:
                    aUnit = CModel::pm2;
                    result = true;
                    break;

                  case - 15:
                    aUnit = CModel::fm2;
                    result = true;
                    break;

                  default:
                    //DebugFile << "Error. This value should never have been reached for the scale of the liter unit." << std::endl;
                    result = false;
                    break;
                }
            }
          else
            {
              result = false;
            }
        }
      else if ((u->getKind() == UNIT_KIND_DIMENSIONLESS))
        {
          double multiplier = u->getMultiplier();
          int scale = u->getScale();

          if (multiplier != 1)
            {
              // check if the multiplier is a multiple of 10
              // so that we might be able to convert it to a scale that makes
              // sense
              double tmp = log10(multiplier);

              if (areApproximatelyEqual(tmp, round(tmp)))
                {
                  scale += (int)round(tmp);
                  multiplier = 1;
                }
            }

          if ((u->getExponent() == 1) &&
              areApproximatelyEqual(multiplier, 1.0) &&
              scale == 0)
            {
              result = true;
              aUnit = CModel::dimensionlessArea;
            }
          else
            {
              result = false;
            }
        }
      else
        {
          result = false;
        }
    }
  else
    {
      result = false;
    }

  return std::make_pair(aUnit, result);
}

/**
 * Returns the copasi VolumeUnit corresponding to the given SBML Volume
 *  UnitDefinition.
 */
std::pair<CModel::VolumeUnit, bool>
SBMLImporter::handleVolumeUnit(const UnitDefinition* uDef)
{
  // simplify the Unitdefiniton first if this normalizes
  // the scale and the multiplier
  bool result = false;
  CModel::VolumeUnit vUnit = CModel::l;

  if (uDef == NULL)
    {
      //DebugFile << "Argument to handleVolumeUnit is NULL pointer." << std::endl;
      fatalError();
    }

  if (uDef->getNumUnits() == 1)
    {
      const Unit* u = uDef->getUnit(0);

      if (u == NULL)
        {
          //DebugFile << "Expected Unit, got NULL pointer." << std::endl;
          CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 54, "volume", uDef->getId().c_str());
        }

      if (((u->getKind() == UNIT_KIND_LITER) || (u->getKind() == UNIT_KIND_LITRE)) && u->getExponent() == 1)
        {
          double multiplier = u->getMultiplier();
          int scale = u->getScale();

          if (multiplier != 1)
            {
              // check if the multiplier is a multiple of 10
              // so that we might be able to convert it to a scale that makes
              // sense
              double tmp = log10(multiplier);

              if (areApproximatelyEqual(tmp, round(tmp)))
                {
                  scale += (int)round(tmp);
                  multiplier = 1;
                }
            }

          if (areApproximatelyEqual(multiplier, 1.0) &&
              ((scale % 3) == 0) &&
              (scale < 1) &&
              (scale > -16))
            {
              switch (scale)
                {
                  case 0:
                    vUnit = CModel::l;
                    result = true;
                    break;

                  case - 3:
                    vUnit = CModel::ml;
                    result = true;
                    break;

                  case - 6:
                    vUnit = CModel::microl;
                    result = true;
                    break;

                  case - 9:
                    vUnit = CModel::nl;
                    result = true;
                    break;

                  case - 12:
                    vUnit = CModel::pl;
                    result = true;
                    break;

                  case - 15:
                    vUnit = CModel::fl;
                    result = true;
                    break;

                  default:
                    //DebugFile << "Error. This value should never have been reached for the scale of the liter unit." << std::endl;
                    result = false;
                    break;
                }
            }
          else
            {
              result = false;
            }
        }
      else if (((u->getKind() == UNIT_KIND_METER) || (u->getKind() == UNIT_KIND_METRE)) && u->getExponent() == 3)
        {
          double multiplier = u->getMultiplier();
          int scale = u->getScale();

          if (multiplier != 1)
            {
              // check if the multiplier is a multiple of 10
              // so that we might be able to convert it to a scale that makes
              // sense
              double tmp = log10(multiplier);

              if (areApproximatelyEqual(tmp, round(tmp)))
                {
                  scale += (int)round(tmp);
                  multiplier = 1;
                }
            }

          if (areApproximatelyEqual(multiplier, 1.0) &&
              (scale == 0))
            {
              vUnit = CModel::m3;
              result = true;
            }
          else
            {
              // try to convert to liter
              Unit* pLitreUnit = convertSBMLCubicmetresToLitres(u);

              if (pLitreUnit != NULL &&
                  pLitreUnit->getExponent() == 1 &&
                  (pLitreUnit->getScale() % 3 == 0) &&
                  (pLitreUnit->getScale() < 1) &&
                  (pLitreUnit->getScale() > -16) &&
                  areApproximatelyEqual(pLitreUnit->getMultiplier(), 1.0))
                {
                  switch (pLitreUnit->getScale())
                    {
                      case 0:
                        vUnit = CModel::l;
                        result = true;
                        break;

                      case - 3:
                        vUnit = CModel::ml;
                        result = true;
                        break;

                      case - 6:
                        vUnit = CModel::microl;
                        result = true;
                        break;

                      case - 9:
                        vUnit = CModel::nl;
                        result = true;
                        break;

                      case - 12:
                        vUnit = CModel::pl;
                        result = true;
                        break;

                      case - 15:
                        vUnit = CModel::fl;
                        result = true;
                        break;

                      default:
                        result = false;
                        break;
                    }
                }
              else
                {
                  result = false;
                }

              delete pLitreUnit;
            }
        }
      else if ((u->getKind() == UNIT_KIND_DIMENSIONLESS))
        {
          double multiplier = u->getMultiplier();
          int scale = u->getScale();

          if (multiplier != 1)
            {
              // check if the multiplier is a multiple of 10
              // so that we might be able to convert it to a scale that makes
              // sense
              double tmp = log10(multiplier);

              if (areApproximatelyEqual(tmp, round(tmp)))
                {
                  scale += (int)round(tmp);
                  multiplier = 1;
                }
            }

          if ((u->getExponent() == 1) &&
              areApproximatelyEqual(multiplier, 1.0) &&
              scale == 0)
            {
              result = true;
              vUnit = CModel::dimensionlessVolume;
            }
          else
            {
              result = false;
            }
        }
      else
        {
          result = false;
        }
    }
  else
    {
      result = false;
    }

  return std::make_pair(vUnit, result);
}

CModelValue* SBMLImporter::createCModelValueFromParameter(const Parameter* sbmlParameter, CModel* copasiModel, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
{
  if (sbmlParameter->isSetUnits())
    {
      /* !!!! create a warning that the units will be ignored. */
      //CCopasiMessage Message(CCopasiMessage::WARNING, MCSBML + 26, sbmlParameter->getId().c_str());
      mIgnoredParameterUnits.push_back(sbmlParameter->getId());
      const_cast<Parameter*>(sbmlParameter)->unsetUnits();
    }

  std::string name = sbmlParameter->getName();

  if (!sbmlParameter->isSetName())
    {
      name = sbmlParameter->getId();
    }

  std::string appendix = "";
  unsigned int counter = 2;
  std::ostringstream numberStream;

  while (copasiModel->getModelValues().getIndex(name + appendix) != C_INVALID_INDEX)
    {
      numberStream.str("");
      numberStream << "_" << counter;
      counter++;
      appendix = numberStream.str();
    }

  std::string sbmlId;

  if (this->mLevel == 1)
    {
      sbmlId = sbmlParameter->getName();
    }
  else
    {
      sbmlId = sbmlParameter->getId();
    }

  CModelValue* pMV = copasiModel->createModelValue(name + appendix, 0.0);
  copasi2sbmlmap[pMV] = const_cast<Parameter*>(sbmlParameter);
  pMV->setSBMLId(sbmlId);
  SBMLImporter::importMIRIAM(sbmlParameter, pMV);
  SBMLImporter::importNotes(pMV, sbmlParameter);

  if (this->mLevel > 2)
    {
      this->mSBMLIdModelValueMap[sbmlId] = pMV;
    }

  return pMV;
}

bool SBMLImporter::sbmlId2CopasiCN(ASTNode* pNode, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap, CCopasiParameterGroup& pParamGroup)
{
  // TODO CRITICAL We need to use a node iterator

  bool success = true;
  unsigned int i, iMax = pNode->getNumChildren();

  if (pNode->getType() == AST_NAME)
    {
      Compartment* pSBMLCompartment = NULL;
      Species* pSBMLSpecies = NULL;
      Reaction* pSBMLReaction = NULL;
      Parameter* pSBMLParameter = NULL;
      std::string sbmlId;
      std::string name = pNode->getName();
      CCopasiParameter* pParam = pParamGroup.getParameter(name);

      std::map<std::string, double>::const_iterator speciesReference = mSBMLSpeciesReferenceIds.find(name);

      if (speciesReference  != mSBMLSpeciesReferenceIds.end())
        {
          // replace the name with the value
          pNode->setType(AST_REAL);
          pNode->setValue(mSBMLSpeciesReferenceIds[name]);
        }
      else if (pParam)
        {
          pNode->setName(pParam->getCN().c_str());
        }
      else
        {
          std::map<CCopasiObject*, SBase*>::iterator it = copasi2sbmlmap.begin();
          std::map<CCopasiObject*, SBase*>::iterator endIt = copasi2sbmlmap.end();
          bool found = false;

          while (it != endIt)
            {
              int type = it->second->getTypeCode();

              switch (type)
                {
                  case SBML_COMPARTMENT:
                    pSBMLCompartment = dynamic_cast<Compartment*>(it->second);

                    if (this->mLevel == 1)
                      {
                        sbmlId = pSBMLCompartment->getName();
                      }
                    else
                      {
                        sbmlId = pSBMLCompartment->getId();
                      }

                    if (sbmlId == pNode->getName())
                      {
                        pNode->setName(dynamic_cast<CCompartment*>(it->first)->getObject(CCopasiObjectName("Reference=InitialVolume"))->getCN().c_str());
                        found = true;
                      }

                    break;

                  case SBML_SPECIES:
                    pSBMLSpecies = dynamic_cast<Species*>(it->second);

                    if (this->mLevel == 1)
                      {
                        sbmlId = pSBMLSpecies->getName();
                      }
                    else
                      {
                        sbmlId = pSBMLSpecies->getId();
                      }

                    if (sbmlId == pNode->getName())
                      {
                        pNode->setName(dynamic_cast<CMetab*>(it->first)->getObject(CCopasiObjectName("Reference=InitialConcentration"))->getCN().c_str());
                        found = true;
                      }

                    break;

                  case SBML_REACTION:
                    pSBMLReaction = dynamic_cast<Reaction*>(it->second);

                    if (this->mLevel == 1)
                      {
                        sbmlId = pSBMLReaction->getName();
                      }
                    else
                      {
                        sbmlId = pSBMLReaction->getId();
                      }

                    if (sbmlId == pNode->getName())
                      {
                        pNode->setName(dynamic_cast<CReaction*>(it->first)->getObject(CCopasiObjectName("Reference=ParticleFlux"))->getCN().c_str());
                        found = true;
                      }

                    break;

                  case SBML_PARAMETER:
                    pSBMLParameter = dynamic_cast<Parameter*>(it->second);

                    if (this->mLevel == 1)
                      {
                        sbmlId = pSBMLParameter->getName();
                      }
                    else
                      {
                        sbmlId = pSBMLParameter->getId();
                      }

                    if (sbmlId == pNode->getName())
                      {
                        pNode->setName(dynamic_cast<CModelValue*>(it->first)->getValueReference()->getCN().c_str());
                        found = true;
                      }

                    break;

                  default:
                    break;
                }

              ++it;
            }

          if (!found) success = false;
        }
    }

  for (i = 0; i < iMax; ++i)
    {
      if (!this->sbmlId2CopasiCN(pNode->getChild(i), copasi2sbmlmap, pParamGroup))
        {
          success = false;
          break;
        }
    }

  return success;
}

#ifdef COPASI_DEBUG
void SBMLImporter::printMap(const std::map<CCopasiObject*, SBase*> & copasi2sbml)
{
  std::map<CCopasiObject*, SBase*>::const_iterator it = copasi2sbml.begin();
  std::map<CCopasiObject*, SBase*>::const_iterator end = copasi2sbml.end();
  std::cout << "Number of elements: " << copasi2sbml.size() << std::endl;

  while (it != end)
    {
      std::cout << "(@" << it->first << ")" << it->first->getObjectName() << " : " << "(@" << it->second << ")" << it->second->getTypeCode() << std::endl;
      ++it;
    }

  std::cout << std::endl;
}
#endif // COPASI_DEBUG

void SBMLImporter::restoreFunctionDB()
{
  // set all the old sbml ids
  std::map<CFunction*, std::string>::iterator it = this->sbmlIdMap.begin();
  std::map<CFunction*, std::string>::iterator endIt = this->sbmlIdMap.end();

  while (it != endIt)
    {
      it->first->setSBMLId(it->second);
      ++it;
    }

  // remove all the functions that were added during import
  std::set<std::string>::iterator it2 = this->mUsedFunctions.begin();
  std::set<std::string>::iterator endIt2 = this->mUsedFunctions.end();

  while (it2 != endIt2)
    {
      CEvaluationTree* pTree = this->functionDB->findFunction(*it2);
      assert(pTree);

      if (pTree->getType() == CEvaluationTree::UserDefined)
        {
          this->functionDB->removeFunction(pTree->getKey());
        }

      ++it2;
    }
}

void SBMLImporter::preprocessNode(ConverterASTNode* pNode, Model* pSBMLModel, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap, Reaction* pSBMLReaction)
{
  // this function goes through the tree several times.
  // this can probably be handled more intelligently
  // maybe we should use some kind of visitor schema
  //
  // check if there are units on pure numbers
  // so that we can create a warning that they have been ignored

  if (this->mLevel > 2 && !this->mUnitOnNumberFound)
    {
      this->mUnitOnNumberFound = SBMLImporter::checkForUnitsOnNumbers(pNode);
    }

  // first replace the calls to explicitly time dependent functions
  this->replaceTimeDependentFunctionCalls(pNode);

  if (!this->mDelayFound || pSBMLReaction != NULL)
    {
      bool result = isDelayFunctionUsed(pNode);

      if (pSBMLReaction != NULL && result)
        {
          // if this is the first delay we find, we have to populate the id set
          // in order to be able to create new global parameters with unique ids
          if (!this->mUsedSBMLIdsPopulated)
            {
              std::map<std::string, const SBase*> idMap;
              std::map<std::string, const SBase*> metaIdMap;
              // this is overkill, but better than writing yet another function to
              // collect ids
              SBMLUtils::collectIds(pSBMLModel, idMap, metaIdMap);
              std::map<std::string, const SBase*>::iterator it = idMap.begin(), endit = idMap.end();

              while (it != endit)
                {
                  this->mUsedSBMLIds.insert(it->first);
                  ++it;
                }

              this->mUsedSBMLIdsPopulated = true;
              CCopasiMessage(CCopasiMessage::WARNING, MCSBML + 86);
            }

          // we need a map to store the replacements for local parameters
          // which occur within delay calls
          std::map<std::string, std::string> replacementMap;
          this->replace_delay_nodes(pNode, pSBMLModel, copasi2sbmlmap, pSBMLReaction, replacementMap);

          if (!replacementMap.empty())
            {
              // replace all local parameter nodes which have been converted to global parameters
              // and that were not in delay calls
              this->replace_name_nodes(pNode, replacementMap);
              // delete the local parameters that have been replaced from the
              // reaction
              std::map<std::string, std::string>::const_iterator it = replacementMap.begin(), endit = replacementMap.end();
              // starting with SBML Level 3, the parameters of a kinetic law are expressed in
              // terms of a new class called LocalParameter instead of Parameter
              // Unfortunatelly libsbml 4.1 uses separate data structures for
              // the Parameters and the LocalParameters which mandates a small
              // code change to be on the safe side
              ListOfParameters* pList = NULL;

              if (this->mLevel > 2)
                {
                  pList = pSBMLReaction->getKineticLaw()->getListOfLocalParameters();
                }
              else
                {
                  pList = pSBMLReaction->getKineticLaw()->getListOfParameters();
                }

              Parameter* pParam = NULL;

              while (it != endit)
                {
                  pParam = pList->remove(it->first);
                  assert(pParam != NULL);
                  pdelete(pParam);
                  ++it;
                }

              this->mReactionsWithReplacedLocalParameters.insert(pSBMLReaction->getId());
            }
        }

      this->mDelayFound = result;
    }

  this->replaceCallNodeNames(pNode);
  this->replaceTimeAndAvogadroNodeNames(pNode);

  if (pSBMLReaction != NULL && !this->mSubstanceOnlySpecies.empty())
    {
      this->multiplySubstanceOnlySpeciesByVolume(pNode);
    }

  if (!this->mSubstanceOnlySpecies.empty() && this->mpCopasiModel->getQuantityUnitEnum() != CModel::number && pSBMLReaction == NULL)
    {
      this->replaceAmountReferences(pNode, pSBMLModel, this->mpCopasiModel->getQuantity2NumberFactor(), copasi2sbmlmap);
    }
}

/**
 * This method replaces references to the id of species which have the
 * hasOnlySubstanceUnits flag set with the reference divided by avogadros
 * number.
 */
void SBMLImporter::replaceAmountReferences(ConverterASTNode* pASTNode, Model* pSBMLModel, double factor, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
{
  CNodeIterator< ConverterASTNode > itNode(pASTNode);
  itNode.setProcessingModes(CNodeIteratorMode::Before);

  while (itNode.next() != itNode.end())
    {
      if (*itNode == NULL)
        {
          continue;
        }

      // check if there is already a multiplication with the avogadro number
      // if so, we just replace that multiplication instead of adding an additional division
      //
      // we only do this, if there are two children to the multiplication because otherwise we would have to do sophisticated checks
      // whether the multiplication is really only between a certain species and the avogadro number
      if (itNode->getType() == AST_TIMES && itNode->getNumChildren() == 2)
        {
          std::set<const Parameter*>::const_iterator avoIt = this->mPotentialAvogadroNumbers.begin();
          std::set<const Parameter*>::const_iterator avoEndit = this->mPotentialAvogadroNumbers.end();
          // check if one of the potantial avogrador numbers is a child to this multiplication
          ASTNode* pChild1 = itNode->getChild(0);
          ASTNode* pChild2 = itNode->getChild(1);
          assert(pChild1 != NULL);
          assert(pChild2 != NULL);

          if (pChild1->getType() == AST_NAME && pChild2->getType() == AST_NAME)
            {
              std::string id1 = pChild1->getName();
              std::string id2 = pChild2->getName();
              ASTNode *pAvogadro = NULL, *pSpecies = NULL;

              if (pChild1 != NULL && pChild2 != NULL)
                {
                  while (avoIt != avoEndit)
                    {
                      if (id1 == (*avoIt)->getId())
                        {
                          pAvogadro = pChild1;
                          pSpecies = pChild2;
                          // store the id of the species in id1 for use below
                          id1 = id2;
                          break;
                        }
                      else if (id2 == (*avoIt)->getId())
                        {
                          pAvogadro = pChild2;
                          pSpecies = pChild1;
                          break;
                        }

                      ++avoIt;
                    }

                  if (pAvogadro != NULL)
                    {
                      // check if the potential speices node is really a substance only species
                      //
                      std::map<Species*, Compartment*>::const_iterator it = this->mSubstanceOnlySpecies.begin();
                      std::map<Species*, Compartment*>::const_iterator endit = this->mSubstanceOnlySpecies.end();

                      while (it != endit)
                        {
                          if (it->first->getId() == id1)
                            {
                              // delete the two children
                              itNode->removeChild(1);
                              itNode->removeChild(0);
                              pdelete(pChild1);
                              pdelete(pChild2);
                              // now we know that we can change the current node
                              // to represent the species
                              itNode->setType(AST_NAME);
                              itNode->setName(id1.c_str());
                              itNode.skipChildren();
                              break;
                            }

                          ++it;
                        }
                    }
                }
            }
          else if ((pChild1->getType() == AST_NAME && (pChild2->getType() == AST_REAL || pChild2->getType() == AST_REAL_E)) ||
                   (pChild2->getType() == AST_NAME && (pChild1->getType() == AST_REAL || pChild1->getType() == AST_REAL_E)))
            {
              double value = 0.0;
              std::string id;

              if (pChild1->getType() == AST_NAME)
                {
                  value = pChild2->getMantissa() * pow(10.0, (double)pChild2->getExponent());
                  id = pChild1->getName();
                }
              else
                {
                  value = pChild1->getMantissa() * pow(10.0, (double)pChild1->getExponent());
                  id = pChild2->getName();
                }

              if (areApproximatelyEqual(factor, value, 1e-3))
                {
                  std::map<Species*, Compartment*>::const_iterator it = this->mSubstanceOnlySpecies.begin();
                  std::map<Species*, Compartment*>::const_iterator endit = this->mSubstanceOnlySpecies.end();

                  while (it != endit)
                    {
                      if (it->first->getId() == id)
                        {
                          // delete the two children
                          itNode->removeChild(1);
                          itNode->removeChild(0);
                          pdelete(pChild1);
                          pdelete(pChild2);
                          // now we know that we can change the current node
                          // to represent the species
                          itNode->setType(AST_NAME);
                          itNode->setName(id.c_str());
                          itNode.skipChildren();
                          break;
                        }

                      ++it;
                    }
                }
            }
        }
      else if (itNode->getType() == AST_NAME)
        {
          // check if pNode is a reference to a hasOnlySubstance species
          std::string id = itNode->getName();
          std::map<Species*, Compartment*>::const_iterator it = this->mSubstanceOnlySpecies.begin();
          std::map<Species*, Compartment*>::const_iterator endit = this->mSubstanceOnlySpecies.end();

          while (it != endit)
            {
              if (it->first->getId() == id)
                {
                  break;
                }

              ++it;
            }

          if (it != endit)
            {
              // replace pNode by a division by the quantity to number
              // factor
              if (this->mPotentialAvogadroNumbers.empty())
                {
                  this->createHasOnlySubstanceUnitFactor(pSBMLModel, factor, copasi2sbmlmap);
                }

              itNode->setType(AST_DIVIDE);
              itNode->setCharacter('/');
              ConverterASTNode* pChild = new ConverterASTNode(AST_NAME);
              pChild->setName(id.c_str());
              itNode->addChild(pChild);
              id = (*this->mPotentialAvogadroNumbers.begin())->getId();
              pChild = new ConverterASTNode(AST_NAME);
              pChild->setName(id.c_str());
              itNode->addChild(pChild);
            }
        }
    }
}

bool SBMLImporter::isDelayFunctionUsed(ConverterASTNode* pASTNode)
{
  bool result = false;
  CNodeIterator< ConverterASTNode > itNode(pASTNode);

  while (itNode.next() != itNode.end())
    {
      if (*itNode == NULL)
        {
          continue;
        }

      if (itNode->getType() == AST_FUNCTION_DELAY)
        {
          result = true;
          break;
        }
    }

  return result;
}

void SBMLImporter::replaceTimeAndAvogadroNodeNames(ASTNode* pASTNode)
{
  CNodeIterator< ASTNode > itNode(pASTNode);

  while (itNode.next() != itNode.end())
    {
      if (*itNode == NULL)
        {
          continue;
        }

      if (itNode->getType() == AST_NAME_TIME)
        {
          itNode->setName(this->mpCopasiModel->getObject(CCopasiObjectName("Reference=Time"))->getCN().c_str());
        }
      else if (itNode->getType() == AST_NAME_AVOGADRO)
        {
          itNode->setName(this->mpCopasiModel->getObject(CCopasiObjectName("Reference=Avogadro Constant"))->getCN().c_str());
        }
    }
}

void SBMLImporter::replaceCallNodeNames(ASTNode* pASTNode)
{
  CNodeIterator< ASTNode > itNode(pASTNode);

  while (itNode.next() != itNode.end())
    {
      if (*itNode == NULL)
        {
          continue;
        }

      if (itNode->getType() == AST_FUNCTION)
        {
          std::map<std::string, std::string>::const_iterator pos = this->mFunctionNameMapping.find(itNode->getName());

          std::map<std::string, std::string>::const_iterator knownPos = mKnownCustomUserDefinedFunctions.find(itNode->getName());

          if (pos == this->mFunctionNameMapping.end())
            {
              CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 47, itNode->getName());
            }

          if (knownPos != mKnownCustomUserDefinedFunctions.end() &&
              (knownPos->second == "RATE" && itNode->getNumChildren() == 1))
            {
              // replace the known function with the correct calls.
              // replace rateOf with the rate of whatever is inside
              std::string symbol = itNode->getChild(0)->getName();
              itNode->removeChild(0);
              itNode->setType(AST_NAME);
              itNode->setName(symbol.c_str());
              itNode->setUserData(strdup("RATE"));
            }
          else
            {
              std::string newName = pos->second;
              itNode->setName(newName.c_str());
              this->mUsedFunctions.insert(newName);
            }
        }
    }
}

/**
 * The methods gets a function where all the parameters have a usage of "PARAMETER".
 * In addition it get the root node of a call to that function which is an expression
 * and contains the actual objects with which the function is called in a certain reaction.
 * From this expression we can determine if there already is a function in the database
 * that does the same. Or we can find out if this function is a Mass Action kinetic.
 */
CFunction* SBMLImporter::findCorrespondingFunction(const CExpression * pExpression, const CReaction* pCopasiReaction)
{
  // We are certain at this point that pExpression is a simple function call.

  // Check whether a function with the name exists in COPASI and whether it is suitable
  CFunction* pCorrespondingFunction = NULL;

  std::string Name = pExpression->getRoot()->getData();
  pCorrespondingFunction = functionDB->findFunction(Name);

  if (pCorrespondingFunction != NULL)
    {
      const CFunctionParameters & Variables = pCorrespondingFunction->getVariables();

      for (size_t i = 0; i < Variables.size(); ++i)
        {
          // This fails for global parameters but those are handled in a second attempt
          if (pCopasiReaction->getParameterIndex(Variables[i]->getObjectName()) == C_INVALID_INDEX)
            {
              pCorrespondingFunction = NULL;
              break;
            }
        }
    }

  if (pCorrespondingFunction != NULL)
    {
      return pCorrespondingFunction;
    }

  // TODO This has never worked
#ifdef XXXX
  // Check for suitable functions which have a different name but have the same call tree as
  std::vector< CFunction * > functions = functionDB->suitableFunctions(pCopasiReaction->getChemEq().getSubstrates().size(),
                                         pCopasiReaction->getChemEq().getProducts().size(),
                                         pCopasiReaction->isReversible() ? TriTrue : TriFalse);
  size_t i, iMax = functions.size();

  for (i = 0; i < iMax; ++i)
    {
      pCorrespondingFunction = (functions[i]);

      // make sure the function is not compared to itself since it can already
      // be in the database if it has been used a call in another function
      // don't compare the mass action kinetics
      if ((pCorrespondingFunction != tree) && (!dynamic_cast<CMassAction*>(pCorrespondingFunction)) && areEqualFunctions(pCorrespondingFunction, tree))
        {
          const CFunctionParameters & Variables = pFun->getVariables();

          for (size_t i = 0; i < Variables.size(); ++i)
            {
              if (pCopasiReaction->getParameterIndex(Variables[i]->getObjectName()) == C_INVALID_INDEX)
                {
                  pCorrespondingFunction = NULL;
                  break;
                }
            }

          pCorrespondingFunction = pFun;
          break;
        }
    }

#endif // XXXX

  return pCorrespondingFunction;
}

// static
bool SBMLImporter::areEqualFunctions(const CFunction* pFun, const CFunction* pFun2)
{
  bool result = true;
  const CFunctionParameters& funParams1 = pFun->getVariables();
  const CFunctionParameters& funParams2 = pFun2->getVariables();

  if (funParams1.size() == funParams2.size())
    {
      size_t i, iMax = funParams1.size();

      for (i = 0; i < iMax; ++i)
        {
          const CFunctionParameter* pFunParam1 = funParams1[i];
          const CFunctionParameter* pFunParam2 = funParams2[i];

          if (pFunParam1->getObjectName() != pFunParam2->getObjectName())
            {
              result = false;
              break;
            }
        }

      if (result == true)
        {
          const CEvaluationNode* pNodeFun1 = static_cast<const CEvaluationNode*>(pFun->getRoot());
          const CEvaluationNode* pNodeFun2 = static_cast<const CEvaluationNode*>(pFun2->getRoot());
          result = areEqualSubtrees(pNodeFun1, pNodeFun2);
        }
    }
  else
    {
      result = false;
    }

  return result;
}

// static
bool SBMLImporter::areEqualSubtrees(const CEvaluationNode* pNode1, const CEvaluationNode* pNode2)
{
  // TODO CRITICAL We need to use a node iterator
  bool result = ((pNode1->getType() == pNode2->getType()) && (pNode1->getData() == pNode2->getData()));
  const CEvaluationNode* pChild1 = static_cast<const CEvaluationNode*>(pNode1->getChild());
  const CEvaluationNode* pChild2 = static_cast<const CEvaluationNode*>(pNode2->getChild());

  while (result && pChild1 && pChild2)
    {
      result = areEqualSubtrees(pChild1, pChild2);
      pChild1 = static_cast<const CEvaluationNode*>(pChild1->getSibling());
      pChild2 = static_cast<const CEvaluationNode*>(pChild2->getSibling());
    }

  result = (result && !pChild1 && !pChild2);
  return result;
}

std::vector<CEvaluationNodeObject*>* SBMLImporter::isMassAction(const CEvaluationTree* pTree, const CChemEq& chemicalEquation, const CEvaluationNodeCall* pCallNode)
{
  // TODO CRITICAL We need to use a node iterator
  CEvaluationTree::Type type = pTree->getType();
  std::vector< std::vector< std::string > > functionArgumentCNs;
  const CEvaluationNode* pChildNode = NULL;
  std::string str;
  std::vector<CEvaluationNodeObject*>* result = NULL;

  switch (type)
    {
      case CEvaluationTree::Function:
      case CEvaluationTree::UserDefined:
      case CEvaluationTree::PreDefined:
        pChildNode = static_cast<const CEvaluationNode*>(pCallNode->getChild());

        while (pChildNode)
          {
            if (pChildNode->getType() == CEvaluationNode::OBJECT)
              {
                str = pChildNode->getData().substr(1, pChildNode->getData().length() - 2);
                functionArgumentCNs.push_back(std::vector<std::string>());
                functionArgumentCNs[functionArgumentCNs.size() - 1].push_back(str);
                pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
              }
            else
              {
                fatalError();
              }
          }

        result = this->isMassActionFunction(dynamic_cast<const CFunction*>(pTree), chemicalEquation, functionArgumentCNs);
        break;

      case CEvaluationTree::Expression:
        result = this->isMassActionExpression(pTree->getRoot(), chemicalEquation);
        break;

      default:
        fatalError();
        break;
    }

  return result;
}

std::vector<CEvaluationNodeObject*>* SBMLImporter::isMassActionExpression(const CEvaluationNode* pRootNode, const CChemEq& chemicalEquation)
{
  bool result = true;
  std::vector<CEvaluationNodeObject*>* v = NULL;

  if (chemicalEquation.getReversibility())
    {
#ifdef COPASI_DEBUG
      CEvaluationNode* pTmpNode = pRootNode->copyBranch();
//      CEvaluationNode* pTmpNode = CEvaluationNodeNormalizer::normalize(pRootNode);
#else
      CEvaluationNode* pTmpNode = pRootNode->copyBranch();
#endif /* COPASI_DEBUG */
      assert(pTmpNode != NULL);
      // the root node must be a minus operator
      // the two children must be irreversible mass action terms
      result = (CEvaluationNode::type(pTmpNode->getType()) == CEvaluationNode::OPERATOR && (CEvaluationNodeOperator::SubType)(CEvaluationNode::subType(pTmpNode->getType())) == CEvaluationNodeOperator::MINUS);

      if (result)
        {
          const CEvaluationNode* pChildNode = static_cast<const CEvaluationNode*>(pTmpNode->getChild());
          result = (pChildNode != NULL);

          if (result)
            {
              CChemEq tmpEq;
              const CCopasiVector<CChemEqElement>* metabolites = &chemicalEquation.getSubstrates();
              size_t i, iMax = metabolites->size();
              result = (iMax > 0);

              if (result)
                {
                  for (i = 0; i < iMax; ++i)
                    {
                      const CChemEqElement* element = (*metabolites)[i];
                      tmpEq.addMetabolite(element->getMetaboliteKey(), element->getMultiplicity(), CChemEq::SUBSTRATE);
                    }

                  v = this->isMassActionExpression(pChildNode, tmpEq);

                  if (!v)
                    {
                      fatalError();
                    }

                  result = !v->empty();

                  if (result)
                    {
                      pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
                      result = (pChildNode != NULL);

                      if (result)
                        {
                          CChemEq tmpEq2;
                          metabolites = &chemicalEquation.getProducts();
                          iMax = metabolites->size();
                          result = (iMax > 0);

                          for (i = 0; i < iMax; ++i)
                            {
                              const CChemEqElement* element = (*metabolites)[i];
                              tmpEq2.addMetabolite(element->getMetaboliteKey(), element->getMultiplicity(), CChemEq::SUBSTRATE);
                            }

                          std::vector<CEvaluationNodeObject*>* v2 = this->isMassActionExpression(pChildNode, tmpEq2);

                          if (!v2)
                            {
                              fatalError();
                            }

                          result = !v2->empty();

                          if (result)
                            {
                              v->push_back((*v2)[0]);
                            }
                          else
                            {
                              v->clear();
                            }

                          pdelete(v2);
                        }
                    }
                }
              else
                {
                  v = new std::vector<CEvaluationNodeObject*>;
                }
            }
        }
      else
        {
          v = new std::vector<CEvaluationNodeObject*>;
        }

      pdelete(pTmpNode);
    }
  else
    {
      // the expression must contain exactly one global or local parameter
      // the expression must contain each substrate in the CChemicalReaction
      std::vector<const CEvaluationNode*> arguments;
      std::map<const CMetab*, C_FLOAT64> multiplicityMap;
      this->separateProductArguments(pRootNode, arguments);
      size_t numParameters = 0, i, iMax = arguments.size();
      v = new std::vector<CEvaluationNodeObject*>;

      if (iMax != 0)
        {
          std::vector<CCopasiContainer*> listOfContainers;
          listOfContainers.push_back(this->mpCopasiModel);

          for (i = 0; (i < iMax) && (numParameters < 2); ++i)
            {
              const CEvaluationNode* pNode = arguments[i];

              // the node can either be an object node
              // or it can be a power function node
              if (pNode->getType() == CEvaluationNode::OBJECT)
                {
                  // it can be a global or a local parameter or an metabolite
                  std::string objectCN = pNode->getData().substr(1, pNode->getData().length() - 2);
                  const CCopasiObject* pObject = mpDataModel->ObjectFromName(listOfContainers, objectCN);

                  if (!pObject)
                    {
                      CCopasiMessage(CCopasiMessage::EXCEPTION, MCSBML + 39, objectCN.c_str());
                    }

                  if (pObject->isReference())
                    {
                      pObject = pObject->getObjectParent();

                      if (dynamic_cast<const CMetab*>(pObject))
                        {
                          const CMetab* pMetab = static_cast<const CMetab*>(pObject);

                          if (multiplicityMap.find(pMetab) != multiplicityMap.end())
                            {
                              multiplicityMap[pMetab] = multiplicityMap[pMetab] + 1.0;
                            }
                          else
                            {
                              multiplicityMap[pMetab] = 1.0;
                            }
                        }
                      else if (dynamic_cast<const CModelValue*>(pObject))
                        {
                          ++numParameters;
                          v->push_back(new CEvaluationNodeObject((CEvaluationNodeObject::SubType)(CEvaluationNode::subType(pNode->getType())), pNode->getData()));
                        }
                      else
                        {
                          result = false;
                          break;
                        }
                    }
                  else if (dynamic_cast<const CCopasiParameter*>(pObject))
                    {
                      ++numParameters;
                      v->push_back(new CEvaluationNodeObject((CEvaluationNodeObject::SubType)(CEvaluationNode::subType(pNode->getType())), pNode->getData()));
                    }
                  else
                    {
                      result = false;
                      break;
                    }
                }
              else if (pNode->getType() == CEvaluationNode::OPERATOR && (CEvaluationNodeOperator::SubType)(CEvaluationNode::subType(pNode->getType())) == CEvaluationNodeOperator::POWER)
                {
                  // the two children must be a metabolite node and a number node in this order
                  const CEvaluationNode* pChildNode = static_cast<const CEvaluationNode*>(pNode->getChild());

                  if (pChildNode->getType() == CEvaluationNode::OBJECT)
                    {
                      std::string objectCN = pChildNode->getData().substr(1, pChildNode->getData().length() - 2);
                      const CCopasiObject* pObject = mpDataModel->ObjectFromName(listOfContainers, objectCN);
                      assert(pObject);

                      if (pObject->isReference())
                        {
                          pObject = pObject->getObjectParent();

                          if (dynamic_cast<const CMetab*>(pObject))
                            {
                              pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
                              assert(pChildNode);

                              CEvaluationNode::Type type = CEvaluationNode::type(pChildNode->getType());

                              if (type == CEvaluationNode::NUMBER)
                                {
                                  const CMetab* pMetab = static_cast<const CMetab*>(pObject);

                                  if (multiplicityMap.find(pMetab) != multiplicityMap.end())
                                    {
                                      multiplicityMap[pMetab] = multiplicityMap[pMetab] + pChildNode->getValue();
                                    }
                                  else
                                    {
                                      multiplicityMap[pMetab] = pChildNode->getValue();
                                    }
                                }
                              else if (type == CEvaluationNode::FUNCTION &&
                                       (CEvaluationNodeFunction::SubType)CEvaluationNode::subType(pChildNode->getType()) == CEvaluationNodeFunction::MINUS &&
                                       CEvaluationNode::type(((CEvaluationNode*)pChildNode->getChild())->getType()) == CEvaluationNode::NUMBER)
                                {
                                  const CMetab* pMetab = static_cast<const CMetab*>(pObject);
                                  multiplicityMap[pMetab] = -1 * ((CEvaluationNodeNumber*)(pChildNode->getChild()))->getValue();
                                }
                              else
                                {
                                  // not a math action
                                }
                            }
                        }
                      else
                        {
                          result = false;
                          break;
                        }
                    }
                  else
                    {
                      result = false;
                      break;
                    }
                }
              else
                {
                  result = false;
                  break;
                }
            }
        }

      if (numParameters != 1)
        {
          result = false;
        }

      if (result)
        {
          const CCopasiVector<CChemEqElement>& metabolites = chemicalEquation.getSubstrates();
          size_t i, iMax = metabolites.size();

          // all metabolites must occur in the muliplicityMap so they have to have the same size
          // and a mass action must have at least one metabolite
          if (iMax == 0 || iMax != multiplicityMap.size()) result = false;

          for (i = 0; i < iMax && result; ++i)
            {
              // the metabolite has to be present in the multiplicityMap, otherwise it is not a mass action
              // the stoichiometry also has to fit
              std::map<const CMetab*, C_FLOAT64>::iterator pos = multiplicityMap.find(metabolites[i]->getMetabolite());

              if (pos == multiplicityMap.end() ||
                  !areApproximatelyEqual(pos->second, metabolites[i]->getMultiplicity(), 0.01))
                {
                  result = false;
                  break;
                }
            }
        }

      if (!result)
        {
          v->clear();
        }
    }

  return v;
}

std::vector<CEvaluationNodeObject*>* SBMLImporter::isMassActionFunction(const CFunction* pFun, const CChemEq& chemicalEquation, const std::vector<std::vector< std::string > >& functionArgumentCNs)
{
  // create an expression from the function and call isMassActionExpression
  CEvaluationTree* pExpressionTree = this->createExpressionFromFunction(pFun, functionArgumentCNs);

  if (!pExpressionTree)
    {
      return NULL;
    }

  std::vector<CEvaluationNodeObject*>* v = this->isMassActionExpression(pExpressionTree->getRoot(), chemicalEquation);

  delete pExpressionTree;
  return v;
}

CEvaluationTree* SBMLImporter::createExpressionFromFunction(const CFunction* pFun, const std::vector<std::vector<std::string > >& functionArgumentCNs)
{
  CEvaluationTree* pTree = NULL;
  const CFunctionParameters& pFunParams = pFun->getVariables();
  std::string str;

  if (pFunParams.size() == functionArgumentCNs.size())
    {
      std::map<std::string , std::string> variable2CNMap;
      size_t i, iMax = pFunParams.size();

      for (i = 0; i < iMax; ++i)
        {
          // vectors should not occur here
          assert(functionArgumentCNs[i].size() == 1);
          variable2CNMap[pFunParams[i]->getObjectName()] = functionArgumentCNs[i][0];
        }

      CEvaluationNode* pTmpNode = this->variables2objects(pFun->getRoot(), variable2CNMap);
      assert(pTmpNode);
      pTree = CEvaluationTree::create(CEvaluationTree::Expression);
      pTree->setRoot(pTmpNode);
    }

  return pTree;
}

void SBMLImporter::separateProductArguments(const CEvaluationNode* pRootNode, std::vector<const CEvaluationNode*>& arguments)
{
  // TODO CRITICAL We need to use a node iterator
  const CEvaluationNodeOperator* pMultiplyNode = dynamic_cast<const CEvaluationNodeOperator*>(pRootNode);

  if (pMultiplyNode && (((CEvaluationNodeOperator::SubType)CEvaluationNode::subType(pMultiplyNode->getType())) == CEvaluationNodeOperator::MULTIPLY))
    {
      // check if one if the children is an object node or a power operator, if so,
      // add the node to the vector
      // the nodes not one of those two are passed to this function recursively.
      const CEvaluationNode* pChildNode = static_cast<const CEvaluationNode*>(pMultiplyNode->getChild());

      while (pChildNode)
        {
          const CEvaluationNodeObject* pObjectNode = dynamic_cast<const CEvaluationNodeObject*>(pChildNode);
          const CEvaluationNodeOperator* pOperatorNode = dynamic_cast<const CEvaluationNodeOperator*>(pChildNode);

          if (pObjectNode)
            {
              arguments.push_back(pObjectNode);
            }
          else if (pOperatorNode && (((CEvaluationNodeOperator::SubType)CEvaluationNode::subType(pOperatorNode->getType())) == CEvaluationNodeOperator::POWER))
            {
              arguments.push_back(pOperatorNode);
            }
          else
            {
              this->separateProductArguments(pChildNode, arguments);

              if (arguments.empty())
                {
                  // it is not a mass action kinetic, so we can stop here
                  break;
                }
            }

          pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
        }
    }
  else
    {
      arguments.clear();
    }
}

void SBMLImporter::setCorrectUsage(CReaction* pCopasiReaction, const CEvaluationNodeCall* pCallNode)
{
  // find out what type each argument in the call node has.
  // it can be a local parameter, a global parameter, a compartment or a metabolite
  // if it is a metabolite, try to find out if it is a substrate, product, or modifier
  if (!pCallNode)
    {
      fatalError();
    }

  const CEvaluationNode* pChildNode = static_cast<const CEvaluationNode*>(pCallNode->getChild());

  std::vector<CCopasiContainer*> listOfContainers = std::vector<CCopasiContainer*>();

  listOfContainers.push_back(this->mpCopasiModel);

  CChemEq& pChemEq = pCopasiReaction->getChemEq();

  std::map<const CChemEqElement*, CChemEq::MetaboliteRole> v;

  const CCopasiVector<CChemEqElement>* pV = &pChemEq.getSubstrates();

  size_t i, iMax = pV->size();

  for (i = 0; i < iMax; ++i)
    {
      v[(*pV)[i]] = CChemEq::SUBSTRATE;
    }

  pV = &pChemEq.getProducts();
  iMax = pV->size();

  for (i = 0; i < iMax; ++i)
    {
      v[(*pV)[i]] = CChemEq::PRODUCT;
    }

  pV = &pChemEq.getModifiers();
  iMax = pV->size();

  for (i = 0; i < iMax; ++i)
    {
      v[(*pV)[i]] = CChemEq::MODIFIER;
    }

  size_t parameterIndex = 0;

  while (pChildNode)
    {
      const CEvaluationNodeObject* pObjectNode = dynamic_cast<const CEvaluationNodeObject*>(pChildNode);

      if (!pObjectNode)
        {
          fatalError();
        }

      const CCopasiObject* object = mpDataModel->ObjectFromName(listOfContainers, pObjectNode->getData().substr(1, pObjectNode->getData().length() - 2));

      if (!object)
        {
          fatalError();
        }

      CFunctionParameter* pFunParam = const_cast<CFunctionParameter*>(pCopasiReaction->getFunction()->getVariables()[parameterIndex]);

      if (dynamic_cast<const CCopasiObjectReference<C_FLOAT64>*>(object))
        {
          object = object->getObjectParent();

          if (dynamic_cast<const CMetab*>(object))
            {
              std::map<const CChemEqElement*, CChemEq::MetaboliteRole>::iterator it = v.begin();
              std::map<const CChemEqElement*, CChemEq::MetaboliteRole>::iterator endIt = v.end();

              while (it != endIt)
                {
                  if (it->first == object)
                    {
                      // get the role of the metabolite
                      switch (it->second)
                        {
                          case CChemEq::SUBSTRATE:
                            // it is a substrate
                            pFunParam->setUsage(CFunctionParameter::SUBSTRATE);
                            break;

                          case CChemEq::PRODUCT:
                            // it is a product
                            pFunParam->setUsage(CFunctionParameter::PRODUCT);
                            break;

                          case CChemEq::MODIFIER:
                            // it is a modifier
                            pFunParam->setUsage(CFunctionParameter::MODIFIER);
                            break;

                          default:
                            fatalError();
                            break;
                        }
                    }

                  ++it;
                }

              if (it == endIt)
                {
                  fatalError();
                }
            }
          else if (dynamic_cast<const CModelValue*>(object))
            {
              // it is a global parameter
              pFunParam->setUsage(CFunctionParameter::PARAMETER);
            }
          else if (dynamic_cast<const CCompartment*>(object))
            {
              // it is a volume
              pFunParam->setUsage(CFunctionParameter::VOLUME);
            }
          else
            {
              fatalError()
            }
        }
      else
        {
          // it is a local parameter
          pFunParam->setUsage(CFunctionParameter::PARAMETER);
        }

      pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
      ++parameterIndex;
    }
}

bool containsKey(const CCopasiVector < CChemEqElement >& list, const std::string& key)
{
  CCopasiVector < CChemEqElement >::const_iterator it = list.begin();
  CCopasiVector < CChemEqElement >::const_iterator end = list.end();

  for (; it != end; ++it)
    {
      if ((*it)->getMetaboliteKey() == key)
        {
          return true;
        }
    }

  return false;
}

void SBMLImporter::doMapping(CReaction* pCopasiReaction, const CEvaluationNodeCall* pCallNode)
{
  // map the first argument of the call node to the first variable of the function of the reaction
  // and so on
  if (!pCallNode)
    {
      fatalError();
    }

  std::vector<CCopasiContainer*> listOfContainers;
  listOfContainers.push_back(this->mpCopasiModel);

  if (dynamic_cast<const CMassAction*>(pCopasiReaction->getFunction()))
    {
      const CEvaluationNodeObject* pChild = dynamic_cast<const CEvaluationNodeObject*>(pCallNode->getChild());
      std::string objectCN = pChild->getData();
      objectCN = objectCN.substr(1, objectCN.length() - 2);
      CCopasiObject* pObject = mpDataModel->ObjectFromName(listOfContainers, objectCN);

      if (!pObject)
        {
          fatalError();
        }

      if (pObject->isReference())
        {
          pObject = pObject->getObjectParent();
        }

      const std::string& objectKey = pObject->getKey();

      pCopasiReaction->setParameterMapping("k1", objectKey);

      const CCopasiVector<CChemEqElement>* metabolites = &pCopasiReaction->getChemEq().getSubstrates();

      size_t i, iMax = metabolites->size();

      size_t j, jMax;

      for (i = 0; i < iMax; ++i)
        for (j = 0, jMax = static_cast<int>(fabs((*metabolites)[i]->getMultiplicity())); j < jMax; j++)
          pCopasiReaction->addParameterMapping("substrate", (*metabolites)[i]->getMetaboliteKey());

      if (pCopasiReaction->isReversible())
        {
          pChild = dynamic_cast<const CEvaluationNodeObject*>(pChild->getSibling());
          std::string objectCN = pChild->getData();
          objectCN = objectCN.substr(1, objectCN.length() - 2);
          CCopasiObject* pObject = mpDataModel->ObjectFromName(listOfContainers, objectCN);

          if (!pObject)
            {
              fatalError();
            }

          if (pObject->isReference())
            {
              pObject = pObject->getObjectParent();
            }

          const std::string& objectKey = pObject->getKey();

          pCopasiReaction->setParameterMapping("k2", objectKey);

          const CCopasiVector<CChemEqElement>* metabolites = &pCopasiReaction->getChemEq().getProducts();

          iMax = metabolites->size();

          for (i = 0; i < iMax; ++i)
            for (j = 0, jMax = static_cast<int>(fabs((*metabolites)[i]->getMultiplicity())); j < jMax; j++)
              pCopasiReaction->addParameterMapping("product", (*metabolites)[i]->getMetaboliteKey());
        }
    }
  else
    {
      const CFunctionParameters & Variables = pCopasiReaction->getFunction()->getVariables();
      size_t i, iMax = Variables.size();
      const CEvaluationNodeObject* pChild = dynamic_cast<const CEvaluationNodeObject*>(pCallNode->getChild());

      for (i = 0; i < iMax; ++i)
        {
          if (!pChild)
            {
              fatalError();
            }

          std::string objectCN = pChild->getData();
          objectCN = objectCN.substr(1, objectCN.length() - 2);
          CCopasiObject* pObject = mpDataModel->ObjectFromName(listOfContainers, objectCN);

          if (!pObject)
            {
              fatalError();
            }

          if (pObject->isReference())
            {
              pObject = pObject->getObjectParent();
            }

          const std::string& objectKey = pObject->getKey();

          pCopasiReaction->setParameterMapping(i, objectKey);

          const CChemEq& eqn = pCopasiReaction->getChemEq();

          bool reversible = eqn.getReversibility();

          // We guess what the role of a variable of newly imported function is:
          if (Variables[i]->getUsage() == CFunctionParameter::VARIABLE)
            {
              CFunctionParameter::Role Role = CFunctionParameter::PARAMETER;

              if (pObject->getObjectType() == "Metabolite")
                {
                  if (containsKey(eqn.getSubstrates(), objectKey))
                    Role = CFunctionParameter::SUBSTRATE;

                  if (Role == CFunctionParameter::PARAMETER &&
                      containsKey(eqn.getProducts(), objectKey))
                    {
                      // Fix for Bug 1882: if not reversible, only mark as product if it does
                      // not appear in the list of modifiers
                      if (!reversible && containsKey(eqn.getModifiers(), objectKey))
                        Role =  CFunctionParameter::MODIFIER;
                      else
                        Role =  CFunctionParameter::PRODUCT;
                    }

                  // It is not a substrate and not a product therefore we must have a modifier
                  if (Role == CFunctionParameter::PARAMETER)
                    {
                      Role = CFunctionParameter::MODIFIER;
                    }
                }
              else if (pObject->getObjectType() == "Model")
                {
                  Role = CFunctionParameter::TIME;
                }
              else if (pObject->getObjectType() == "Compartment")
                {
                  Role = CFunctionParameter::VOLUME;
                }

              const_cast< CFunctionParameter * >(Variables[i])->setUsage(Role);
            }

          pChild = dynamic_cast<const CEvaluationNodeObject*>(pChild->getSibling());
        }
    }
}

bool SBMLImporter::isSimpleFunctionCall(const CEvaluationNode* pRootNode)
{
  // it is a simple function call if it is a CEvaluationNodeCall object and all
  // its arguments are object nodes.
  bool result = true;

  if (dynamic_cast<const CEvaluationNodeCall*>(pRootNode))
    {
      const CEvaluationNode* pChildNode = static_cast<const CEvaluationNode*>(pRootNode->getChild());

      // I guess it must have at least one child to qualify.
      if (!pChildNode)
        {
          result = false;
        }

      while (pChildNode)
        {
          if (!dynamic_cast<const CEvaluationNodeObject*>(pChildNode))
            {
              result = false;
              break;
            }

          pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
        }
    }
  else
    {
      result = false;
    }

  return result;
}

ConverterASTNode* SBMLImporter::isMultipliedByVolume(const ASTNode* node, const std::string& compartmentSBMLId)
{
  ConverterASTNode* result = NULL;

  if (node->getType() == AST_TIMES || node->getType() == AST_DIVIDE)
    {
      ConverterASTNode* pTmpResultNode = new ConverterASTNode(node->getType());
      unsigned int i, iMax = node->getNumChildren();
      bool found = false;

      for (i = 0; i < iMax; ++i)
        {
          const ASTNode* child = node->getChild(i);

          if (node->getType() == AST_TIMES && child->getType() == AST_NAME && child->getName() == compartmentSBMLId)
            {
              found = true;
            }
          else if ((!found) && (child->getType() == AST_TIMES || child->getType() == AST_DIVIDE))
            {
              ASTNode* pSubResult = this->isMultipliedByVolume(child, compartmentSBMLId);

              if (pSubResult)
                {
                  found = true;

                  if (pSubResult->getNumChildren() > 1)
                    {
                      pTmpResultNode->addChild(pSubResult);
                    }
                  else if (pSubResult->getNumChildren() == 1)
                    {
                      pTmpResultNode->addChild(static_cast<ASTNode*>(static_cast<ConverterASTNode*>(pSubResult)->removeChild(0)));
                      delete pSubResult;
                    }
                  else
                    {
                      delete pSubResult;
                    }
                }
              else
                {
                  pTmpResultNode->addChild(new ConverterASTNode(*child));
                }
            }
          else
            {
              // bail in case of hOSU
              if (child->getType() == AST_NAME &&  mSubstanceOnlySpecies.size() > 0)
                {
                  std::map<Species*, Compartment*>::iterator it = this->mSubstanceOnlySpecies.begin();
                  std::map<Species*, Compartment*>::iterator endIt = this->mSubstanceOnlySpecies.end();

                  while (it != endIt)
                    {
                      if (it->first->getId() == child->getName()) return NULL;

                      it++;
                    }
                }

              pTmpResultNode->addChild(new ConverterASTNode(*child));
            }
        }

      if (found)
        {
          result = pTmpResultNode;
        }
      else
        {
          delete pTmpResultNode;
        }
    }

  return result;
}

CEvaluationNode* SBMLImporter::variables2objects(const CEvaluationNode* pOrigNode, const std::map<std::string, std::string>& replacementMap)
{
  // TODO CRITICAL We need to use a node iterator
  CEvaluationNode* pResultNode = NULL;

  if (dynamic_cast<const CEvaluationNodeVariable*>(pOrigNode))
    {
      std::map<std::string , std::string>::const_iterator pos = replacementMap.find(pOrigNode->getData());

      if (pos == replacementMap.end()) fatalError();

      pResultNode = new CEvaluationNodeObject(CEvaluationNodeObject::CN, "<" + pos->second + ">");
    }
  else
    {
      pResultNode = CEvaluationNode::create(pOrigNode->getType(), pOrigNode->getData());
      const CEvaluationNode* pChildNode = static_cast<const CEvaluationNode*>(pOrigNode->getChild());

      while (pChildNode)
        {
          pResultNode->addChild(this->variables2objects(pChildNode, replacementMap));
          pChildNode = static_cast<const CEvaluationNode*>(pChildNode->getSibling());
        }
    }

  return pResultNode;
}

void SBMLImporter::renameMassActionParameters(CEvaluationNodeCall* pCallNode)
{
  std::vector<CCopasiContainer*> v;
  v.push_back(this->mpCopasiModel);
  CEvaluationNodeObject* pObjectNode = dynamic_cast<CEvaluationNodeObject*>(pCallNode->getChild());
  assert(pObjectNode);
  CCopasiObjectName objectName = CCopasiObjectName(pObjectNode->getData().substr(1, pObjectNode->getData().length() - 2));
  CCopasiObject* pObject = mpDataModel->ObjectFromName(v, objectName);
  assert(pObject);

  if (dynamic_cast<CCopasiParameter*>(pObject))
    {
      pObject->setObjectName("k1");
      pObjectNode->setData("<" + pObject->getCN() + ">");
    }

  pObjectNode = dynamic_cast<CEvaluationNodeObject*>(pObjectNode->getSibling());

  if (pObjectNode)
    {
      objectName = CCopasiObjectName(pObjectNode->getData().substr(1, pObjectNode->getData().length() - 2));
      pObject = mpDataModel->ObjectFromName(v, objectName);
      assert(pObject);

      if (dynamic_cast<CCopasiParameter*>(pObject))
        {
          pObject->setObjectName("k2");
          pObjectNode->setData("<" + pObject->getCN() + ">");
        }
    }
}

bool SBMLImporter::containsVolume(const ASTNode* pNode, const std::string& compartmentSBMLId)
{
  bool result = false;
  unsigned int i, iMax = pNode->getNumChildren();

  for (i = 0; i < iMax; ++i)
    {
      if (pNode->getChild(i)->getType() == AST_NAME &&
          pNode->getChild(i)->getName() == compartmentSBMLId)
        {
          result = true;
          break;
        }
    }

  return result;
}

void SBMLImporter::setImportHandler(CProcessReport* pHandler)
{mpImportHandler = pHandler;}

CProcessReport* SBMLImporter::getImportHandlerAddr()
{return mpImportHandler;}

bool SBMLImporter::removeUnusedFunctions(CFunctionDB* pTmpFunctionDB, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap)
{
  if (pTmpFunctionDB)
    {
      unsigned C_INT32 step = 0, totalSteps = 0;
      size_t hStep = C_INVALID_INDEX;
      size_t i, iMax = this->mpCopasiModel->getReactions().size();

      if (mpImportHandler)
        {
          step = 0;
          totalSteps =
            (unsigned C_INT32)(iMax + this->mpCopasiModel->getCompartments().size() + this->mpCopasiModel->getMetabolites().size() + this->mpCopasiModel->getModelValues().size());
          hStep = mpImportHandler->addItem("Searching used functions...",
                                           step,
                                           &totalSteps);
        }

      std::set<std::string> functionNameSet;

      for (i = 0; i < iMax; ++i)
        {
          const CFunction* pTree = this->mpCopasiModel->getReactions()[i]->getFunction();

          // this is a hack, but if we changed stoichiometries in reactions, we can no longer be sure that the
          // kinetic function fits the reaction, so in this case it is probably better to set all functions used
          // in directions to general
          // If we don't do this, the user will not be able to use the function in the reaction in the GUI if he
          // changed it once.
          // TODO only set those functions to general that really contain modified stoichiometries
          if (this->mConversionFactorFound)
            {
              // TODO another const_cast that should disappear
              // TODO because there is probably a better place to do this
              // TODO but this means rewriting the code that modifies the
              // TODO stoichiometries of reactions and put it all in one place.
              const_cast<CFunction*>(pTree)->setReversible(TriUnspecified);
            }

          if (functionNameSet.find(pTree->getObjectName()) == functionNameSet.end())
            {
              functionNameSet.insert(pTree->getObjectName());
              this->findFunctionCalls(pTree->getRoot(), functionNameSet);
            }

          ++step;

          if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return false;
        }

      iMax = this->mpCopasiModel->getCompartments().size();

      for (i = 0; i < iMax; ++i)
        {
          CModelEntity* pME = this->mpCopasiModel->getCompartments()[i];

          if (pME->getStatus() != CModelEntity::FIXED)
            {
              const CEvaluationTree* pTree = pME->getExpressionPtr();

              if (pTree != NULL)
                {
                  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
                }
            }

          if (pME->getStatus() != CModelEntity::ASSIGNMENT)
            {
              const CEvaluationTree* pTree = pME->getInitialExpressionPtr();

              if (pTree != NULL)
                {
                  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
                }
            }

          ++step;
        }

      iMax = this->mpCopasiModel->getMetabolites().size();

      for (i = 0; i < iMax; ++i)
        {
          CModelEntity* pME = this->mpCopasiModel->getMetabolites()[i];

          if (pME->getStatus() != CModelEntity::FIXED)
            {
              const CEvaluationTree* pTree = pME->getExpressionPtr();

              if (pTree != NULL)
                {
                  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
                }
            }

          if (pME->getStatus() != CModelEntity::ASSIGNMENT)
            {
              const CEvaluationTree* pTree = pME->getInitialExpressionPtr();

              if (pTree != NULL)
                {
                  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
                }
            }

          ++step;
        }

      iMax = this->mpCopasiModel->getModelValues().size();

      for (i = 0; i < iMax; ++i)
        {
          CModelEntity* pME = this->mpCopasiModel->getModelValues()[i];

          if (pME->getStatus() != CModelEntity::FIXED)
            {
              const CEvaluationTree* pTree = pME->getExpressionPtr();

              if (pTree != NULL)
                {
                  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
                }
            }

          if (pME->getStatus() != CModelEntity::ASSIGNMENT)
            {
              const CEvaluationTree* pTree = pME->getInitialExpressionPtr();

              if (pTree != NULL)
                {
                  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
                }
            }

          ++step;
        }

      // find the used function in events
      iMax = this->mpCopasiModel->getEvents().size();

      for (i = 0; i < iMax; ++i)
        {
          CEvent* pEvent = this->mpCopasiModel->getEvents()[i];
          assert(pEvent != NULL);
          const CEvaluationTree* pTree = pEvent->getTriggerExpressionPtr();
          // an event has to have a trigger
          assert(pTree != NULL);

          if (pTree != NULL)
            {
              this->findFunctionCalls(pTree->getRoot(), functionNameSet);
            }

          // handle the delay
          pTree = pEvent->getDelayExpressionPtr();

          if (pTree != NULL)
            {
              this->findFunctionCalls(pTree->getRoot(), functionNameSet);
            }

          // handle all assignments
          size_t j, jMax = pEvent->getAssignments().size();

          for (j = 0; j < jMax; ++j)
            {
              CEventAssignment* pEventAssignment = pEvent->getAssignments()[j];
              assert(pEventAssignment != NULL);

              if (pEventAssignment != NULL)
                {
                  pTree = pEventAssignment->getExpressionPtr();
                  // each event assignment has to have an expression
                  assert(pTree != NULL);

                  if (pTree != NULL)
                    {
                      this->findFunctionCalls(pTree->getRoot(), functionNameSet);
                    }
                }
            }
        }

      CFunctionDB* pFunctionDB = CCopasiRootContainer::getFunctionList();

      if (mpImportHandler)
        {
          mpImportHandler->finishItem(hStep);
          step = 0;
          totalSteps = (unsigned C_INT32) pTmpFunctionDB->loadedFunctions().size();
          hStep = mpImportHandler->addItem("Removing unused functions...",
                                           step,
                                           &totalSteps);
        }

      // here we could have a dialog asking the user if unused functions should
      // be removed.

      CCopasiVectorN < CFunction >::const_iterator it = pTmpFunctionDB->loadedFunctions().begin();
      CCopasiVectorN < CFunction >::const_iterator end = pTmpFunctionDB->loadedFunctions().end();

      for (; it != end; ++it)
        {
          CEvaluationTree * pTree = *it;

          if (functionNameSet.find(pTree->getObjectName()) == functionNameSet.end())
            {
              mUsedFunctions.erase(pTree->getObjectName());
              pFunctionDB->loadedFunctions().remove(pTree->getObjectName());
              // delete the entry from the copasi2sbmlmap.
              std::map<CCopasiObject*, SBase*>::iterator pos = copasi2sbmlmap.find(pTree);
              assert(pos != copasi2sbmlmap.end());
              copasi2sbmlmap.erase(pos);
            }

          ++step;

          if (mpImportHandler && !mpImportHandler->progressItem(hStep)) return false;
        }

      if (mpImportHandler)
        {
          mpImportHandler->finishItem(hStep);
        }
    }

  return true;
}

void SBMLImporter::findFunctionCalls(const CEvaluationNode* pNode, std::set<std::string>& functionNameSet)
{
  if (pNode)
    {
      CFunctionDB* pFunctionDB = CCopasiRootContainer::getFunctionList();
      CCopasiTree<const CEvaluationNode>::iterator treeIt = pNode;

      while (treeIt != NULL)
        {
          if (CEvaluationNode::type((*treeIt).getType()) == CEvaluationNode::CALL)
            {
              // unQuote not necessary since getIndex in CCopasiVector takes care of this.
              CEvaluationTree* pTree = pFunctionDB->findFunction((*treeIt).getData());

              if (functionNameSet.find(pTree->getObjectName()) == functionNameSet.end())
                {
                  functionNameSet.insert(pTree->getObjectName());
                  this->findFunctionCalls(pTree->getRoot(), functionNameSet);
                }
            }

          ++treeIt;
        }
    }
}

bool SBMLImporter::isStochasticModel(const Model* pSBMLModel)
{
  bool stochastic = false;
  unsigned int i;
  const UnitDefinition* pUD = pSBMLModel->getUnitDefinition("substance");

  if (pUD)
    {
      stochastic = (pUD->getNumUnits() == 1 &&
                    pUD->getUnit(0)->getKind() == UNIT_KIND_ITEM);

      for (i = 0; (stochastic == true) && (i < pSBMLModel->getNumReactions()); ++i)
        {
          stochastic = !pSBMLModel->getReaction(i)->getReversible();
        }
    }

  return stochastic;
}

void SBMLImporter::importSBMLRule(const Rule* sbmlRule, std::map<CCopasiObject*, SBase*>& copasi2sbmlmap, Model* pSBMLModel)
{
  // so far we only support assignment rules and rate rules
  int type = sbmlRule->getTypeCode();

  if (type == SBML_ASSIGNMENT_RULE)
    {
      const AssignmentRule* pAssignmentRule = dynamic_cast<const AssignmentRule*>(sbmlRule);

      if (pAssignmentRule && pAssignmentRule->isSetVariable())
        {
          this->importRule(pAssignmentRule, CModelEntity::ASSIGNMENT, copasi2sbmlmap, pSBMLModel);
        }
      else
        {
          fatalError();
        }
    }
  else if (type == SBML_RATE_RULE)
    {
      const RateRule* pRateRule = dynamic_cast<const RateRule*