<commit_before>/* Begin CVS Header
   $Source: /Volumes/Home/Users/shoops/cvs/copasi_dev/copasi/trajectory/CStochMethod.cpp,v $
   $Revision: 1.53 $
   $Name:  $
   $Author: ssahle $ 
   $Date: 2006/04/12 14:35:01 $
   End CVS Header */

#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

#include <mathematics.h>
#include <vector>
#include <numeric>
#include <limits.h>
#include <set>

#include "copasi.h"
#include "CStochMethod.h"
#include "CStochDirectMethod.h"
#include "CStochNextReactionMethod.h"
#include "utilities/CCopasiVector.h"
#include "function/CFunction.h"
#include "randomGenerator/CRandom.h"
#include "CTrajectoryMethod.h"
#include "CTrajectoryProblem.h"
#include "model/CState.h"
#include "model/CCompartment.h"
#include "model/CModel.h"

C_INT32 CStochMethod::checkModel(CModel * C_UNUSED(pmodel))
{
  // Here several checks will be performed to validate the model
  return 2; // suggest next reaction method

  // TODO check if stoich is integer
}

CStochMethod *
CStochMethod::createStochMethod(CTrajectoryProblem * pProblem)
{
  C_INT32 result = 2; // next reaction method as default
  if (pProblem && pProblem->getModel())
    {
      result = checkModel(pProblem->getModel());
    }

  CStochMethod * method = NULL;

  if (result < 0)
    {
      //error. stochastic simulation not possible
    }
  else if (result == 1)
    {
      method = new CStochDirectMethod();
    }
  else if (result == 2)
    {
      method = new CStochNextReactionMethod();
    }
  else
    {
      //should not happen
    }

  return method;
}

CStochMethod::CStochMethod(const CCopasiContainer * pParent):
    CTrajectoryMethod(CCopasiMethod::stochastic, pParent)
{
  // Max number of doSingleStep() per step()
  addParameter("STOCH.MaxSteps",
               CCopasiParameter::INT, (C_INT32) 1000000);
  // 0: default; 1: direct method; 2: next reaction method
  addParameter("STOCH.Subtype",
               CCopasiParameter::UINT, (unsigned C_INT32) 2);

  addParameter("STOCH.UseRandomSeed",
               CCopasiParameter::BOOL, false);
  addParameter("STOCH.RandomSeed",
               CCopasiParameter::UINT, (unsigned C_INT32) 1);

  mpRandomGenerator = CRandom::createGenerator(CRandom::r250);
}

CStochMethod::~CStochMethod()
{
  delete mpRandomGenerator;
  mpRandomGenerator = NULL;
}

void CStochMethod::step(const double & deltaT)
{
  // write the current state to the model:
  //mpProblem->getModel()->setState(mpCurrentState); //?

  // check for possible overflows:
  unsigned C_INT32 i;
  unsigned C_INT32 imax;

  for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++)
    if (mpProblem->getModel()->getMetabolitesX()[i]->getValue() >= mMaxIntBeforeStep)
      {
        CCopasiMessage(CCopasiMessage::EXCEPTION, "at least one particle number got to big.");
        // TODO:throw exception or something like that
      }

  // do several steps:
  C_FLOAT64 time = mpCurrentState->getTime();
  C_FLOAT64 endtime = time + deltaT;

  for (i = 0; ((i < (unsigned C_INT32) mMaxSteps) && (time < endtime)); i++)
    {
      time = doSingleStep(time, endtime);
    }
  mpCurrentState->setTime(time);

  if ((i >= (unsigned C_INT32) mMaxSteps) && (!mMaxStepsReached))
    {
      mMaxStepsReached = true; //only report this message once
      CCopasiMessage(CCopasiMessage::WARNING, "maximum number of reaction events was reached in at least one simulation step.\nThat means time intervals in the output may not be what you requested.");
    }

  // get back the particle numbers:

  /* Set the variable Metabolites */
  C_FLOAT64 * Dbl = mpCurrentState->beginIndependent();
  for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++, Dbl++)
    *Dbl = mpProblem->getModel()->getMetabolitesX()[i]->getValue();

  return;
}

