COPASI API  4.16.103
CTimeSeries.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2013 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2004 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include <limits>
16 
17 #include "copasi.h"
18 
19 #include "CTimeSeries.h"
20 
24 #include "model/CModel.h"
25 #include "report/CKeyFactory.h"
27 
28 #include "sbml/SBase.h"
29 #include "sbml/Compartment.h"
30 #include "sbml/Species.h"
31 #include "sbml/Parameter.h"
32 #include "sbml/Model.h"
33 
34 // static
35 std::string CTimeSeries::mDummyString("");
36 
37 // static
39 
42  CMatrix< C_FLOAT64 >(),
43  mAllocatedSteps(0),
44  mRecordedSteps(0),
45  mpIt(mArray),
46  mpEnd(mArray + size()),
47  mpState(NULL),
48  mTitles(),
49  mCompartment(),
50  mPivot(),
51  mKeys(),
52  mNumberToQuantityFactor(0.0)
53 {}
54 
56  COutputInterface(src),
57  CMatrix< C_FLOAT64 >(src),
58  mAllocatedSteps(src.mAllocatedSteps),
59  mRecordedSteps(src.mRecordedSteps),
60  mpIt(mArray + mRecordedSteps * mCols),
61  mpEnd(mArray + size()),
62  mpState(src.mpState),
63  mTitles(src.mTitles),
64  mCompartment(src.mCompartment),
65  mPivot(src.mPivot),
66  mKeys(src.mKeys),
67  mNumberToQuantityFactor(src.mNumberToQuantityFactor)
68 {}
69 
71 {}
72 
73 void CTimeSeries::allocate(const size_t & steps)
74 {
75  // The actual allocation is deferred to compile
76  mAllocatedSteps = steps;
77 }
78 
80 {
81  size_t diff;
82  diff = mAllocatedSteps / 4;
83 
84  if (diff < 10)
85  diff = 10;
86  else if (diff > 10000)
87  diff = 10000;
88 
89  mAllocatedSteps += diff;
91 
93  mpEnd = mArray + size();
94 }
95 
97 {
98  mObjects.clear();
101  mRecordedSteps = 0;
102  mpIt = mArray;
103  mpEnd = mArray + size();
104  mpState = NULL;
105  mTitles.clear();
106  mCompartment.resize(0);
107  mPivot.resize(0);
108  mKeys.clear();
110 }
111 
112 // virtual
113 bool CTimeSeries::compile(std::vector< CCopasiContainer * > listOfContainer,
114  const CCopasiDataModel* pDataModel)
115 {
116  const CModel * pModel =
117  dynamic_cast< const CModel * >(pDataModel->ObjectFromName(listOfContainer, pDataModel->getModel()->getCN()));
118 
119  if (pModel == NULL)
120  return false;
121 
122  mpState = & pModel->getState();
123 
124  const CStateTemplate & Template = pModel->getStateTemplate();
125 
126  // We store all variables of the system.
127  // The reason for this is that events will be able to change even fixed values.
128  CModelEntity *const* it = Template.getEntities();
129  CModelEntity *const* end = Template.endFixed();
130 
131  size_t i, imax = end - it;
132 
134 
135  mObjects.clear();
136 
137  mPivot.resize(imax);
138  mTitles.resize(imax);
139  mCompartment.resize(imax);
140  mKeys.resize(imax);
141 
142  mRecordedSteps = 0;
143  mpIt = mArray;
144  mpEnd = mArray + size();
146 
148 
149  const CMetab * pMetab;
150 
151  for (i = 0; it != end; ++i, ++it)
152  {
153  if ((pMetab = dynamic_cast< const CMetab *>(*it)) != NULL)
154  {
155  mTitles[i] = CMetabNameInterface::getDisplayName(pModel, *pMetab, false);
156  mCompartment[i] = Template.getIndex(pMetab->getCompartment());
157  }
158  else
159  {
160  mTitles[i] = (*it)->getObjectDisplayName();
161  }
162 
163  mKeys[i] = (*it)->getKey();
164 
165  mObjects.insert((*it)->getValueReference());
166  }
167 
168  mTitles[0] = "Time";
169  mKeys[0] = pModel->getKey();
170 
171  const size_t * pUserOrder = Template.getUserOrder().array();
172  const size_t * pUserOrderEnd = pUserOrder + Template.getUserOrder().size();
173  it = Template.getEntities();
174 
175  for (i = 0; pUserOrder != pUserOrderEnd; ++pUserOrder)
176  mPivot[i++] = *pUserOrder;
177 
178  return true;
179 }
180 
181 // virtual
183 {
184  if (activity != DURING)
185  return;
186 
187  // We may have to reallocate due to additional output caused from events
188  if (mpIt == mpEnd)
189  {
191  }
192 
193  if (mpIt != mpEnd)
194  {
195  memcpy(mpIt, &mpState->getTime(), mCols * sizeof(C_FLOAT64));
196  mpIt += mCols;
197  mRecordedSteps++;
198  }
199 }
200 
201 // virtual
203 {
204  // We may have to reallocate due to additional output caused from events
205  if (mpIt == mpEnd)
206  {
208  }
209 
210  if (mpIt != mpEnd)
211  {
212  C_FLOAT64 * pIt = mpIt;
213  mpIt += mCols;
214  mRecordedSteps++;
215 
216  // We copy NaN to indicate separation, which is similar to plotting.
217  for (; pIt != mpIt; ++pIt)
218  *pIt = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
219  }
220 }
221 
222 // virtual
224 {}
225 
226 //*** the methods to retrieve data from the CTimeSeries *******
227 
228 const size_t & CTimeSeries::getRecordedSteps() const
229 {return mRecordedSteps;}
230 
231 const size_t & CTimeSeries::getNumVariables() const
232 {return mCols;}
233 
234 const C_FLOAT64 & CTimeSeries::getData(const size_t & step,
235  const size_t & var) const
236 {
237  if (step < mRecordedSteps && var < mCols)
238  return *(mArray + step * mCols + mPivot[var]);
239 
240  return mDummyFloat;
241 }
242 
244  const size_t & var) const
245 {
246  if (step < mRecordedSteps && var < mCols)
247  {
248  const size_t & Col = mPivot[var];
249 
250  if (mCompartment[Col] != C_INVALID_INDEX)
251  return *(mArray + step * mCols + Col) * mNumberToQuantityFactor / *(mArray + step * mCols + mCompartment[Col]);
252  else
253  return *(mArray + step * mCols + Col);
254  }
255 
256  return mDummyFloat;
257 }
258 
259 const std::string & CTimeSeries::getTitle(const size_t & var) const
260 {
261  if (var < mCols)
262  return mTitles[mPivot[var]];
263 
264  return mDummyString;
265 }
266 
267 const std::string & CTimeSeries::getKey(const size_t & var) const
268 {
269  if (var < mCols)
270  return mKeys[mPivot[var]];
271 
272  return mDummyString;
273 }
274 
275 std::string CTimeSeries::getSBMLId(const size_t & var, const CCopasiDataModel* pDataModel) const
276 {
277  std::string key = getKey(var);
278  std::string result("");
279 
280  if (key != mDummyString)
281  {
282  const CCopasiObject* pObject = CCopasiRootContainer::getKeyFactory()->get(key);
283 
284  if (pObject != NULL)
285  {
286  std::map<CCopasiObject*, SBase*>::const_iterator pos = const_cast<CCopasiDataModel*>(pDataModel)->getCopasi2SBMLMap().find(const_cast<CCopasiObject*>(pObject));
287 
288  if (pos != const_cast<CCopasiDataModel*>(pDataModel)->getCopasi2SBMLMap().end())
289  {
290  const SBase* pSBMLObject = pos->second;
291  const Compartment* pSBMLCompartment = NULL;
292  const Species* pSBMLSpecies = NULL;
293  const Parameter* pSBMLParameter = NULL;
294  const Model* pSBMLModel = NULL;
295 
296  switch (pSBMLObject->getTypeCode())
297  {
298  case SBML_COMPARTMENT:
299  pSBMLCompartment = dynamic_cast<const Compartment*>(pSBMLObject);
300 
301  if (pSBMLCompartment && pSBMLCompartment->isSetId())
302  {
303  result = pSBMLCompartment->getId();
304  }
305 
306  break;
307 
308  case SBML_SPECIES:
309  pSBMLSpecies = dynamic_cast<const Species*>(pSBMLObject);
310 
311  if (pSBMLSpecies && pSBMLSpecies->isSetId())
312  {
313  result = pSBMLSpecies->getId();
314  }
315 
316  break;
317 
318  case SBML_PARAMETER:
319  pSBMLParameter = dynamic_cast<const Parameter*>(pSBMLObject);
320 
321  if (pSBMLParameter && pSBMLParameter->isSetId())
322  {
323  result = pSBMLParameter->getId();
324  }
325 
326  break;
327 
328  case SBML_MODEL:
329  pSBMLModel = dynamic_cast<const Model*>(pSBMLObject);
330 
331  if (pSBMLModel && pSBMLModel->isSetId())
332  {
333  result = pSBMLModel->getId();
334  }
335 
336  break;
337 
338  default:
339  break;
340  }
341  }
342  }
343  }
344 
345  return result;
346 }
347 
348 int CTimeSeries::save(const std::string& fileName, bool writeParticleNumbers, const std::string& separator) const
349 {
350  std::ofstream fileStream(CLocaleString::fromUtf8(fileName).c_str());
351  std::ostringstream* stringStream = new std::ostringstream();
352  (*stringStream) << "# ";
353  size_t counter2;
354  size_t maxCount2 = this->getNumVariables();
355 
356  for (counter2 = 0; counter2 < maxCount2; ++counter2)
357  {
358  (*stringStream) << this->getTitle(counter2) << separator;
359  }
360 
361  (*stringStream) << std::endl;
362  fileStream << stringStream->str();
363 
364  if (!fileStream.good()) return 1;
365 
366  size_t counter;
367  size_t maxCount = mRecordedSteps;
368 
369  for (counter = 0; counter < maxCount; ++counter)
370  {
371  delete stringStream;
372  stringStream = new std::ostringstream();
373 
374  for (counter2 = 0; counter2 < maxCount2; ++counter2)
375  {
376  C_FLOAT64 value;
377 
378  if (writeParticleNumbers)
379  {
380  value = this->getData(counter, counter2);
381  }
382  else
383  {
384  value = this->getConcentrationData(counter, counter2);
385  }
386 
387  (*stringStream) << value << separator;
388  }
389 
390  (*stringStream) << std::endl;
391  fileStream << stringStream->str();
392 
393  if (!fileStream.good()) return 1;
394  }
395 
396  fileStream.close();
397  delete stringStream;
398  return 0;
399 }
const size_t & getRecordedSteps() const
C_FLOAT64 getConcentrationData(const size_t &step, const size_t &variable) const
std::string getSBMLId(const size_t &variable, const CCopasiDataModel *pDataModel) const
virtual bool compile(std::vector< CCopasiContainer * > listOfContainer, const CCopasiDataModel *pDataModel)
const size_t & getNumVariables() const
C_FLOAT64 * mpEnd
Definition: CTimeSeries.h:182
virtual CCopasiObjectName getCN() const
int save(const std::string &fileName, bool writeConcentrations=false, const std::string &separator="\t") const
void clear()
Definition: CTimeSeries.cpp:96
std::vector< std::string > mTitles
Definition: CTimeSeries.h:192
CCopasiObject * get(const std::string &key)
size_t mRecordedSteps
Definition: CTimeSeries.h:172
const std::string & getTitle(const size_t &variable) const
void allocate(const size_t &steps)
Definition: CTimeSeries.cpp:73
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
#define C_INVALID_INDEX
Definition: copasi.h:222
size_t getIndex(const CModelEntity *entity) const
Definition: CState.cpp:231
C_FLOAT64 * mArray
Definition: CMatrix.h:87
static std::string getDisplayName(const CModel *model, const std::string &key, const bool &quoted)
Definition: CMetab.h:178
const CState * mpState
Definition: CTimeSeries.h:187
const std::string & getKey(const size_t &variable) const
static std::string mDummyString
Definition: CTimeSeries.h:218
virtual void finish()
const C_FLOAT64 & getNumber2QuantityFactor() const
Definition: CModel.cpp:2357
CModelEntity ** endFixed()
Definition: CState.cpp:213
CVector< size_t > mPivot
Definition: CTimeSeries.h:203
static C_FLOAT64 mDummyFloat
Definition: CTimeSeries.h:223
const std::string & getKey() const
Definition: CModel.cpp:1142
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
std::set< const CCopasiObject * > mObjects
size_t size() const
Definition: CVector.h:100
const C_FLOAT64 & getData(const size_t &step, const size_t &variable) const
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
static CKeyFactory * getKeyFactory()
#define C_FLOAT64
Definition: copasi.h:92
C_FLOAT64 mNumberToQuantityFactor
Definition: CTimeSeries.h:213
CType * array()
Definition: CVector.h:139
const CStateTemplate & getStateTemplate() const
Definition: CModel.cpp:1172
virtual size_t size() const
Definition: CMatrix.h:132
C_FLOAT64 * mpIt
Definition: CTimeSeries.h:177
Definition: CModel.h:50
const CState & getState() const
Definition: CModel.cpp:1771
static CLocaleString fromUtf8(const std::string &utf8)
const CVector< size_t > & getUserOrder() const
Definition: CState.cpp:201
void increaseAllocation()
Definition: CTimeSeries.cpp:79
size_t mAllocatedSteps
Definition: CTimeSeries.h:167
virtual void separate(const Activity &activity)
CCopasiObject * ObjectFromName(const std::vector< CCopasiContainer * > &listOfContainer, const CCopasiObjectName &CN) const
virtual void output(const Activity &activity)
const CCompartment * getCompartment() const
Definition: CMetab.cpp:222
CVector< size_t > mCompartment
Definition: CTimeSeries.h:198
CModelEntity ** getEntities()
Definition: CState.cpp:204
std::vector< std::string > mKeys
Definition: CTimeSeries.h:208