COPASI API  4.16.103
CExperiment.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2015 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2005 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include <fstream>
16 #include <limits>
17 #include <cmath>
18 
19 #include "copasi.h"
20 
21 #include "CExperiment.h"
22 #include "CExperimentObjectMap.h"
23 #include "CFitTask.h"
24 
27 #include "model/CModel.h"
29 #include "report/CKeyFactory.h"
30 #include "utilities/CTableCell.h"
31 #include "utilities/CSort.h"
32 #include "utilities/CDirEntry.h"
33 #include "utilities/utility.h"
35 
36 std::istream & skipLine(std::istream & in);
37 
38 #define InvalidIndex std::numeric_limits< unsigned C_INT32 >::max()
39 
40 const std::string CExperiment::TypeName[] =
41 {
42  "ignored",
43  "independent",
44  "dependent",
45  "Time",
46  ""
47 };
48 
49 const char* CExperiment::XMLType[] =
50 {
51  "ignored",
52  "independent",
53  "dependent",
54  "time",
55  NULL
56 };
57 
58 const std::string CExperiment::WeightMethodName[] =
59 {
60  "Mean",
61  "Mean Square",
62  "Standard Deviation",
63  "Value Scaling",
64  ""
65 };
66 
67 const char* CExperiment::WeightMethodType[] =
68 {
69  "Mean",
70  "MeanSquare",
71  "StandardDeviation",
72  "ValueScaling",
73  NULL
74 };
75 
77  const std::string & name):
78  CCopasiParameterGroup(name, pParent),
79  mpFileName(NULL),
80  mpFirstRow(NULL),
81  mpLastRow(NULL),
82  mpTaskType(NULL),
83  mpNormalizeWeightsPerExperiment(NULL),
84  mpSeparator(NULL),
85  mpWeightMethod(NULL),
86  mpRowOriented(NULL),
87  mpHeaderRow(NULL),
88  mpNumColumns(NULL),
89  mColumnName(),
90  mpObjectMap(NULL),
91  mDataTime(0),
92  mDataIndependent(0, 0),
93  mDataDependent(0, 0),
94  mScale(0, 0),
95  mMissingData(false),
96  mMeans(0),
97  mColumnScale(0),
98  mDefaultColumnScale(0),
99  mDependentValues(0),
100  mIndependentUpdateMethods(0),
101  mRefreshMethods(),
102  mIndependentObjects(),
103  mIndependentValues(0),
104  mNumDataRows(0),
105  mpDataDependentCalculated(NULL),
106  mMean(0),
107  mMeanSD(0),
108  mObjectiveValue(0),
109  mRMS(0),
110  mRowObjectiveValue(0),
111  mRowRMS(0),
112  mColumnObjectiveValue(0),
113  mColumnRMS(0),
114  mColumnValidValueCount(0),
115  mDependentObjects(),
116  mFittingPoints("Fitted Points", this),
117  mExtendedTimeSeries(0),
118  mStorageIt(NULL),
119  mExtendedTimeSeriesSize(0)
120 {
122 
124 }
125 
127  const CCopasiContainer * pParent):
128  CCopasiParameterGroup(src, static_cast< const CCopasiContainer * >((pParent != NULL) ? pParent : src.getObjectDataModel())),
129  mpFileName(NULL),
130  mpFirstRow(NULL),
131  mpLastRow(NULL),
132  mpTaskType(NULL),
133  mpNormalizeWeightsPerExperiment(NULL),
134  mpSeparator(NULL),
135  mpWeightMethod(NULL),
136  mpRowOriented(NULL),
137  mpHeaderRow(NULL),
138  mpNumColumns(NULL),
139  mColumnName(src.mColumnName),
140  mpObjectMap(NULL),
141  mDataTime(src.mDataTime),
142  mDataIndependent(src.mDataIndependent),
143  mDataDependent(src.mDataDependent),
144  mScale(src.mScale),
145  mMissingData(src.mMissingData),
146  mMeans(src.mMeans),
147  mColumnScale(src.mColumnScale),
148  mDefaultColumnScale(src.mDefaultColumnScale),
149  mDependentValues(src.mDependentValues),
150  mIndependentUpdateMethods(src.mIndependentUpdateMethods),
151  mRefreshMethods(src.mRefreshMethods),
152  mIndependentObjects(src.mIndependentObjects),
153  mIndependentValues(src.mIndependentValues),
154  mNumDataRows(src.mNumDataRows),
155  mpDataDependentCalculated(src.mpDataDependentCalculated),
156  mMean(src.mMean),
157  mMeanSD(src.mMeanSD),
158  mObjectiveValue(0),
159  mRMS(src.mRMS),
160  mRowObjectiveValue(src.mRowObjectiveValue),
161  mRowRMS(src.mRowRMS),
162  mColumnObjectiveValue(src.mColumnObjectiveValue),
163  mColumnRMS(src.mColumnRMS),
164  mColumnValidValueCount(src.mColumnValidValueCount),
165  mDependentObjects(src.mDependentObjects),
166  mFittingPoints(src.mFittingPoints, this),
167  mExtendedTimeSeries(src.mExtendedTimeSeries),
168  mStorageIt(src.mStorageIt),
169  mExtendedTimeSeriesSize(src.mExtendedTimeSeriesSize)
170 
171 {
173 
175 }
176 
178  const CCopasiContainer * pParent):
179  CCopasiParameterGroup(group, static_cast< const CCopasiContainer * >((pParent != NULL) ? pParent : group.getObjectDataModel())),
180  mpFileName(NULL),
181  mpFirstRow(NULL),
182  mpLastRow(NULL),
183  mpTaskType(NULL),
184  mpNormalizeWeightsPerExperiment(NULL),
185  mpSeparator(NULL),
186  mpWeightMethod(NULL),
187  mpRowOriented(NULL),
188  mpHeaderRow(NULL),
189  mpNumColumns(NULL),
190  mColumnName(),
191  mpObjectMap(NULL),
192  mDataTime(0),
193  mDataIndependent(0, 0),
194  mDataDependent(0, 0),
195  mScale(0, 0),
196  mMissingData(false),
197  mMeans(0),
198  mColumnScale(0),
199  mDefaultColumnScale(0),
200  mDependentValues(0),
201  mIndependentUpdateMethods(0),
202  mRefreshMethods(),
203  mIndependentObjects(),
204  mIndependentValues(0),
205  mNumDataRows(0),
206  mpDataDependentCalculated(NULL),
207  mMean(0),
208  mMeanSD(0),
209  mObjectiveValue(0),
210  mRMS(0),
211  mRowObjectiveValue(0),
212  mRowRMS(0),
213  mColumnObjectiveValue(0),
214  mColumnRMS(0),
215  mColumnValidValueCount(0),
216  mDependentObjects(),
217  mFittingPoints("Fitted Points", this),
218  mExtendedTimeSeries(0),
219  mStorageIt(NULL),
220  mExtendedTimeSeriesSize(0)
221 {
223 
225 }
226 
228 
230 {
231  std::string Key = *getValue("Key").pKEY;
232 
233  clear();
234 
235  *static_cast<CCopasiParameterGroup *>(this) =
236  *static_cast<const CCopasiParameterGroup *>(&rhs);
237 
238  setValue("Key", Key);
239 
240  mpFileName = getValue("File Name").pFILE;
241  mpFirstRow = getValue("First Row").pUINT;
242  mpLastRow = getValue("Last Row").pUINT;
243  mpTaskType = (CCopasiTask::Type *) getValue("Experiment Type").pUINT;
244  mpNormalizeWeightsPerExperiment = getValue("Normalize Weights per Experiment").pBOOL;
245  mpSeparator = getValue("Separator").pSTRING;
246  mpWeightMethod = (WeightMethod *) getValue("Weight Method").pUINT;
247  mpRowOriented = getValue("Data is Row Oriented").pBOOL;
248  mpHeaderRow = getValue("Row containing Names").pUINT;
249  mpNumColumns = getValue("Number of Columns").pUINT;
250 
251  elevateChildren();
252 
253  return *this;
254 }
255 
257 {
259  mKey = CCopasiRootContainer::getKeyFactory()->add("Experiment", this);
260 
262 
263  mpFileName =
264  assertParameter("File Name", CCopasiParameter::FILE, std::string(""))->getValue().pFILE;
265  mpFirstRow =
267  mpLastRow =
268  assertParameter("Last Row", CCopasiParameter::UINT, (unsigned C_INT32) InvalidIndex)->getValue().pUINT;
270  assertParameter("Experiment Type", CCopasiParameter::UINT, (unsigned C_INT32) CCopasiTask::unset)->getValue().pUINT;
271 
273  assertParameter("Normalize Weights per Experiment", CCopasiParameter::BOOL, true)->getValue().pBOOL;
274 
275  mpSeparator =
276  assertParameter("Separator", CCopasiParameter::STRING, std::string("\t"))->getValue().pSTRING;
278  assertParameter("Weight Method", CCopasiParameter::UINT, (unsigned C_INT32) MEAN_SQUARE)->getValue().pUINT;
279  mpRowOriented =
280  assertParameter("Data is Row Oriented", CCopasiParameter::BOOL, (bool) true)->getValue().pBOOL;
281  mpHeaderRow =
282  assertParameter("Row containing Names", CCopasiParameter::UINT, (unsigned C_INT32) InvalidIndex)->getValue().pUINT;
283  mpNumColumns =
284  assertParameter("Number of Columns", CCopasiParameter::UINT, (unsigned C_INT32) 0)->getValue().pUINT;
285 
286  assertGroup("Object Map");
287 
288  // Old files have a typo in the name of the separator parameter.
289  CCopasiParameter * pParameter = getParameter("Seperator");
290 
291  if (pParameter != NULL)
292  {
293  *mpSeparator = *pParameter->getValue().pSTRING;
294  removeParameter("Seperator");
295  }
296 
297  elevateChildren();
298 }
299 
301 {
302  mpObjectMap =
303  elevate<CExperimentObjectMap, CCopasiParameterGroup>(getGroup("Object Map"));
304 
305  if (!mpObjectMap) return false;
306 
307  CCopasiParameterGroup *pGroup = getGroup("Column Role");
308 
309  if (pGroup) // We have an old data format
310  {
311  size_t i, imax = pGroup->size();
312  CExperimentObjectMap Roles;
313  Roles.setNumCols(imax);
314 
315  for (i = 0; i < imax; i++)
316  {
317  Roles.setRole(i, *(Type *) pGroup->getValue(StringPrint("%d", i)).pUINT);
318  Roles.setObjectCN(i, mpObjectMap->getObjectCN(i));
319  }
320 
321  mpObjectMap->clear();
322  *mpObjectMap = Roles;
323  mpObjectMap =
324  elevate<CExperimentObjectMap, CCopasiParameterGroup>(getGroup("Object Map"));
325 
326  removeParameter("Column Role");
327 
328  *mpWeightMethod = SD;
329  }
330 
332 
333  return true;
334 }
335 
337 {
338  size_t i, imax = mpObjectMap->size();
339 
341  CFittingPoint * pPoint;
342 
343  for (i = 0; i < imax; i++)
344  if (mpObjectMap->getRole(i) == dependent)
345  {
346  pPoint = new CFittingPoint(mpObjectMap->getObjectCN(i));
347  mFittingPoints.add(pPoint, true);
348  }
349 }
350 
351 void CExperiment::updateFittedPointValues(const size_t & index, bool includeSimulation)
352 {
355 
356  if (index >= mNumDataRows ||
358  {
359  for (; it != end; ++it)
360  (*it)->setValues(std::numeric_limits<C_FLOAT64>::quiet_NaN(),
361  std::numeric_limits<C_FLOAT64>::quiet_NaN(),
362  std::numeric_limits<C_FLOAT64>::quiet_NaN(),
363  std::numeric_limits<C_FLOAT64>::quiet_NaN());
364 
365  return;
366  }
367 
368  C_FLOAT64 Independent;
369 
371  Independent = mDataTime[index];
372  else
373  Independent = (C_FLOAT64) index;
374 
375  C_FLOAT64 Residual;
376 
377  C_FLOAT64 * pDataDependentCalculated =
379  C_FLOAT64 * pDataDependent = mDataDependent[index];
380  C_FLOAT64 * pScale = mScale[index];
381 
382  for (; it != end; ++it, ++pScale, ++pDataDependentCalculated, ++pDataDependent)
383  {
384 #ifdef COPASI_PARAMETERFITTING_RESIDUAL_SCALING
385  Residual = (*pDataDependentCalculated - *pDataDependent) / (*pDataDependentCalculated > 1 ? *pDataDependentCalculated : 1.0);
386 #else
387  Residual = (*pDataDependentCalculated - *pDataDependent) * *pScale;
388 #endif
389 
390  (*it)->setValues(Independent,
391  *pDataDependent,
392  includeSimulation ? *pDataDependentCalculated : std::numeric_limits<C_FLOAT64>::quiet_NaN(),
393  Residual);
394  }
395 
396  return;
397 }
398 
400 {
403 
404  if (index >= extendedTimeSeriesSize())
405  {
406  for (; it != end; ++it)
407  (*it)->setValues(std::numeric_limits<C_FLOAT64>::quiet_NaN(),
408  std::numeric_limits<C_FLOAT64>::quiet_NaN(),
409  std::numeric_limits<C_FLOAT64>::quiet_NaN(),
410  std::numeric_limits<C_FLOAT64>::quiet_NaN());
411 
412  return;
413  }
414 
415  size_t i;
416 
417  for (i = 1; it != end; ++it, ++i)
418  {
419  (*it)->setValues(mExtendedTimeSeries[index * (mDataDependent.numCols() + 1)],
420  std::numeric_limits<C_FLOAT64>::quiet_NaN(),
421  mExtendedTimeSeries[index * (mDataDependent.numCols() + 1) + i] ,
422  std::numeric_limits<C_FLOAT64>::quiet_NaN());
423  }
424 }
425 
427  C_FLOAT64 *& residuals) const
428 {
429  C_FLOAT64 Residual;
430  C_FLOAT64 s = 0.0;
431 
432  C_FLOAT64 const * pDataDependent = mDataDependent[index];
433  C_FLOAT64 const * pEnd = pDataDependent + mDataDependent.numCols();
434  C_FLOAT64 * const * ppDependentValues = mDependentValues.array();
435  C_FLOAT64 const * pScale = mScale[index];
436 
437  std::vector< Refresh * >::const_iterator it = mRefreshMethods.begin();
438  std::vector< Refresh * >::const_iterator end = mRefreshMethods.end();
439 
440  for (; it != end; ++it)
441  (**it)();
442 
443  if (mMissingData)
444  {
445  if (residuals)
446  for (; pDataDependent != pEnd;
447  pDataDependent++, ppDependentValues++, pScale++, residuals++)
448  {
449  if (isnan(*pDataDependent))
450  {
451  // We ignore missing data, i.e., the residual is 0.
452  *residuals = 0.0;
453  continue;
454  }
455 
456 #ifdef COPASI_PARAMETERFITTING_RESIDUAL_SCALING
457  *residuals = (*pDataDependent - **ppDependentValues) / std::max(1.0, **ppDependentValues);
458 #else
459  *residuals = (*pDataDependent - **ppDependentValues) * *pScale;
460 #endif
461 
462  s += *residuals * *residuals;
463  }
464  else
465  for (; pDataDependent != pEnd;
466  pDataDependent++, ppDependentValues++, pScale++)
467  {
468  if (isnan(*pDataDependent)) continue;
469 
470 #ifdef COPASI_PARAMETERFITTING_RESIDUAL_SCALING
471  Residual = (*pDataDependent - **ppDependentValues) / std::max(1.0, **ppDependentValues);
472 #else
473  Residual = (*pDataDependent - **ppDependentValues) * *pScale;
474 #endif
475 
476  s += Residual * Residual;
477  }
478  }
479  else
480  {
481  if (residuals)
482  for (; pDataDependent != pEnd;
483  pDataDependent++, ppDependentValues++, pScale++, residuals++)
484  {
485 #ifdef COPASI_PARAMETERFITTING_RESIDUAL_SCALING
486  *residuals = (*pDataDependent - **ppDependentValues) / std::max(1.0, **ppDependentValues);
487 #else
488  *residuals = (*pDataDependent - **ppDependentValues) * *pScale;
489 #endif
490 
491  s += *residuals * *residuals;
492  }
493  else
494  for (; pDataDependent != pEnd;
495  pDataDependent++, ppDependentValues++, pScale++)
496  {
497 
498 #ifdef COPASI_PARAMETERFITTING_RESIDUAL_SCALING
499  Residual = (*pDataDependent - **ppDependentValues) / std::max(1.0, **ppDependentValues);
500 #else
501  Residual = (*pDataDependent - **ppDependentValues) * *pScale;
502 #endif
503 
504  s += Residual * Residual;
505  }
506  }
507 
508  return s;
509 }
510 
512  C_FLOAT64 *& dependentValues)
513 {
514  if (index == 0)
515  mpDataDependentCalculated = dependentValues;
516 
517  C_FLOAT64 Residual;
518  C_FLOAT64 s = 0.0;
519 
520  C_FLOAT64 const * pDataDependent = mDataDependent[index];
521  C_FLOAT64 const * pEnd = pDataDependent + mDataDependent.numCols();
522  C_FLOAT64 * const * ppDependentValues = mDependentValues.array();
523  C_FLOAT64 const * pScale = mScale[index];
524 
525  std::vector< Refresh * >::const_iterator it = mRefreshMethods.begin();
526  std::vector< Refresh * >::const_iterator end = mRefreshMethods.end();
527 
528  for (; it != end; ++it)
529  (**it)();
530 
531  if (mMissingData)
532  {
533  for (; pDataDependent != pEnd;
534  pDataDependent++, ppDependentValues++, pScale++, dependentValues++)
535  {
536  *dependentValues = **ppDependentValues;
537 
538  if (isnan(*pDataDependent)) continue;
539 
540 #ifdef COPASI_PARAMETERFITTING_RESIDUAL_SCALING
541  Residual = (*pDataDependent - *dependentValues) / std::max(1.0, *dependentValues);
542 #else
543  Residual = (*pDataDependent - *dependentValues) * *pScale;
544 #endif
545 
546  s += Residual * Residual;
547  }
548  }
549  else
550  {
551  for (; pDataDependent != pEnd;
552  pDataDependent++, ppDependentValues++, pScale++, dependentValues++)
553  {
554  *dependentValues = **ppDependentValues;
555 
556 #ifdef COPASI_PARAMETERFITTING_RESIDUAL_SCALING
557  Residual = (*pDataDependent - *dependentValues) / std::max(1.0, *dependentValues);
558 #else
559  Residual = (*pDataDependent - *dependentValues) * *pScale;
560 #endif
561 
562  s += Residual * Residual;
563  }
564  }
565 
566  return s;
567 }
568 
570 {
572  mExtendedTimeSeries.resize(s * (this->getDependentData().numCols() + 1)); //+1 for time
574 }
575 
577 {
578  //first store time
579  *mStorageIt = time; ++mStorageIt;
580 
581  //do all necessary refreshs
582  std::vector< Refresh * >::const_iterator it = mRefreshMethods.begin();
583  std::vector< Refresh * >::const_iterator end = mRefreshMethods.end();
584 
585  for (; it != end; ++it)
586  (**it)();
587 
588  //store the calculated data
589  C_FLOAT64 * const * ppDependentValues = mDependentValues.array();
590  size_t i, imax = mDataDependent.numCols();
591 
592  for (i = 0; i < imax; ++i, ++ppDependentValues, ++mStorageIt)
593  *mStorageIt = **ppDependentValues;
594 }
595 
597 {
599 }
600 
602 {
603  C_FLOAT64 * pTime = NULL;
604  C_FLOAT64 SavedTime = 0.0;
605 
607  {
608  pTime = const_cast<C_FLOAT64 *>(&getObjectDataModel()->getModel()->getTime());
609  SavedTime = *pTime;
610  }
611 
612  size_t numRows = mDataDependent.numRows();
613  size_t numCols = mDataDependent.numCols();
614 
615  // Overall statistic;
616  mMean = 0.0;
617  mMeanSD = 0.0;
618  mObjectiveValue = 0.0;
619  mRMS = 0.0;
620  size_t ValidValueCount = 0;
621 
622  // per row statistic;
623  mRowObjectiveValue.resize(numRows);
624  mRowObjectiveValue = 0.0;
625  mRowRMS.resize(numRows);
626  mRowRMS = 0.0;
627  CVector< size_t > RowCount;
628  RowCount.resize(numRows);
629  RowCount = 0;
630 
631  // per column statistic;
632  mColumnObjectiveValue.resize(numCols);
633  mColumnObjectiveValue = 0.0;
634  mColumnRMS.resize(numCols);
635  mColumnRMS = 0.0;
638 
639  size_t i, j;
640  C_FLOAT64 Residual;
641 
642  if (mpDataDependentCalculated == NULL)
643  return false;
644 
645  C_FLOAT64 * pDataDependentCalculated = mpDataDependentCalculated;
646  C_FLOAT64 * pDataDependent = mDataDependent.array();
647  C_FLOAT64 * pScale = mScale.array();
648 
649  for (i = 0; i < numRows; i++)
650  {
651  for (j = 0; j < numCols; j++, pDataDependentCalculated++, pDataDependent++, pScale++)
652  {
653 
654 #ifdef COPASI_PARAMETERFITTING_RESIDUAL_SCALING
655  Residual = (*pDataDependentCalculated - *pDataDependent) * std::max(1.0, *pDataDependentCalculated);
656 #else
657  Residual = (*pDataDependentCalculated - *pDataDependent) * *pScale;
658 #endif
659 
660  if (isnan(Residual)) continue;
661 
662  mMean += Residual;
663 
664  Residual = Residual * Residual;
665 
666  mObjectiveValue += Residual;
667  ValidValueCount++;
668 
669  mRowObjectiveValue[i] += Residual;
670  RowCount[i]++;
671 
672  mColumnObjectiveValue[j] += Residual;
674  }
675  }
676 
677  if (ValidValueCount)
678  {
679  mMean /= ValidValueCount;
680  mRMS = sqrt(mObjectiveValue / ValidValueCount);
681  }
682  else
683  {
684  mMean = std::numeric_limits<C_FLOAT64>::quiet_NaN();
685  mRMS = std::numeric_limits<C_FLOAT64>::quiet_NaN();
686  }
687 
688  for (i = 0; i < numRows; i++)
689  {
690  if (RowCount[i])
691  mRowRMS[i] = sqrt(mRowObjectiveValue[i] / RowCount[i]);
692  else
693  mRowRMS[i] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
694  }
695 
696  for (j = 0; j < numCols; j++)
697  {
698  if (mColumnValidValueCount[j])
700  else
701  mColumnRMS[j] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
702  }
703 
704  pDataDependentCalculated = mpDataDependentCalculated;
705  pDataDependent = mDataDependent.array();
706  pScale = mScale.array();
707 
708  for (i = 0; i < numRows; i++)
709  {
710  for (j = 0; j < numCols; j++, pDataDependentCalculated++, pDataDependent++, pScale++)
711  {
712 
713 #ifdef COPASI_PARAMETERFITTING_RESIDUAL_SCALING
714  Residual = mMean - (*pDataDependentCalculated - *pDataDependent) / std::max(1.0, *pDataDependentCalculated);
715 #else
716  Residual = mMean - (*pDataDependentCalculated - *pDataDependent) * *pScale;
717 #endif
718 
719  if (isnan(Residual)) continue;
720 
721  mMeanSD += Residual * Residual;
722  }
723  }
724 
725  if (ValidValueCount)
726  mMeanSD = sqrt(mMeanSD / ValidValueCount);
727  else
728  mMeanSD = std::numeric_limits<C_FLOAT64>::quiet_NaN();
729 
730  if (*mpTaskType == CCopasiTask::timeCourse) *pTime = SavedTime;
731 
732  return true;
733 }
734 
735 bool CExperiment::compile(const std::vector< CCopasiContainer * > listOfContainer)
736 {
737  bool success = true;
738 
739  if (!mpObjectMap->compile(listOfContainer))
740  success = false;
741 
742  size_t LastMappedColumn = mpObjectMap->getLastColumn();
744 
745  size_t i, imax = mpObjectMap->getLastNotIgnoredColumn();
746 
747  if (*mpNumColumns < imax)
748  *mpNumColumns = imax;
749 
750  //if (*mpNumColumns <= imax)
751  // {
752  // CCopasiMessage(CCopasiMessage::ERROR, MCFitting + 4, imax + 1, *mpNumColumns + 1);
753  // return false; // More column types specified than we have data columns
754  // }
755 
756  if (LastMappedColumn < imax || LastMappedColumn == C_INVALID_INDEX)
757  {
759  return false; // More column types specified than we have mapped columns
760  }
761 
762  size_t IndependentCount = mDataIndependent.numCols();
763  size_t DependentCount = mDataDependent.numCols();
764 
765  mDependentValues.resize(DependentCount);
766  mIndependentUpdateMethods.resize(IndependentCount);
767  mIndependentValues.resize(IndependentCount);
768  mIndependentObjects.clear();
769  mDependentObjects.clear();
770  std::set< const CCopasiObject * > Dependencies;
771 
772  IndependentCount = 0;
773  DependentCount = 0;
774 
775  bool TimeFound = false;
776 
777  for (i = 0; i <= imax; i++)
778  switch (mpObjectMap->getRole(i))
779  {
780  case ignore:
781  break;
782 
783  case independent:
784 
785  if (!Objects[i]) // Object not found
786  {
788  return false;
789  }
790 
791  if (!Objects[i]->isValueDbl())
792  {
793  CCopasiMessage(CCopasiMessage::ERROR, MCFitting + 6, Objects[i]->getObjectDisplayName().c_str(), i + 1);
794  return false;
795  }
796 
797  mIndependentObjects.insert(Objects[i]);
798  mIndependentUpdateMethods[IndependentCount] =
799  Objects[i]->getUpdateMethod();
800  mIndependentValues[IndependentCount] =
801  *(C_FLOAT64 *)Objects[i]->getValuePointer();
802  // :TODO: do we have to check if getValuePointer() return a valid pointer?
803 
804  IndependentCount++;
805  break;
806 
807  case dependent:
808 
809  if (!Objects[i]) // Object not found
810  {
812  return false;
813  }
814 
815  if (!Objects[i]->isValueDbl())
816  {
817  CCopasiMessage(CCopasiMessage::ERROR, MCFitting + 6, Objects[i]->getObjectDisplayName().c_str(), i + 1);
818  return false;
819  }
820 
821  mDependentValues[DependentCount] =
822  (C_FLOAT64 *) Objects[i]->getValuePointer();
823  // :TODO: do we have to check if getValuePointer() return a valid pointer?
824  mDependentObjects[Objects[i]] = DependentCount;
825  mColumnScale[DependentCount] = mpObjectMap->getScale(i);
826  Dependencies.insert(Objects[i]->getValueObject());
827 
828  DependentCount++;
829  break;
830 
831  case time:
832  TimeFound = true;
833  break;
834  }
835 
836  /* We need to check whether a column is mapped to time */
837  if (!TimeFound && *mpTaskType == CCopasiTask::timeCourse)
838  success = false;
839 
840  // Allocation and initialization of statistical information
841  size_t numRows = mDataDependent.numRows();
842  size_t numCols = mDataDependent.numCols();
843 
844  // Overall statistic;
845  mMean = std::numeric_limits<C_FLOAT64>::quiet_NaN();
846  mMeanSD = std::numeric_limits<C_FLOAT64>::quiet_NaN();
847  mObjectiveValue = std::numeric_limits<C_FLOAT64>::quiet_NaN();
848  mRMS = std::numeric_limits<C_FLOAT64>::quiet_NaN();
849 
850  // per row statistic;
851  mRowObjectiveValue.resize(numRows);
852  mRowObjectiveValue = std::numeric_limits<C_FLOAT64>::quiet_NaN();
853  mRowRMS.resize(numRows);
854  mRowRMS = std::numeric_limits<C_FLOAT64>::quiet_NaN();
855 
856  // per column statistic;
857  mColumnObjectiveValue.resize(numCols);
858  mColumnObjectiveValue = std::numeric_limits<C_FLOAT64>::quiet_NaN();
859  mColumnRMS.resize(numCols);
860  mColumnRMS = std::numeric_limits<C_FLOAT64>::quiet_NaN();
862  mColumnValidValueCount = std::numeric_limits<size_t>::quiet_NaN();
863 
864  CModel * pModel =
865  dynamic_cast< CModel * >(getObjectDataModel()->ObjectFromName(listOfContainer, CCopasiObjectName("Model=" + CCopasiObjectName::escape(getObjectDataModel()->getModel()->getObjectName()))));
866 
868 
870 
871  return success;
872 }
873 
874 bool CExperiment::read(std::istream & in,
875  size_t & currentLine)
876 {
877  // Allocate for reading
878  size_t i, imax = mpObjectMap->size();
879 
880  if (*mpNumColumns < imax)
881  *mpNumColumns = imax;
882  //if (*mpNumColumns < imax)
883  // {
884  // CCopasiMessage(CCopasiMessage::ERROR, MCFitting + 4, imax, *mpNumColumns);
885  // return false; // More column types specified than we have data columns
886  // }
887 
888  size_t IndependentCount = 0;
889  size_t DependentCount = 0;
890  size_t TimeCount = 0;
891  size_t IgnoreCount = 0;
892 
893  for (i = 0; i < imax; i++)
894  switch (mpObjectMap->getRole(i))
895  {
896  case ignore:
897  IgnoreCount++;
898  break;
899 
900  case independent:
901  IndependentCount++;
902  break;
903 
904  case dependent:
905  DependentCount++;
906  break;
907 
908  case time:
909  TimeCount++;
910  break;
911  }
912 
914  ((*mpHeaderRow < *mpFirstRow || *mpLastRow < *mpHeaderRow) ? 1 : 0);
915 
916  mDataTime.resize(TimeCount ? mNumDataRows : 0);
917  mDataIndependent.resize(mNumDataRows, IndependentCount);
918  mDataDependent.resize(mNumDataRows, DependentCount);
920  mColumnName.resize(IndependentCount + DependentCount + TimeCount);
921 
922  // resize the vectors for the statistics
923  mMeans.resize(DependentCount);
924  mColumnScale.resize(DependentCount);
925  mDefaultColumnScale.resize(DependentCount);
926  mColumnValidValueCount.resize(DependentCount);
927 
928  if (!TimeCount && *mpTaskType == CCopasiTask::timeCourse)
929  {
931  return false;
932  }
933 
934  if (DependentCount == 0)
935  {
937  return false;
938  }
939 
940  if (mNumDataRows == 0)
941  {
943  return false;
944  }
945 
946  CTableRow Row(*mpNumColumns, (*mpSeparator)[0]);
947  const std::vector< CTableCell > & Cells = Row.getCells();
948 
949  size_t j;
950 
951  if (currentLine > *mpFirstRow) return false; // We are past our first line
952 
953  // forwind to our first line
954  for (j = currentLine; j < *mpFirstRow && !in.fail(); j++)
955  {
956  skipLine(in);
957  currentLine++;
958  }
959 
960  for (j = 0; j < mNumDataRows && !in.fail(); j++, currentLine++)
961  {
962  in >> Row;
963 
964  if (currentLine == *mpHeaderRow)
965  {
966  j--;
967 
968  size_t Column = 0;
969 
970  for (i = 0; i < *mpNumColumns; i++)
971  if (mpObjectMap->getRole(i) != ignore)
972  mColumnName[Column++] = Cells[i].getName();
973 
974  continue;
975  }
976 
977  IndependentCount = 0;
978  DependentCount = 0;
979 
980  for (i = 0; i < imax; i++)
981  {
982  switch (mpObjectMap->getRole(i))
983  {
984  case ignore:
985  break;
986 
987  case independent:
988 
989  if (!Cells[i].isValue())
990  {
992  getObjectName().c_str(), currentLine, i + 1);
993  return false;
994  }
995 
996  mDataIndependent[j][IndependentCount++] =
997  Cells[i].getValue();
998  break;
999 
1000  case dependent:
1001  mDataDependent[j][DependentCount++] =
1002  Cells[i].getValue();
1003  break;
1004 
1005  case time:
1006 
1007  if (!Cells[i].isValue())
1008  {
1010  getObjectName().c_str(), currentLine, i + 1);
1011  return false;
1012  }
1013 
1014  mDataTime[j] = Cells[i].getValue();
1015  break;
1016  }
1017  }
1018  }
1019 
1020  if ((in.fail() && !in.eof()))
1021  {
1023  return false;
1024  }
1025 
1026  if (j != mNumDataRows)
1027  {
1028  CCopasiMessage(CCopasiMessage::ERROR, MCFitting + 7, mNumDataRows, j - 1);
1029  return false;
1030  }
1031 
1032  // If it is a time course this is the place to assert that it is sorted.
1034  {
1035  CVector<size_t> Pivot;
1037 
1038  mDataTime.applyPivot(Pivot);
1040  mDataDependent.applyPivot(Pivot);
1041 
1042  for (mNumDataRows--; mNumDataRows != C_INVALID_INDEX; mNumDataRows--)
1043  if (!isnan(mDataTime[mNumDataRows])) break;
1044 
1045  mNumDataRows++;
1046  }
1047 
1048  return calculateWeights();
1049 }
1050 
1052 {
1053  // We need to calculate the means and the weights
1055 
1056  size_t DependentCount = mMeans.size();
1057  CVector< C_FLOAT64 > MeanSquares(DependentCount);
1058  CVector< C_FLOAT64 > ColumnEpsilons(DependentCount);
1059  C_FLOAT64 * pColumnEpsilon;
1060  size_t i, j;
1061 
1062  mMeans = 0.0;
1063  MeanSquares = 0.0;
1064  ColumnEpsilons = std::numeric_limits< C_FLOAT64 >::infinity();
1065 
1067  mMissingData = false;
1068 
1069  // Calculate the means
1070  for (i = 0; i < mNumDataRows; i++)
1071  for (j = 0, pColumnEpsilon = ColumnEpsilons.array(); j < DependentCount; j++, pColumnEpsilon++)
1072  {
1073  C_FLOAT64 & Data = mDataDependent[i][j];
1074 
1075  if (!isnan(Data))
1076  {
1078  mMeans[j] += Data;
1079  MeanSquares[j] += Data * Data;
1080 
1081  if (Data != 0.0 && fabs(Data) < *pColumnEpsilon)
1082  {
1083  *pColumnEpsilon = fabs(Data);
1084  }
1085  }
1086  else
1087  mMissingData = true;
1088  }
1089 
1090  // calculate the means;
1091  for (j = 0, pColumnEpsilon = ColumnEpsilons.array(); j < DependentCount; j++, pColumnEpsilon++)
1092  {
1093  if (*pColumnEpsilon == std::numeric_limits< C_FLOAT64 >::infinity())
1094  {
1095  // TODO CRITICAL We need to create a preference for this.
1096  *pColumnEpsilon = 1e8 * std::numeric_limits< C_FLOAT64 >::epsilon();
1097  }
1098 
1099  if (mColumnValidValueCount[j])
1100  {
1101  mMeans[j] /= mColumnValidValueCount[j];
1102  MeanSquares[j] /= mColumnValidValueCount[j];
1103  }
1104  else
1105  {
1106  // :TODO: We need to create an error message here to instruct the user
1107  // :TODO: to mark this column as ignored as it contains no data.
1108  mMeans[j] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1109  MeanSquares[j] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1110  }
1111  }
1112 
1113  for (i = 0; i < DependentCount; i++)
1114  {
1115  C_FLOAT64 & DefaultColumScale = mDefaultColumnScale[i];
1116 
1117  switch (*mpWeightMethod)
1118  {
1119  case SD:
1120  DefaultColumScale = MeanSquares[i] - mMeans[i] * mMeans[i];
1121  break;
1122 
1123  case MEAN:
1124  DefaultColumScale = mMeans[i] * mMeans[i];
1125  break;
1126 
1127  case MEAN_SQUARE:
1128  DefaultColumScale = MeanSquares[i];
1129  break;
1130 
1131  case VALUE_SCALING:
1132  DefaultColumScale = ColumnEpsilons[i] * ColumnEpsilons[i] * 1e-12;
1133  break;
1134  }
1135 
1136  if (DefaultColumScale < MinWeight) MinWeight = DefaultColumScale;
1137  }
1138 
1140  {
1141  MinWeight = 1.0;
1142  }
1143 
1144  if (*mpWeightMethod != VALUE_SCALING)
1145  {
1146  // We have to calculate the default weights
1147  for (i = 0; i < DependentCount; i++)
1148  mDefaultColumnScale[i] =
1149  (MinWeight + sqrt(std::numeric_limits< C_FLOAT64 >::epsilon()))
1150  / (mDefaultColumnScale[i] + sqrt(std::numeric_limits< C_FLOAT64 >::epsilon()));
1151  }
1152 
1153  return true;
1154 }
1155 
1156 const std::map< CCopasiObject *, size_t > & CExperiment::getDependentObjects() const
1157 {return mDependentObjects;}
1158 
1160 {
1161  mColumnName.resize(*mpNumColumns);
1162 
1163  if (*mpHeaderRow == InvalidIndex) return false;
1164 
1165  // Open the file
1166  std::ifstream in;
1167  in.open(CLocaleString::fromUtf8(getFileName()).c_str(), std::ios::binary);
1168 
1169  if (in.fail()) return false;
1170 
1171  // Forwind to header row.
1172  size_t i;
1173 
1174  for (i = 1; i < *mpHeaderRow && !in.fail(); i++)
1175  skipLine(in);
1176 
1177  // Read row
1178  CTableRow Row(*mpNumColumns, (*mpSeparator)[0]);
1179  const std::vector< CTableCell > & Cells = Row.getCells();
1180  in >> Row;
1181 
1182  if (in.fail() && !in.eof()) return false;
1183 
1184  for (i = 0; i < *mpNumColumns; i++)
1185  mColumnName[i] = Cells[i].getName();
1186 
1187  return true;
1188 }
1189 
1191 {
1192  size_t tmp, count = 0;
1193 
1194  std::ifstream in;
1195  in.open(CLocaleString::fromUtf8(getFileName()).c_str(), std::ios::binary);
1196 
1197  if (in.fail()) return false;
1198 
1199  // Forwind to first row.
1200  size_t i;
1201 
1202  for (i = 1; i < *mpFirstRow && !in.fail(); i++)
1203  skipLine(in);
1204 
1205  CTableRow Row(0, (*mpSeparator)[0]);
1206 
1207  for (i--; i < *mpLastRow; i++)
1208  if ((tmp = Row.guessColumnNumber(in, false)) > count)
1209  count = tmp;
1210 
1211  return count;
1212 }
1213 
1214 const std::vector< std::string > & CExperiment::getColumnNames() const
1215 {return mColumnName;}
1216 
1218 {
1219  size_t i, imax = mIndependentUpdateMethods.size();
1220 
1221  for (i = 0; i < imax; i++)
1223 
1224  return true;
1225 }
1226 
1228 {
1229  size_t i, imax = mIndependentUpdateMethods.size();
1230 
1231  for (i = 0; i < imax; i++)
1233 
1234  return true;
1235 }
1236 
1238 {return *mpTaskType;}
1239 
1241 {
1242  switch (type)
1243  {
1246  *mpTaskType = type;
1247  return true;
1248  break;
1249 
1250  default:
1251  break;
1252  }
1253 
1254  return false;
1255 }
1256 
1258 {
1260 }
1261 
1263 {
1266  else
1267  return true;
1268 }
1269 
1271 {return mDataTime;}
1272 
1274 {return mDataIndependent;}
1275 
1277 {return mDataDependent;}
1278 
1279 const std::string & CExperiment::getFileName() const
1280 {
1281  std::string * pFileName = const_cast<CExperiment *>(this)->mpFileName;
1282 
1283  if (CDirEntry::isRelativePath(*pFileName) &&
1284  !CDirEntry::makePathAbsolute(*pFileName,
1286  *pFileName = CDirEntry::fileName(*pFileName);
1287 
1288  return *mpFileName;
1289 }
1290 
1291 bool CExperiment::setFileName(const std::string & fileName)
1292 {
1293  *mpFileName = fileName;
1294  return true;
1295 }
1296 
1298 {return * mpObjectMap;}
1299 
1301 {return mFittingPoints;}
1302 
1303 const unsigned C_INT32 & CExperiment::getFirstRow() const
1304 {return *mpFirstRow;}
1305 
1306 bool CExperiment::setFirstRow(const unsigned C_INT32 & first)
1307 {
1308  if (first > *mpLastRow ||
1309  (first == *mpLastRow && first == *mpHeaderRow)) return false;
1310 
1311  *mpFirstRow = first;
1312  return true;
1313 }
1314 
1315 const unsigned C_INT32 & CExperiment::getLastRow() const
1316 {return *mpLastRow;}
1317 
1318 bool CExperiment::setLastRow(const unsigned C_INT32 & last)
1319 {
1320  if (*mpFirstRow > last ||
1321  (*mpFirstRow == last && last == *mpHeaderRow)) return false;
1322 
1323  *mpLastRow = last;
1324  return true;
1325 }
1326 
1327 const unsigned C_INT32 & CExperiment::getHeaderRow() const
1328 {return *mpHeaderRow;}
1329 
1330 bool CExperiment::setHeaderRow(const unsigned C_INT32 & header)
1331 {
1332  if (header == *mpFirstRow && header == *mpLastRow) return false;
1333 
1334  *mpHeaderRow = header;
1335  return true;
1336 }
1337 
1338 const unsigned C_INT32 & CExperiment::getNumColumns() const
1339 {return *mpNumColumns;}
1340 
1341 bool CExperiment::setNumColumns(const unsigned C_INT32 & cols)
1342 {
1343  *mpNumColumns = cols;
1344  return true;
1345 }
1346 
1348 {return mNumDataRows;}
1349 
1350 const std::string & CExperiment::getSeparator() const
1351 {return *mpSeparator;}
1352 
1353 bool CExperiment::setSeparator(const std::string & separator)
1354 {
1355  *mpSeparator = separator;
1356  return true;
1357 }
1358 
1360 {return *mpWeightMethod;}
1361 
1363 {
1364  if (*mpWeightMethod == weightMethod) return true;
1365 
1366  // Reset to default weights
1367  *mpWeightMethod = weightMethod;
1368  std::vector< CCopasiParameter * >::iterator it = mpObjectMap->CCopasiParameter::getValue().pGROUP->begin();
1369  std::vector< CCopasiParameter * >::iterator end = mpObjectMap->CCopasiParameter::getValue().pGROUP->end();
1370 
1371  for (; it != end; ++ it)
1372  static_cast< CExperimentObjectMap::CDataColumn * >(*it)->setScale(std::numeric_limits<C_FLOAT64>::quiet_NaN());
1373 
1374  return true;
1375 }
1376 
1377 const bool & CExperiment::isRowOriented() const
1378 {return *mpRowOriented;}
1379 
1380 bool CExperiment::setIsRowOriented(const bool & isRowOriented)
1381 {
1383  return true;
1384 }
1385 
1387  const CExperiment * rhs)
1388 {
1389  return (*lhs->mpFileName < *rhs->mpFileName ||
1390  (*lhs->mpFileName == *rhs->mpFileName &&
1391  *lhs->mpFirstRow < *rhs->mpFirstRow));
1392 }
1393 
1394 bool operator == (const CExperiment & lhs,
1395  const CExperiment & rhs)
1396 {
1397  std::string Key = *lhs.getValue("Key").pKEY;
1398  const_cast<CExperiment *>(&lhs)->setValue("Key", *rhs.getValue("Key").pKEY);
1399 
1400  bool Result =
1401  (*static_cast<const CCopasiParameterGroup *>(&lhs) ==
1402  *static_cast<const CCopasiParameterGroup *>(&rhs));
1403 
1404  const_cast<CExperiment *>(&lhs)->setValue("Key", Key);
1405 
1406  return Result;
1407 }
1408 
1409 void CExperiment::printResult(std::ostream * ostream) const
1410 {
1411  std::ostream & os = *ostream;
1412 
1413  os << "File Name:\t" << getFileName() << std::endl;
1414  os << "Experiment:\t" << getObjectName() << std::endl;
1415 
1416  os << "Mean:\t" << mMean << std::endl;
1417  os << "Objective Value:\t" << mObjectiveValue << std::endl;
1418  os << "Root Mean Square:\t" << mRMS << std::endl;
1419 
1420  size_t i, imax = mNumDataRows;
1421  size_t j, jmax = mDataDependent.numCols();
1422  size_t k, kmax = mpObjectMap->getLastNotIgnoredColumn() + 1;
1423 
1424  const CVector<CCopasiObject *> & Objects =
1426 
1427  os << "Row\t";
1428 
1430  os << "Time\t";
1431 
1432  for (k = 0; k < kmax; k++)
1434  {
1435  std::string Name;
1436 
1437  if (k < Objects.size() && Objects[k] != NULL)
1438  Name = Objects[k]->getObjectDisplayName();
1439  else
1440  Name = "unknown";
1441 
1442  os << Name << "(Data)\t";
1443  os << Name << "(Fit)\t";
1444  os << Name << "(Weighted Error)\t";
1445  }
1446 
1447  os << "Objective Value\tRoot Mean Square" << std::endl << std::endl;
1448 
1449  C_FLOAT64 * pDataDependentCalculated = mpDataDependentCalculated;
1450 
1451  if (pDataDependentCalculated)
1452  for (i = 0; i < imax; i++)
1453  {
1454  os << i + 1 << ".\t";
1455 
1457  os << mDataTime[i] << "\t";
1458 
1459  for (j = 0; j < jmax; j++, pDataDependentCalculated++)
1460  {
1461  os << mDataDependent(i, j) << "\t";
1462  os << *pDataDependentCalculated << "\t";
1463  os << mScale(i, j) *(*pDataDependentCalculated - mDataDependent(i, j)) << "\t";
1464  }
1465 
1466  os << mRowObjectiveValue[i] << "\t" << mRowRMS[i] << std::endl;
1467  }
1468  else
1469  for (i = 0; i < imax; i++)
1470  {
1471  os << i + 1 << ".\t";
1472 
1474  os << mDataTime[i] << "\t";
1475 
1476  for (j = 0; j < jmax; j++)
1477  {
1478  os << mDataDependent(i, j) << "\tNaN\tNaN\t";
1479  }
1480 
1481  if (i < mRowObjectiveValue.size())
1482  {
1483  os << mRowObjectiveValue[i] << "\t";
1484  }
1485  else
1486  {
1487  os << "NaN\t";
1488  }
1489 
1490  if (i < mRowRMS.size())
1491  {
1492  os << mRowRMS[i];
1493  }
1494  else
1495  {
1496  os << "NaN";
1497  }
1498 
1499  os << std::endl;
1500  }
1501 
1502  os << "Objective Value";
1503 
1505  os << "\t";
1506 
1507  for (j = 0; j < jmax; j++)
1508  {
1509  if (j < mColumnObjectiveValue.size())
1510  os << "\t\t\t" << mColumnObjectiveValue[j];
1511  else
1512  os << "\t\t\tNaN";
1513  }
1514 
1515  os << std::endl;
1516 
1517  os << "Root Mean Square";
1518 
1520  os << "\t";
1521 
1522  for (j = 0; j < jmax; j++)
1523  {
1524  if (j < mColumnRMS.size())
1525  os << "\t\t\t" << mColumnRMS[j];
1526  else
1527  os << "\t\t\tNaN";
1528  }
1529 
1530  os << std::endl;
1531 
1532  os << "Weight";
1533 
1535  os << "\t";
1536 
1537  for (j = 0; j < jmax; j++)
1538  {
1539  if (j < mColumnScale.size())
1540  os << "\t\t\t" << mColumnScale[j];
1541  else
1542  os << "\t\t\tNaN";
1543  }
1544 
1545  os << std::endl;
1546 
1547  return;
1548 }
1549 
1551 {return mObjectiveValue;}
1552 
1554 {return mRMS;}
1555 
1557 {return mMean;}
1558 
1560 {return mMeanSD;}
1561 
1563 {
1564  std::map< CCopasiObject *, size_t >::const_iterator it
1565  = mDependentObjects.find(pObject);
1566 
1567  if (it != mDependentObjects.end())
1568  return mColumnObjectiveValue[it->second];
1569  else
1570  return std::numeric_limits<C_FLOAT64>::quiet_NaN();
1571 }
1572 
1574 {
1575  std::map< CCopasiObject *, size_t>::const_iterator it
1576  = mDependentObjects.find(const_cast<CCopasiObject*>(pObject));
1577 
1578  if (it == mDependentObjects.end())
1579  return std::numeric_limits<C_FLOAT64>::quiet_NaN();
1580 
1581  return mDefaultColumnScale[it->second];
1582 }
1583 
1585 {
1586  std::map< CCopasiObject *, size_t>::const_iterator it
1587  = mDependentObjects.find(pObject);
1588 
1589  if (it != mDependentObjects.end())
1590  return mColumnRMS[it->second];
1591  else
1592  return std::numeric_limits<C_FLOAT64>::quiet_NaN();
1593 }
1594 
1596 {
1597  std::map< CCopasiObject *, size_t>::const_iterator it
1598  = mDependentObjects.find(pObject);
1599 
1600  if (it == mDependentObjects.end() ||
1601  mpDataDependentCalculated == NULL)
1602  return std::numeric_limits<C_FLOAT64>::quiet_NaN();
1603 
1604  C_FLOAT64 Mean = 0;
1605  C_FLOAT64 Residual;
1606  size_t numRows = mDataDependent.numRows();
1607  size_t numCols = mDataDependent.numCols();
1608 
1609  const C_FLOAT64 *pDataDependentCalculated = mpDataDependentCalculated + it->second;
1610  const C_FLOAT64 *pEnd = pDataDependentCalculated + numRows * numCols;
1611  const C_FLOAT64 *pDataDependent = mDataDependent.array() + it->second;
1612  const C_FLOAT64 & Weight = sqrt(mColumnScale[it->second]);
1613 
1614  for (; pDataDependentCalculated != pEnd;
1615  pDataDependentCalculated += numCols, pDataDependent += numCols)
1616  {
1617  Residual = Weight * (*pDataDependentCalculated - *pDataDependent);
1618 
1619  if (isnan(Residual)) continue;
1620 
1621  Mean += Residual;
1622  }
1623 
1624  return Mean;
1625 }
1626 
1628  const C_FLOAT64 & errorMean) const
1629 {
1630  std::map< CCopasiObject *, size_t>::const_iterator it
1631  = mDependentObjects.find(pObject);
1632 
1633  if (it == mDependentObjects.end() ||
1634  mpDataDependentCalculated == NULL)
1635  return std::numeric_limits<C_FLOAT64>::quiet_NaN();
1636 
1637  C_FLOAT64 MeanSD = 0;
1638  C_FLOAT64 Residual;
1639  size_t numRows = mDataDependent.numRows();
1640  size_t numCols = mDataDependent.numCols();
1641 
1642  const C_FLOAT64 *pDataDependentCalculated = mpDataDependentCalculated + it->second;
1643  const C_FLOAT64 *pEnd = pDataDependentCalculated + numRows * numCols;
1644  const C_FLOAT64 *pDataDependent = mDataDependent.array() + it->second;
1645  const C_FLOAT64 *pScale = mScale.array() + it->second;
1646 
1647  for (; pDataDependentCalculated != pEnd;
1648  pDataDependentCalculated += numCols, pDataDependent += numCols)
1649  {
1650 
1651 #ifdef COPASI_PARAMETERFITTING_RESIDUAL_SCALING
1652  Residual = errorMean - (*pDataDependentCalculated - *pDataDependent) / std::max(1.0, *pDataDependentCalculated);
1653 #else
1654  Residual = errorMean - (*pDataDependentCalculated - *pDataDependent) * *pScale;
1655 #endif
1656 
1657  if (isnan(Residual)) continue;
1658 
1659  MeanSD += Residual * Residual;
1660  }
1661 
1662  return MeanSD;
1663 }
1664 
1666 {
1667  std::map< CCopasiObject *, size_t>::const_iterator it
1668  = mDependentObjects.find(pObject);
1669 
1670  if (it != mDependentObjects.end())
1671  return mColumnValidValueCount[it->second];
1672  else
1673  return 0;
1674 }
1675 
1676 const std::set< const CCopasiObject * > & CExperiment::getIndependentObjects() const
1677 {
1678  return mIndependentObjects;
1679 }
1680 
1682 {
1684 
1685  // We need to initialize the scaling factors for the residuals
1686  C_FLOAT64 * pData = mDataDependent.array();
1687  C_FLOAT64 * pScale = mScale.array();
1688  C_FLOAT64 * pScaleEnd = pScale + mScale.size();
1689 
1690  C_FLOAT64 * pColumnScale = mColumnScale.array();
1691  C_FLOAT64 * pColumnScaleEnd = pColumnScale + mColumnScale.size();
1692 
1693  for (; pScale < pScaleEnd;)
1694  {
1695  for (pColumnScale = mColumnScale.array(); pColumnScale < pColumnScaleEnd; ++pColumnScale, ++pScale, ++pData)
1696  {
1697  switch (*mpWeightMethod)
1698  {
1699  case VALUE_SCALING:
1700  *pScale = 1.0 / std::max(fabs(*pData), *pColumnScale);
1701  break;
1702 
1703  default:
1704  *pScale = sqrt(*pColumnScale);
1705  break;
1706  }
1707  }
1708  }
1709 }
1710 
1712 {
1714 }
1715 
1716 /* CFittingPoint Implementation */
1717 
1718 CFittingPoint::CFittingPoint(const std::string & name,
1719  const CCopasiContainer * pParent):
1720  CCopasiContainer("Fitting Point", pParent, "Fitted Point"),
1721  mModelObjectCN(name),
1722  mIndependentValue(std::numeric_limits<C_FLOAT64>::quiet_NaN()),
1723  mMeasuredValue(std::numeric_limits<C_FLOAT64>::quiet_NaN()),
1724  mFittedValue(std::numeric_limits<C_FLOAT64>::quiet_NaN()),
1725  mWeightedError(std::numeric_limits<C_FLOAT64>::quiet_NaN())
1726 {initObjects();}
1727 
1729  const CCopasiContainer * pParent):
1730  CCopasiContainer(src, pParent),
1731  mModelObjectCN(src.mModelObjectCN),
1732  mIndependentValue(src.mIndependentValue),
1733  mMeasuredValue(src.mMeasuredValue),
1734  mFittedValue(src.mFittedValue),
1735  mWeightedError(src.mWeightedError)
1736 {initObjects();}
1737 
1739 
1740 // virtual
1741 std::string CFittingPoint::getObjectDisplayName(bool regular, bool richtext) const
1742 {
1743  const CCopasiDataModel * pDataModel = this->getObjectDataModel();
1744 
1745  if (pDataModel == NULL)
1746  {
1747  return CCopasiContainer::getObjectDisplayName(regular, richtext);
1748  }
1749 
1750  const CCopasiObject * pObject = dynamic_cast< const CCopasiObject * >(pDataModel->getObject(mModelObjectCN));
1751 
1752  if (pObject == NULL)
1753  {
1754  return CCopasiContainer::getObjectDisplayName(regular, richtext);
1755  }
1756 
1757  return pObject->getObjectDisplayName(regular, richtext);
1758 }
1759 
1760 const std::string & CFittingPoint::getModelObjectCN() const
1761 {
1762  return mModelObjectCN;
1763 }
1764 
1765 void CFittingPoint::setValues(const C_FLOAT64 & independent,
1766  const C_FLOAT64 & measured,
1767  const C_FLOAT64 & fitted,
1768  const C_FLOAT64 & weightedError)
1769 {
1770  mIndependentValue = independent;
1771  mMeasuredValue = measured;
1772  mFittedValue = fitted;
1773  mWeightedError = weightedError;
1774 }
1775 
1777 {
1782 }
1783 
1784 std::istream & skipLine(std::istream & in)
1785 {
1786  char c;
1787 
1788  for (in.get(c); c != 0x0a && c != 0x0d; in.get(c))
1789  {
1790  if (in.fail() || in.eof()) break;
1791  }
1792 
1793  // Eat additional line break characters appearing on DOS and Mac text format;
1794  if ((c == 0x0d && in.peek() == 0x0a) || // DOS
1795  (c == 0x0a && in.peek() == 0x0d)) // Mac
1796  in.ignore(1);
1797 
1798  return in;
1799 }
CCopasiDataModel * getObjectDataModel()
bool remove(const std::string &key)
void sortWithPivot(RandomAccessIterator first, RandomAccessIterator last, CVector< size_t > &pivot)
Definition: CSort.h:77
const CMatrix< C_FLOAT64 > & getIndependentData() const
const std::string & getFileName() const
C_FLOAT64 getDefaultScale(const CCopasiObject *const &pObject) const
std::vector< std::string > mColumnName
Definition: CExperiment.h:610
virtual bool elevateChildren()
bool setExperimentType(const CCopasiTask::Type &type)
C_FLOAT64 getScale(const size_t &index) const
virtual std::string getObjectDisplayName(bool regular=true, bool richtext=false) const
virtual const std::string & getName(const size_t &index) const
C_FLOAT64 mWeightedError
Definition: CExperiment.h:58
bool setValue(const CType &value)
bool setScale(const size_t &index, const C_FLOAT64 &scale)
C_FLOAT64 mMean
Definition: CExperiment.h:665
virtual ~CExperiment()
const unsigned C_INT32 & getFirstRow() const
const C_FLOAT64 & getObjectiveValue() const
const CMatrix< C_FLOAT64 > & getDependentData() const
const std::string & getObjectName() const
static bool isRelativePath(const std::string &path)
Definition: CDirEntry.cpp:414
bool compile(const std::vector< CCopasiContainer * > listOfContainer=CCopasiContainer::EmptyList)
CVector< C_FLOAT64 > mDefaultColumnScale
Definition: CExperiment.h:646
CVector< C_FLOAT64 > mColumnRMS
Definition: CExperiment.h:674
CExperiment::Type getRole(const size_t &index) const
const std::vector< CTableCell > & getCells() const
Definition: CTableCell.cpp:141
unsigned C_INT32 * mpFirstRow
Definition: CExperiment.h:565
CFittingPoint(const std::string &name="unknown", const CCopasiContainer *pParent=NULL)
CMatrix< C_FLOAT64 > mDataDependent
Definition: CExperiment.h:630
CVector< C_FLOAT64 > mRowObjectiveValue
Definition: CExperiment.h:670
CVector< C_FLOAT64 * > mDependentValues
Definition: CExperiment.h:648
bool mMissingData
Definition: CExperiment.h:640
iterator begin()
virtual size_t numRows() const
Definition: CMatrix.h:138
bool setIsRowOriented(const bool &isRowOriented)
CVector< size_t > mColumnValidValueCount
Definition: CExperiment.h:675
std::string * mpSeparator
Definition: CExperiment.h:585
bool setHeaderRow(const unsigned C_INT32 &headerRow)
size_t mNumDataRows
Definition: CExperiment.h:658
std::string * mpFileName
Definition: CExperiment.h:560
std::istream & skipLine(std::istream &in)
virtual void * getValuePointer() const
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
static std::string fileName(const std::string &path)
Definition: CDirEntry.cpp:119
bool setFileName(const std::string &fileName)
const std::set< const CCopasiObject * > & getUptoDateObjects() const
Definition: CModel.cpp:1178
#define C_INVALID_INDEX
Definition: copasi.h:222
const std::map< CCopasiObject *, size_t > & getDependentObjects() const
const WeightMethod & getWeightMethod() const
unsigned C_INT32 * mpNumColumns
Definition: CExperiment.h:605
#define InvalidIndex
Definition: CExperiment.cpp:38
#define C_INT32
Definition: copasi.h:90
void initExtendedTimeSeries(size_t s)
static std::string escape(const std::string &name)
const std::set< const CCopasiObject * > & getIndependentObjects() const
size_t guessColumnNumber() const
C_FLOAT64 mObjectiveValue
Definition: CExperiment.h:667
C_FLOAT64 * mpDataDependentCalculated
Definition: CExperiment.h:663
C_FLOAT64 * mStorageIt
Definition: CExperiment.h:688
C_FLOAT64 mFittedValue
Definition: CExperiment.h:57
void setNormalizeWeightsPerExperiment(bool flag)
const CCopasiVector< CFittingPoint > & getFittingPoints() const
CCopasiVector< CFittingPoint > mFittingPoints
Definition: CExperiment.h:682
std::vector< Refresh * > mRefreshMethods
Definition: CExperiment.h:652
size_t getNumDataRows() const
static const char * XMLType[]
Definition: CExperiment.h:84
unsigned C_INT32 * mpLastRow
Definition: CExperiment.h:570
bool removeParameter(const std::string &name)
const std::vector< std::string > & getColumnNames() const
const C_FLOAT64 & getRMS() const
CMatrix< C_FLOAT64 > mDataIndependent
Definition: CExperiment.h:625
iterator end()
const CVector< CCopasiObject * > & getMappedObjects() const
CExperimentObjectMap * mpObjectMap
Definition: CExperiment.h:615
void storeExtendedTimeSeriesData(C_FLOAT64 time)
virtual std::string getObjectDisplayName(bool regular=true, bool richtext=false) const
const C_FLOAT64 & getErrorMean() const
const size_t & getLastColumn() const
static bool compare(const CExperiment *lhs, const CExperiment *rhs)
const std::string & getModelObjectCN() const
bool restoreModelIndependentData()
CVector< C_FLOAT64 > mRowRMS
Definition: CExperiment.h:671
virtual const CCopasiObject * getValueObject() const
static const char * WeightMethodType[]
Definition: CExperiment.h:106
#define MCFitting
bool setSeparator(const std::string &seperator)
std::set< const CCopasiObject * > mIndependentObjects
Definition: CExperiment.h:654
CVector< C_FLOAT64 > mMeans
Definition: CExperiment.h:642
C_FLOAT64 mMeanSD
Definition: CExperiment.h:666
void updateFittedPointValues(const size_t &index, bool includeSimulation)
const C_FLOAT64 & getErrorMeanSD() const
virtual bool add(const CType &src)
bool setRole(const size_t &index, const CExperiment::Type &role)
const Value & getValue() const
CCopasiTask::Type * mpTaskType
Definition: CExperiment.h:575
C_FLOAT64 getErrorSum(CCopasiObject *const &pObject) const
CExperiment(const CCopasiContainer *pParent, const std::string &name="Experiment")
Definition: CExperiment.cpp:76
CExperimentObjectMap & getObjectMap()
bool setValue(const std::string &name, const CType &value)
CRegisteredObjectName mModelObjectCN
Definition: CExperiment.h:54
CVector< C_FLOAT64 > mColumnScale
Definition: CExperiment.h:644
unsigned C_INT32 * pUINT
const unsigned C_INT32 & getHeaderRow() const
CCopasiParameterGroup * assertGroup(const std::string &name)
C_FLOAT64 sumOfSquares(const size_t &index, C_FLOAT64 *&residuals) const
CCopasiParameter * getParameter(const std::string &name)
static const std::string WeightMethodName[]
Definition: CExperiment.h:101
std::string add(const std::string &prefix, CCopasiObject *pObject)
C_FLOAT64 sumOfSquaresStore(const size_t &index, C_FLOAT64 *&dependentValues)
std::map< CCopasiObject *, size_t > mDependentObjects
Definition: CExperiment.h:680
size_t mExtendedTimeSeriesSize
Definition: CExperiment.h:691
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
const C_FLOAT64 & getTime() const
Definition: CModel.cpp:1190
const bool & isRowOriented() const
size_t getColumnValidValueCount(CCopasiObject *const &pObject) const
size_t size() const
Definition: CVector.h:100
long int flag
Definition: f2c.h:52
size_t extendedTimeSeriesSize() const
bool * mpRowOriented
Definition: CExperiment.h:595
const CCopasiTask::Type & getExperimentType() const
virtual std::string getObjectDisplayName(bool regular=true, bool richtext=false) const
bool calculateWeights()
const CVector< C_FLOAT64 > & getTimeData() const
static CKeyFactory * getKeyFactory()
unsigned C_INT32 * mpHeaderRow
Definition: CExperiment.h:600
CVector< C_FLOAT64 > mDataTime
Definition: CExperiment.h:620
CVector< C_FLOAT64 > mColumnObjectiveValue
Definition: CExperiment.h:673
C_FLOAT64 mIndependentValue
Definition: CExperiment.h:55
void initializeScalingMatrix()
#define C_FLOAT64
Definition: copasi.h:92
void updateFittedPoints()
void updateFittedPointValuesFromExtendedTimeSeries(const size_t &index)
CType * array()
Definition: CVector.h:139
bool applyPivot(const CVector< size_t > &pivot)
Definition: CMatrix.h:351
void setValues(const C_FLOAT64 &independent, const C_FLOAT64 &measured, const C_FLOAT64 &fitted, const C_FLOAT64 &weightedError)
std::string StringPrint(const char *format,...)
Definition: utility.cpp:87
bool compile(const std::vector< CCopasiContainer * > listOfContainer=CCopasiContainer::EmptyList)
WeightMethod * mpWeightMethod
Definition: CExperiment.h:590
const std::string & getSeparator() const
virtual size_t size() const
Definition: CMatrix.h:132
static std::vector< Refresh * > buildUpdateSequence(const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())
bool isValueDbl() const
bool setNumColumns(const unsigned C_INT32 &cols)
virtual void clear()
const CCopasiParameter::Value & getValue(const std::string &name) const
bool * mpNormalizeWeightsPerExperiment
Definition: CExperiment.h:580
bool calculateStatistics()
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
size_t guessColumnNumber(std::istream &is, const bool &rewind)
Definition: CTableCell.cpp:163
Definition: CModel.h:50
std::string getObjectCN(const size_t &index) const
bool setFirstRow(const unsigned C_INT32 &firstRow)
CCopasiParameterGroup * getGroup(const std::string &name)
virtual const CObjectInterface * getObject(const CCopasiObjectName &cn) const
static CLocaleString fromUtf8(const std::string &utf8)
bool setNumCols(const size_t &numCols)
C_FLOAT64 mRMS
Definition: CExperiment.h:668
CVector< C_FLOAT64 > mExtendedTimeSeries
Definition: CExperiment.h:685
static bool makePathAbsolute(std::string &relativePath, const std::string &absoluteTo)
Definition: CDirEntry.cpp:481
virtual void printResult(std::ostream *ostream) const
virtual size_t numCols() const
Definition: CMatrix.h:144
const unsigned C_INT32 & getNumColumns() const
const unsigned C_INT32 & getLastRow() const
bool updateModelWithIndependentData(const size_t &index)
CVector< C_FLOAT64 > mIndependentValues
Definition: CExperiment.h:656
bool setLastRow(const unsigned C_INT32 &lastRow)
CVector< UpdateMethod * > mIndependentUpdateMethods
Definition: CExperiment.h:650
bool operator==(const CExperiment &lhs, const CExperiment &rhs)
CCopasiObject * ObjectFromName(const std::vector< CCopasiContainer * > &listOfContainer, const CCopasiObjectName &CN) const
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
bool getNormalizeWeightsPerExperiment() const
bool setObjectCN(const size_t &index, const std::string &objectCN)
bool readColumnNames()
virtual CType * array()
Definition: CMatrix.h:337
CExperiment & operator=(const CExperiment &rhs)
static const std::string TypeName[]
Definition: CExperiment.h:79
C_FLOAT64 mMeasuredValue
Definition: CExperiment.h:56
bool setWeightMethod(const WeightMethod &weightMethod)
CMatrix< C_FLOAT64 > mScale
Definition: CExperiment.h:635
bool applyPivot(const CVectorCore< size_t > &pivot)
Definition: CVector.h:172
size_t getLastNotIgnoredColumn() const
bool read(std::istream &in, size_t &currentLine)
void initializeParameter()
void fixBuild55()
#define max(a, b)
Definition: f2c.h:176