void CStochMethod::start(const CState * initialState)
{
  /* get configuration data */
  mMaxSteps = * getValue("STOCH.MaxSteps").pINT;

  bool useRandomSeed = * getValue("STOCH.UseRandomSeed").pBOOL;
  unsigned C_INT32 randomSeed = * getValue("STOCH.RandomSeed").pUINT;
  if (useRandomSeed) mpRandomGenerator->initialize(randomSeed);

  *mpCurrentState = *initialState; //TODO seem to be identical

  unsigned C_INT32 i, imax;

  mNumNumbers = mpCurrentState->getNumVariable();
  mNumbers.resize(mNumNumbers);
  C_FLOAT64 * Dbl = mpCurrentState->beginIndependent();
  for (i = 0; i < mNumNumbers; ++i, Dbl++)
    {
      mNumbers[i] = (C_INT64) * Dbl;
      *Dbl = floor(*Dbl);
    }
  //TODO also put fixed variables here

  imax = mpCurrentState->getNumFixed();
  for (i = 0; i < imax; ++i, Dbl++)
    *Dbl = floor(*Dbl);

  mpModel = mpProblem->getModel();
  mpModel->setState(*mpCurrentState);

  mNumReactions = mpModel->getReactions().size();

  mAmu.clear(); mAmuOld.clear();
  for (i = 0; i < mNumReactions; i++)
    {
      mAmu.push_back(0);
      mAmuOld.push_back(0);
    }

  setupDependencyGraphAndBalances();
  //std::cout << mDG;
  updatePropensities();

  // call init of the specific simulation method
  initMethod(mpCurrentState->getTime());

  mMaxStepsReached = false;
  return;
}

C_INT32 CStochMethod::updatePropensities()<commit_msg>fix #490 now only metabolite particle numbers are rounded to integer values (not volumes, ...)<commit_after>/* Begin CVS Header
   $Source: /Volumes/Home/Users/shoops/cvs/copasi_dev/copasi/trajectory/CStochMethod.cpp,v $
   $Revision: 1.54 $
   $Name:  $
   $Author: ssahle $ 
   $Date: 2006/04/19 07:55:11 $
   End CVS Header */

#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

#include <mathematics.h>
#include <vector>
#include <numeric>
#include <limits.h>
#include <set>

#include "copasi.h"
#include "CStochMethod.h"
#include "CStochDirectMethod.h"
#include "CStochNextReactionMethod.h"
#include "utilities/CCopasiVector.h"
#include "function/CFunction.h"
#include "randomGenerator/CRandom.h"
#include "CTrajectoryMethod.h"
#include "CTrajectoryProblem.h"
#include "model/CState.h"
#include "model/CCompartment.h"
#include "model/CModel.h"

C_INT32 CStochMethod::checkModel(CModel * C_UNUSED(pmodel))
{
  // Here several checks will be performed to validate the model
  return 2; // suggest next reaction method

  // TODO check if stoich is integer
}

CStochMethod *
CStochMethod::createStochMethod(CTrajectoryProblem * pProblem)
{
  C_INT32 result = 2; // next reaction method as default
  if (pProblem && pProblem->getModel())
    {
      result = checkModel(pProblem->getModel());
    }

  CStochMethod * method = NULL;

  if (result < 0)
    {
      //error. stochastic simulation not possible
    }
  else if (result == 1)
    {
      method = new CStochDirectMethod();
    }
  else if (result == 2)
    {
      method = new CStochNextReactionMethod();
    }
  else
    {
      //should not happen
    }

  return method;
}

CStochMethod::CStochMethod(const CCopasiContainer * pParent):
    CTrajectoryMethod(CCopasiMethod::stochastic, pParent)
{
  // Max number of doSingleStep() per step()
  addParameter("STOCH.MaxSteps",
               CCopasiParameter::INT, (C_INT32) 1000000);
  // 0: default; 1: direct method; 2: next reaction method
  addParameter("STOCH.Subtype",
               CCopasiParameter::UINT, (unsigned C_INT32) 2);

  addParameter("STOCH.UseRandomSeed",
               CCopasiParameter::BOOL, false);
  addParameter("STOCH.RandomSeed",
               CCopasiParameter::UINT, (unsigned C_INT32) 1);

  mpRandomGenerator = CRandom::createGenerator(CRandom::r250);
}

CStochMethod::~CStochMethod()
{
  delete mpRandomGenerator;
  mpRandomGenerator = NULL;
}

void CStochMethod::step(const double & deltaT)
{
  // write the current state to the model:
  //mpProblem->getModel()->setState(mpCurrentState); //?

  // check for possible overflows:
  unsigned C_INT32 i;
  unsigned C_INT32 imax;

  for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++)
    if (mpProblem->getModel()->getMetabolitesX()[i]->getValue() >= mMaxIntBeforeStep)
      {
        CCopasiMessage(CCopasiMessage::EXCEPTION, "at least one particle number got to big.");
        // TODO:throw exception or something like that
      }

  // do several steps:
  C_FLOAT64 time = mpCurrentState->getTime();
  C_FLOAT64 endtime = time + deltaT;

  for (i = 0; ((i < (unsigned C_INT32) mMaxSteps) && (time < endtime)); i++)
    {
      time = doSingleStep(time, endtime);
    }
  mpCurrentState->setTime(time);

  if ((i >= (unsigned C_INT32) mMaxSteps) && (!mMaxStepsReached))
    {
      mMaxStepsReached = true; //only report this message once
      CCopasiMessage(CCopasiMessage::WARNING, "maximum number of reaction events was reached in at least one simulation step.\nThat means time intervals in the output may not be what you requested.");
    }

  // get back the particle numbers:

  /* Set the variable Metabolites */
  C_FLOAT64 * Dbl = mpCurrentState->beginIndependent();
  for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++, Dbl++)
    *Dbl = mpProblem->getModel()->getMetabolitesX()[i]->getValue();

  return;
}

void CStochMethod::start(const CState * initialState)
{
  /* get configuration data */
  mMaxSteps = * getValue("STOCH.MaxSteps").pINT;

  bool useRandomSeed = * getValue("STOCH.UseRandomSeed").pBOOL;
  unsigned C_INT32 randomSeed = * getValue("STOCH.RandomSeed").pUINT;
  if (useRandomSeed) mpRandomGenerator->initialize(randomSeed);

  *mpCurrentState = *initialState; //TODO seem to be identical

  mpModel = mpProblem->getModel();

  unsigned C_INT32 i, imax;

  mNumNumbers = mpCurrentState->getNumVariable();
  mNumbers.resize(mNumNumbers);
  C_FLOAT64 * Dbl = mpCurrentState->beginIndependent();
  CModelEntity*const* obj = mpModel->getStateTemplate().beginIndependent();
  for (i = 0; i < mNumNumbers; ++i, Dbl++, ++obj)
    {
      mNumbers[i] = (C_INT64) * Dbl;
      *Dbl = floor(*Dbl);
      //std::cout << (*obj)->getObjectName() << std::endl;
      //obj can later be used to handle variables differently
    }

  imax = mpCurrentState->getNumFixed();
  for (i = 0; i < imax; ++i, Dbl++, ++obj)
    {
      //std::cout << (*obj)->getObjectName() << std::endl;
      if (dynamic_cast<const CMetab*>(*obj))
        *Dbl = floor(*Dbl);
    }

  mpModel->setState(*mpCurrentState);

  mNumReactions = mpModel->getReactions().size();

  mAmu.clear(); mAmuOld.clear();
  for (i = 0; i < mNumReactions; i++)
    {
      mAmu.push_back(0);
      mAmuOld.push_back(0);
    }

  setupDependencyGraphAndBalances();
  //std::cout << mDG;
  updatePropensities();

  // call init of the specific simulation method
  initMethod(mpCurrentState->getTime());

  mMaxStepsReached = false;
  return;
}

C_INT32 CStochMethod::updatePropensities()
{'fraction_non_alphanumeric': 0.1351845085860431, 'fraction_numerical': 0.01835951772013153, 'mean_word_length': 3.7542336083369516, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 71, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2957085', 'n_tokens_mistral': 3900, 'n_tokens_neox': 3668, 'n_words': 878}