COPASI API  4.16.103
CFitProblem.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 <cmath>
16 
17 #include "copasi.h"
18 
19 #include "CFitProblem.h"
20 #include "CFitItem.h"
21 #include "CFitTask.h"
22 #include "CExperimentSet.h"
23 #include "CExperiment.h"
24 
27 #include "model/CModel.h"
28 #include "model/CState.h"
30 #include "report/CKeyFactory.h"
37 
38 #include "lapack/blaswrap.h" //use blas
39 #include "lapack/lapackwrap.h" //use CLAPACK
40 
41 // Default constructor
43  const CCopasiContainer * pParent):
44  COptProblem(type, pParent),
45  mpParmSteadyStateCN(NULL),
46  mpParmTimeCourseCN(NULL),
47  mpExperimentSet(NULL),
48  mpSteadyState(NULL),
49  mpTrajectory(NULL),
50  mExperimentUpdateMethods(0, 0),
51  mExperimentUndoMethods(0, 0),
52  mExperimentConstraints(0, 0),
53  mExperimentDependentValues(0),
54  mpCrossValidationSet(NULL),
55  mCrossValidationUpdateMethods(0, 0),
56  mCrossValidationConstraints(0, 0),
57  mCrossValidationDependentValues(0),
58  mCrossValidationSolutionValue(mWorstValue),
59  mCrossValidationRMS(std::numeric_limits<C_FLOAT64>::quiet_NaN()),
60  mCrossValidationSD(std::numeric_limits<C_FLOAT64>::quiet_NaN()),
61  mCrossValidationObjective(mWorstValue),
62  mThresholdCounter(0),
63  mpTrajectoryProblem(NULL),
64  mpInitialState(NULL),
65  mResiduals(0),
66  mRMS(std::numeric_limits<C_FLOAT64>::quiet_NaN()),
67  mSD(std::numeric_limits<C_FLOAT64>::quiet_NaN()),
68  mParameterSD(0),
69  mFisher(0, 0),
70  mpFisherMatrixInterface(NULL),
71  mpFisherMatrix(NULL),
72  mFisherEigenvalues(0, 0),
73  mpFisherEigenvaluesMatrixInterface(NULL),
74  mpFisherEigenvaluesMatrix(NULL),
75  mFisherEigenvectors(0, 0),
76  mpFisherEigenvectorsMatrixInterface(NULL),
77  mpFisherEigenvectorsMatrix(NULL),
78  mFisherScaled(0, 0),
79  mpFisherScaledMatrixInterface(NULL),
80  mpFisherScaledMatrix(NULL),
81  mFisherScaledEigenvalues(0, 0),
82  mpFisherScaledEigenvaluesMatrixInterface(NULL),
83  mpFisherScaledEigenvaluesMatrix(NULL),
84  mFisherScaledEigenvectors(0, 0),
85  mpFisherScaledEigenvectorsMatrixInterface(NULL),
86  mpFisherScaledEigenvectorsMatrix(NULL),
87  mCorrelation(0, 0),
88  mpCorrelationMatrixInterface(NULL),
89  mpCorrelationMatrix(NULL),
90  mpCreateParameterSets(NULL),
91  mTrajectoryUpdate(false)
92 
93 {
94  initObjects();
96 }
97 
98 // copy constructor
100  const CCopasiContainer * pParent):
101  COptProblem(src, pParent),
102  mpParmSteadyStateCN(NULL),
103  mpParmTimeCourseCN(NULL),
104  mpExperimentSet(NULL),
105  mpSteadyState(NULL),
106  mpTrajectory(NULL),
107  mExperimentUpdateMethods(0, 0),
108  mExperimentUndoMethods(0, 0),
109  mExperimentConstraints(0, 0),
110  mExperimentDependentValues(src.mExperimentDependentValues),
111  mpCrossValidationSet(NULL),
112  mCrossValidationUpdateMethods(0, 0),
113  mCrossValidationConstraints(0, 0),
114  mCrossValidationDependentValues(src.mCrossValidationDependentValues),
115  mCrossValidationSolutionValue(mWorstValue),
116  mCrossValidationRMS(std::numeric_limits<C_FLOAT64>::quiet_NaN()),
117  mCrossValidationSD(std::numeric_limits<C_FLOAT64>::quiet_NaN()),
118  mCrossValidationObjective(mWorstValue),
119  mThresholdCounter(0),
120  mpTrajectoryProblem(NULL),
121  mpInitialState(NULL),
122  mResiduals(src.mResiduals),
123  mRMS(src.mRMS),
124  mSD(src.mSD),
125  mParameterSD(src.mParameterSD),
126  mFisher(src.mFisher),
127  mpFisherMatrixInterface(NULL),
128  mpFisherMatrix(NULL),
129  mFisherEigenvalues(src.mFisherEigenvalues),
130  mpFisherEigenvaluesMatrixInterface(NULL),
131  mpFisherEigenvaluesMatrix(NULL),
132  mFisherEigenvectors(src.mFisherEigenvectors),
133  mpFisherEigenvectorsMatrixInterface(NULL),
134  mpFisherEigenvectorsMatrix(NULL),
135  mFisherScaled(src.mFisherScaled),
136  mpFisherScaledMatrixInterface(NULL),
137  mpFisherScaledMatrix(NULL),
138  mFisherScaledEigenvalues(src.mFisherScaledEigenvalues),
139  mpFisherScaledEigenvaluesMatrixInterface(NULL),
140  mpFisherScaledEigenvaluesMatrix(NULL),
141  mFisherScaledEigenvectors(src.mFisherScaledEigenvectors),
142  mpFisherScaledEigenvectorsMatrixInterface(NULL),
143  mpFisherScaledEigenvectorsMatrix(NULL),
144  mCorrelation(src.mCorrelation),
145  mpCorrelationMatrixInterface(NULL),
146  mpCorrelationMatrix(NULL),
147  mpCreateParameterSets(NULL),
148  mTrajectoryUpdate(false)
149 {
150  initObjects();
152 }
153 
154 // Destructor
156 {
173 }
174 
176 {
179 
181  mpFisherMatrix = new CArrayAnnotation("Fisher Information Matrix", this, mpFisherMatrixInterface, false);
182  mpFisherMatrix->setDescription("Fisher Information Matrix");
183  mpFisherMatrix->setDimensionDescription(0, "Parameters");
184  mpFisherMatrix->setDimensionDescription(1, "Parameters");
185  mpFisherMatrix->setMode(CArrayAnnotation::STRINGS);
186 
189  mpFisherEigenvaluesMatrix->setDescription("FIM Eigenvalues");
190  mpFisherEigenvaluesMatrix->setDimensionDescription(0, "Eigenvalues");
191  mpFisherEigenvaluesMatrix->setDimensionDescription(1, "Result");
192  mpFisherEigenvaluesMatrix->setMode(CArrayAnnotation::NUMBERS);
193 
196  mpFisherEigenvectorsMatrix->setDescription("FIM Eigenvectors");
197  mpFisherEigenvectorsMatrix->setDimensionDescription(0, "Eigenvectors");
198  mpFisherEigenvectorsMatrix->setDimensionDescription(1, "Parameters");
199  mpFisherEigenvectorsMatrix->setMode(0, CArrayAnnotation::NUMBERS);
200  mpFisherEigenvectorsMatrix->setMode(1, CArrayAnnotation::STRINGS);
201 
203  mpFisherScaledMatrix = new CArrayAnnotation("Fisher Information Matrix (scaled)", this, mpFisherScaledMatrixInterface, false);
204  mpFisherScaledMatrix->setDescription("Fisher Information Matrix (scaled)");
205  mpFisherScaledMatrix->setDimensionDescription(0, "Parameters");
206  mpFisherScaledMatrix->setDimensionDescription(1, "Parameters");
207  mpFisherScaledMatrix->setMode(CArrayAnnotation::STRINGS);
208 
210  mpFisherScaledEigenvaluesMatrix = new CArrayAnnotation("FIM Eigenvalues (scaled)", this, mpFisherScaledEigenvaluesMatrixInterface, false);
211  mpFisherScaledEigenvaluesMatrix->setDescription("FIM Eigenvalues (scaled)");
212  mpFisherScaledEigenvaluesMatrix->setDimensionDescription(0, "Eigenvalues");
213  mpFisherScaledEigenvaluesMatrix->setDimensionDescription(1, "Result");
214  mpFisherScaledEigenvaluesMatrix->setMode(CArrayAnnotation::NUMBERS);
215 
218  mpFisherScaledEigenvectorsMatrix->setDescription("FIM Eigenvectors (scaled)");
219  mpFisherScaledEigenvectorsMatrix->setDimensionDescription(0, "Eigenvectors");
220  mpFisherScaledEigenvectorsMatrix->setDimensionDescription(1, "Parameters");
221  mpFisherScaledEigenvectorsMatrix->setMode(0, CArrayAnnotation::NUMBERS);
222  mpFisherScaledEigenvectorsMatrix->setMode(1, CArrayAnnotation::STRINGS);
223 
225  mpCorrelationMatrix = new CArrayAnnotation("Correlation Matrix", this, mpCorrelationMatrixInterface, false);
226  mpCorrelationMatrix->setDescription("Correlation Matrix");
227  mpCorrelationMatrix->setDimensionDescription(0, "Parameters");
228  mpCorrelationMatrix->setDimensionDescription(1, "Parameters");
229  mpCorrelationMatrix->setMode(CArrayAnnotation::STRINGS);
230 }
231 
233 {
234  removeParameter("Subtask");
235  mpParmSubtaskCN = NULL;
236  removeParameter("ObjectiveExpression");
238  *mpParmMaximize = false;
239 
244 
246  assertParameter("Create Parameter Sets", CCopasiParameter::BOOL, false)-> getValue().pBOOL;
247 
248  assertGroup("Experiment Set");
249 
250  assertGroup("Validation Set");
251 
252  elevateChildren();
253 }
254 
255 void CFitProblem::setCreateParameterSets(const bool & create)
256 {*mpCreateParameterSets = create;}
257 
259 {return *mpCreateParameterSets;}
260 
262 {
263  // This call is necessary since CFitProblem is derived from COptProblem.
264  if (!COptProblem::elevateChildren()) return false;
265 
266  // Due to a naming conflict the following parameters may have been overwritten during
267  // the load of a CopasiML file we replace them with default values if that was the case.
272 
273  CCopasiVectorN< CCopasiTask > * pTasks = NULL;
274  CCopasiDataModel* pDataModel = getObjectDataModel();
275 
276  if (pDataModel)
277  pTasks = pDataModel->getTaskList();
278 
279  if (pTasks == NULL)
280  pTasks = dynamic_cast<CCopasiVectorN< CCopasiTask > *>(getObjectAncestor("Vector"));
281 
282  if (pTasks)
283  {
284  size_t i, imax = pTasks->size();
285 
286  if (!mpParmSteadyStateCN->compare(0, 5 , "Task_") ||
287  *mpParmSteadyStateCN == "")
288  for (i = 0; i < imax; i++)
289  if ((*pTasks)[i]->getType() == CCopasiTask::steadyState)
290  {
291  *mpParmSteadyStateCN = (*pTasks)[i]->getCN();
292  break;
293  }
294 
295  if (!mpParmTimeCourseCN->compare(0, 5 , "Task_") ||
296  *mpParmTimeCourseCN == "")
297  for (i = 0; i < imax; i++)
298  if ((*pTasks)[i]->getType() == CCopasiTask::timeCourse)
299  {
300  *mpParmTimeCourseCN = (*pTasks)[i]->getCN();
301  break;
302  }
303  }
304 
305  std::map<std::string, std::string> ExperimentMap;
306  CCopasiParameterGroup * pGroup;
307  CExperiment * pExperiment;
308 
309  std::vector<CCopasiParameter *> * pExperiments =
310  getGroup("Experiment Set")->CCopasiParameter::getValue().pGROUP;
311  std::vector<CCopasiParameter *>::iterator itExp;
312  std::vector<CCopasiParameter *>::iterator endExp;
313 
314  for (itExp = pExperiments->begin(), endExp = pExperiments->end(); itExp != endExp; ++itExp)
315  if ((pGroup = dynamic_cast< CCopasiParameterGroup * >(*itExp)) != NULL &&
316  pGroup->getParameter("Key") != NULL)
317  ExperimentMap[*pGroup->getValue("Key").pKEY] = (*itExp)->getObjectName();
318 
320  elevate<CExperimentSet, CCopasiParameterGroup>(getGroup("Experiment Set"));
321 
322  if (!mpExperimentSet) return false;
323 
324  std::map<std::string, std::string>::iterator itMap;
325  std::map<std::string, std::string>::iterator endMap;
326 
327  for (itMap = ExperimentMap.begin(), endMap = ExperimentMap.end(); itMap != endMap; ++itMap)
328  {
329  pExperiment = mpExperimentSet->getExperiment(itMap->second);
330  itMap->second = pExperiment->CCopasiParameter::getKey();
331  pExperiment->setValue("Key", itMap->second);
332  }
333 
334  std::map<std::string, std::string> CrossValidationMap;
335 
336  // We have old CopasiML files (only snapshots and test releases), which use
337  // "Cross Validation Set"
338  pGroup = getGroup("Cross Validation Set");
339 
340  if (pGroup != NULL)
341  {
342  removeParameter("Validation Set");
343  pGroup->setObjectName("Validation Set");
344  }
345 
346  pExperiments =
347  getGroup("Validation Set")->CCopasiParameter::getValue().pGROUP;
348 
349  for (itExp = pExperiments->begin(), endExp = pExperiments->end(); itExp != endExp; ++itExp)
350  if ((pGroup = dynamic_cast< CCopasiParameterGroup * >(*itExp)) != NULL &&
351  pGroup->getParameter("Key") != NULL)
352  CrossValidationMap[*pGroup->getValue("Key").pKEY] = (*itExp)->getObjectName();
353 
354  // This intermediate elevation step should be not needed but Viusal C fails when directly
355  // going to the CCrossValidationSet.
356  elevate< CExperimentSet, CCopasiParameterGroup >(getGroup("Validation Set"));
358  elevate< CCrossValidationSet, CCopasiParameterGroup >(getGroup("Validation Set"));
359 
360  if (!mpCrossValidationSet) return false;
361 
362  for (itMap = CrossValidationMap.begin(), endMap = CrossValidationMap.end(); itMap != endMap; ++itMap)
363  {
364  pExperiment = mpCrossValidationSet->getExperiment(itMap->second);
365  itMap->second = pExperiment->CCopasiParameter::getKey();
366  pExperiment->setValue("Key", itMap->second);
367  }
368 
369  std::vector<COptItem * >::iterator it = mpOptItems->begin();
370  std::vector<COptItem * >::iterator end = mpOptItems->end();
371 
372  for (; it != end; ++it)
373  {
374  if (!((*it) = elevate<CFitItem, COptItem>(*it)))
375  return false;
376 
377  pExperiments =
378  (*it)->getParameter("Affected Experiments")->getValue().pGROUP;
379 
380  for (itExp = pExperiments->begin(), endExp = pExperiments->end(); itExp != endExp; ++itExp)
381  (*itExp)->setValue(ExperimentMap[*(*itExp)->getValue().pKEY]);
382 
383  pExperiments =
384  (*it)->getParameter("Affected Cross Validation Experiments")->getValue().pGROUP;
385 
386  for (itExp = pExperiments->begin(), endExp = pExperiments->end(); itExp != endExp; ++itExp)
387  (*itExp)->setValue(CrossValidationMap[*(*itExp)->getValue().pKEY]);
388  }
389 
390  it = mpConstraintItems->begin();
391  end = mpConstraintItems->end();
392 
393  for (; it != end; ++it)
394  {
395  if (!((*it) = elevate<CFitConstraint, COptItem>(*it)))
396  return false;
397 
398  pExperiments =
399  (*it)->getParameter("Affected Experiments")->getValue().pGROUP;
400 
401  for (itExp = pExperiments->begin(), endExp = pExperiments->end(); itExp != endExp; ++itExp)
402  (*itExp)->setValue(ExperimentMap[*(*itExp)->getValue().pKEY]);
403 
404  pExperiments =
405  (*it)->getParameter("Affected Cross Validation Experiments")->getValue().pGROUP;
406 
407  for (itExp = pExperiments->begin(), endExp = pExperiments->end(); itExp != endExp; ++itExp)
408  (*itExp)->setValue(CrossValidationMap[*(*itExp)->getValue().pKEY]);
409  }
410 
411  return true;
412 }
413 
415 {return COptProblem::setModel(pModel);}
416 
418 {return COptProblem::setCallBack(pCallBack);}
419 
421 {
422  bool success = true;
423 
424  mHaveStatistics = false;
425  mStoreResults = false;
426 
428  {
429  while (CCopasiMessage::peekLastMessage().getNumber() == MCOptimization + 5 ||
430  CCopasiMessage::peekLastMessage().getNumber() == MCOptimization + 7)
432 
435 
436  success = false;
437  }
438 
439  std::vector< CCopasiContainer * > ContainerList;
440  ContainerList.push_back(getObjectAncestor("Vector"));
441  CCopasiDataModel* pDataModel = getObjectDataModel();
442  assert(pDataModel != NULL);
443 
444  // We only need to initialize the steady-state task if steady-state data is present.
446  {
447  mpSteadyState =
448  dynamic_cast< CSteadyStateTask * >(pDataModel->ObjectFromName(ContainerList, *mpParmSteadyStateCN));
449 
450  if (mpSteadyState == NULL)
451  {
452  mpSteadyState =
453  static_cast<CSteadyStateTask *>((*pDataModel->getTaskList())["Steady-State"]);
454  }
455 
456  if (mpSteadyState == NULL) fatalError();
457 
460  }
461  else
462  {
463  mpSteadyState = NULL;
464  }
465 
467 
468  // We only need to initialize the trajectory task if time course data is present.
470  {
471  mpTrajectory =
472  dynamic_cast< CTrajectoryTask * >(pDataModel->ObjectFromName(ContainerList, *mpParmTimeCourseCN));
473 
474  if (mpTrajectory == NULL)
475  {
476  mpTrajectory =
477  static_cast<CTrajectoryTask *>((*pDataModel->getTaskList())["Time-Course"]);
478  }
479 
480  if (mpTrajectory == NULL) fatalError();
481 
483 
484  // do not update initial values when running fit
487 
489 
491  new CTrajectoryProblem(*static_cast<CTrajectoryProblem *>(mpTrajectory->getProblem()));
492  }
493  else
494  {
495  mpTrajectory = NULL;
496  }
497 
498  ContainerList.clear();
499 
500  ContainerList.push_back(mpModel);
501 
502  CFitTask * pTask = dynamic_cast<CFitTask *>(getObjectParent());
503 
504  if (pTask)
505  {
506  ContainerList.push_back(pTask);
507 
508  if (mpSteadyState != NULL)
509  {
510  ContainerList.push_back(mpSteadyState);
511  }
512 
513  if (mpTrajectory != NULL)
514  {
515  ContainerList.push_back(mpTrajectory);
516  }
517  }
518 
519  if (!mpExperimentSet->compile(ContainerList)) return false;
520 
521  // Initialize the object set which the experiment independent objects.
522  std::vector< std::set< const CCopasiObject * > > ObjectSet;
523  ObjectSet.resize(mpExperimentSet->getExperimentCount());
524  size_t i, imax;
525 
526  for (i = 0, imax = mpExperimentSet->getExperimentCount(); i < imax; i++)
527  {
529  }
530 
531  // Build a matrix of experiment and experiment local items.
533  mpOptItems->size());
535 
536  // Build a matrix of experiment and experiment local undo items.
538  mpOptItems->size());
539  mExperimentUndoMethods = NULL;
540 
542 
543  std::vector<COptItem * >::iterator it = mpOptItems->begin();
544  std::vector<COptItem * >::iterator end = mpOptItems->end();
545 
546  std::vector<COptItem * >::iterator itTmp;
547 
548  CFitItem * pItem;
549  size_t j;
550  size_t Index;
551 
552  imax = mSolutionVariables.size();
553 
554  mFisher.resize(imax, imax);
556  mFisherEigenvalues.resize(imax, 1);
558  mFisherEigenvectors.resize(imax, imax);
560  mFisherScaled.resize(imax, imax);
564  mFisherScaledEigenvectors.resize(imax, imax);
566  mCorrelation.resize(imax, imax);
568 
569  for (j = 0; it != end; ++it, j++)
570  {
571  pItem = static_cast<CFitItem *>(*it);
572  pItem->updateBounds(mpOptItems->begin());
573 
574  std::string Annotation = pItem->getObjectDisplayName();
575 
576  imax = pItem->getExperimentCount();
577 
578  if (imax == 0)
579  {
580  for (i = 0, imax = mpExperimentSet->getExperimentCount(); i < imax; i++)
581  {
582  mExperimentUpdateMethods(i, j) = pItem->COptItem::getUpdateMethod();
583  ObjectSet[i].insert(pItem->getObject());
584  }
585  }
586  else
587  {
588  Annotation += "; {" + pItem->getExperiments() + "}";
589 
590  for (i = 0; i < imax; i++)
591  {
592  if ((Index = mpExperimentSet->keyToIndex(pItem->getExperiment(i))) == C_INVALID_INDEX)
593  return false;
594 
595  mExperimentUpdateMethods(Index, j) = pItem->COptItem::getUpdateMethod();
596  ObjectSet[Index].insert(pItem->getObject());
597 
598  // We need to undo the changes for all non affected experiments.
599  // We can do that by adding the update method with the current model value to
600  mExperimentUndoMethods(Index, j) = pItem->COptItem::getUpdateMethod();
601  };
602  }
603 
604  mpFisherMatrix->setAnnotationString(0, j, Annotation);
605  mpFisherMatrix->setAnnotationString(1, j, Annotation);
607  mpFisherScaledMatrix->setAnnotationString(0, j, Annotation);
608  mpFisherScaledMatrix->setAnnotationString(1, j, Annotation);
610  mpCorrelationMatrix->setAnnotationString(0, j, Annotation);
611  mpCorrelationMatrix->setAnnotationString(1, j, Annotation);
612  }
613 
614  // Create a joined sequence of update methods for parameters and independent values.
615  for (i = 0, imax = mpExperimentSet->getExperimentCount(); i < imax; i++)
616  {
617  ObjectSet[i].erase(NULL);
619  }
620 
621  // Build a matrix of experiment and constraint items;
623  mpConstraintItems->size());
624  mExperimentConstraints = NULL;
626  ObjectSet.clear();
627  ObjectSet.resize(mpExperimentSet->getExperimentCount());
628 
629  it = mpConstraintItems->begin();
630  end = mpConstraintItems->end();
631 
632  CFitConstraint * pConstraint;
633  std::set< const CCopasiObject * >::const_iterator itDepend;
634  std::set< const CCopasiObject * >::const_iterator endDepend;
635 
636  for (j = 0; it != end; ++it, j++)
637  {
638  pConstraint = static_cast<CFitConstraint *>(*it);
639  itDepend = pConstraint->getDirectDependencies().begin();
640  endDepend = pConstraint->getDirectDependencies().end();
641  imax = pConstraint->getExperimentCount();
642 
643  if (imax == 0)
644  {
645  for (i = 0, imax = mpExperimentSet->getExperimentCount(); i < imax; i++)
646  {
647  mExperimentConstraints(i, j) = pConstraint;
648  ObjectSet[i].insert(itDepend, endDepend);
649  }
650  }
651  else
652  {
653  for (i = 0; i < imax; i++)
654  {
655  if ((Index = mpExperimentSet->keyToIndex(pConstraint->getExperiment(i))) == C_INVALID_INDEX)
656  return false;
657 
658  mExperimentConstraints(Index, j) = pConstraint;
659  ObjectSet[Index].insert(itDepend, endDepend);
660  };
661  }
662  }
663 
664  for (i = 0, imax = mpExperimentSet->getExperimentCount(); i < imax; i++)
666 
668 
669  if (!mpCrossValidationSet->compile(ContainerList)) return false;
670 
671  // Initialize the object set which the experiment independent objects.
672  ObjectSet.clear();
673  ObjectSet.resize(mpCrossValidationSet->getExperimentCount());
674 
675  for (i = 0, imax = mpCrossValidationSet->getExperimentCount(); i < imax; i++)
676  {
678  }
679 
680  // Build a matrix of cross validation experiments and local items.
682  mpOptItems->size());
684 
685  // Build a matrix of experiment and experiment local undo items.
687  mpOptItems->size());
689 
691 
692  it = mpOptItems->begin();
693  end = mpOptItems->end();
694 
695  for (j = 0; it != end; ++it, j++)
696  {
697  pItem = static_cast<CFitItem *>(*it);
698  pItem->updateBounds(mpOptItems->begin());
699 
700  imax = pItem->getCrossValidationCount();
701 
702  if (imax == 0)
703  {
704  for (i = 0, imax = mpCrossValidationSet->getExperimentCount(); i < imax; i++)
705  {
706  mCrossValidationUpdateMethods(i, j) = pItem->COptItem::getUpdateMethod();
707  ObjectSet[i].insert(pItem->getObject());
708  }
709  }
710  else
711  {
712  for (i = 0; i < imax; i++)
713  {
715  return false;
716 
717  mCrossValidationUpdateMethods(Index, j) = pItem->COptItem::getUpdateMethod();
718  ObjectSet[Index].insert(pItem->getObject());
719 
720  // We need to undo the changes for all non affected cross validations.
721  // We can do that by adding the update method with the current model value to
722  mCrossValidationUndoMethods(Index, j) = pItem->COptItem::getUpdateMethod();
723  };
724  }
725  }
726 
727  // Create a joined sequence of update methods for parameters and independent values.
728  for (i = 0, imax = mpCrossValidationSet->getExperimentCount(); i < imax; i++)
729  {
731  }
732 
733  // Build a matrix of cross validation experiments and constraint items;
735  mpConstraintItems->size());
738  ObjectSet.clear();
739  ObjectSet.resize(mpCrossValidationSet->getExperimentCount());
740 
741  it = mpConstraintItems->begin();
742  end = mpConstraintItems->end();
743 
744  for (j = 0; it != end; ++it, j++)
745  {
746  pConstraint = static_cast<CFitConstraint *>(*it);
747 
748  imax = pConstraint->getCrossValidationCount();
749 
750  if (imax == 0)
751  {
752  for (i = 0, imax = mpCrossValidationSet->getExperimentCount(); i < imax; i++)
753  {
754  mCrossValidationConstraints(i, j) = pConstraint;
755  ObjectSet[i].insert(pConstraint->getObject());
756  }
757  }
758  else
759  {
760  for (i = 0; i < imax; i++)
761  {
762  if ((Index = mpCrossValidationSet->keyToIndex(pConstraint->getCrossValidation(i))) == C_INVALID_INDEX)
763  return false;
764 
765  mCrossValidationConstraints(Index, j) = pConstraint;
766  ObjectSet[Index].insert(pConstraint->getObject());
767  };
768  }
769  }
770 
771  for (i = 0, imax = mpCrossValidationSet->getExperimentCount(); i < imax; i++)
773 
775 
777  mThresholdCounter = 0;
778 
779  setResidualsRequired(false);
780 
783 
784  return true;
785 }
786 
788 {
789  std::vector< COptItem * >::const_iterator it = mpConstraintItems->begin();
790  std::vector< COptItem * >::const_iterator end = mpConstraintItems->end();
791 
793 
794  for (; it != end; ++it)
795  if (static_cast<CFitConstraint *>(*it)->getConstraintViolation() > 0.0)
796  {
798  return false;
799  }
800 
801  return true;
802 }
803 
804 /**
805  * Utility function creating a parameter set for each experiment
806  */
808 {
809  if (pExp == NULL) return;
810 
811  CModel* model = pExp->getObjectDataModel()->getModel();
812  std::string origname = "PE: " + UTCTimeStamp() + " Exp: " + pExp->getObjectName();
813  std::string name = origname;
814  int count = 0;
815 
816  while (model->getModelParameterSets().getIndex(name) != C_INVALID_INDEX)
817  {
818  std::stringstream str; str << origname << " (" << ++count << ")";
819  name = str.str();
820  }
821 
822  CModelParameterSet* set = new CModelParameterSet(name);
823  model->getModelParameterSets().add(set, true);
824  set->createFromModel();
825 }
826 
828 {
829  mCounter += 1;
830  bool Continue = true;
831 
832  size_t i, imax = mpExperimentSet->getExperimentCount();
833  size_t j;
834  size_t kmax;
835  mCalculateValue = 0.0;
836 
837  CExperiment * pExp = NULL;
838 
839  C_FLOAT64 * Residuals = mResiduals.array();
840  C_FLOAT64 * DependentValues = mExperimentDependentValues.array();
843  std::vector<COptItem *>::iterator itItem;
844  std::vector<COptItem *>::iterator endItem = mpOptItems->end();
845  std::vector<COptItem *>::iterator itConstraint;
846  std::vector<COptItem *>::iterator endConstraint = mpConstraintItems->end();
847 
848  std::vector< Refresh *>::const_iterator itRefresh;
849  std::vector< Refresh *>::const_iterator endRefresh;
850 
851  // Reset the constraints memory
852  for (itConstraint = mpConstraintItems->begin(); itConstraint != endConstraint; ++itConstraint)
853  static_cast<CFitConstraint *>(*itConstraint)->resetConstraintViolation();
854 
855  CFitConstraint **ppConstraint = mExperimentConstraints.array();
856  CFitConstraint **ppConstraintEnd;
857 
858  try
859  {
860  for (i = 0; i < imax && Continue; i++) // For each experiment
861  {
862  pExp = mpExperimentSet->getExperiment(i);
863 
864  // Set the model to its original state.
865  // TODO CRITICAL This is the incorrect state if we are running as part of a scan.
868 
869  // set the global and experiment local fit item values.
870  for (itItem = mpOptItems->begin(); itItem != endItem; itItem++, pUpdate++)
871  if (*pUpdate)
872  (**pUpdate)(static_cast<CFitItem *>(*itItem)->getLocalValue());
873 
874  kmax = pExp->getNumDataRows();
875 
876  switch (pExp->getExperimentType())
877  {
879 
880  // set independent data
881  for (j = 0; j < kmax && Continue; j++) // For each data row;
882  {
884 
885  // We need to apply the parameter and independent
886  // value updates as one unit.
887  itRefresh = mExperimentInitialRefreshes[i].begin();
888  endRefresh = mExperimentInitialRefreshes[i].end();
889 
890  while (itRefresh != endRefresh)
891  (**itRefresh++)();
892 
893  Continue = mpSteadyState->process(true);
894 
895  if (!Continue)
896  {
897  mFailedCounter++;
899  break;
900  }
901 
902  // We check after each simulation whether the constraints are violated.
903  // Make sure the constraint values are up to date.
904  itRefresh = mExperimentConstraintRefreshes[i].begin();
905  endRefresh = mExperimentConstraintRefreshes[i].end();
906 
907  for (; itRefresh != endRefresh; ++itRefresh)
908  (**itRefresh)();
909 
910  ppConstraint = mExperimentConstraints[i];
911  ppConstraintEnd = ppConstraint + mExperimentConstraints.numCols();
912 
913  for (; ppConstraint != ppConstraintEnd; ++ppConstraint)
914  if (*ppConstraint)(*ppConstraint)->calculateConstraintViolation();
915 
916  if (mStoreResults)
917  mCalculateValue += pExp->sumOfSquaresStore(j, DependentValues);
918  else
919  mCalculateValue += pExp->sumOfSquares(j, Residuals);
920  }
921 
923  {
925  }
926 
927  break;
928 
930 
931  size_t numIntermediateSteps;
932 
933  if (mStoreResults)
934  {
935  //calculate a reasonable number of intermediate points
936  numIntermediateSteps = 4; //TODO
937  //resize the storage for the extended time series
938  pExp->initExtendedTimeSeries(numIntermediateSteps * kmax - numIntermediateSteps + 1);
939  }
940 
941  for (j = 0; j < kmax && Continue; j++) // For each data row;
942  {
943  if (j)
944  {
945  if (mStoreResults)
946  {
947  //do additional intermediate steps for nicer display
948  C_FLOAT64 ttt;
949  size_t ic;
950 
951  for (ic = 1; ic < numIntermediateSteps; ++ic)
952  {
953  ttt = pExp->getTimeData()[j - 1] + (pExp->getTimeData()[j] - pExp->getTimeData()[j - 1]) * (C_FLOAT64(ic) / numIntermediateSteps);
955  //save the simulation results in the experiment
956  pExp->storeExtendedTimeSeriesData(ttt);
957  }
958  }
959 
960  //do the regular step
961  mpTrajectory->processStep(pExp->getTimeData()[j]);
962  }
963  else
964  {
965  // Set independent data. A time course only has one set of
966  // independent data.
968 
969  // We need to apply the parameter and independent
970  // value updates as one unit.
971  itRefresh = mExperimentInitialRefreshes[i].begin();
972  endRefresh = mExperimentInitialRefreshes[i].end();
973 
974  while (itRefresh != endRefresh)
975  (**itRefresh++)();
976 
977  static_cast< CTrajectoryProblem * >(mpTrajectory->getProblem())->setStepNumber(1);
978  mpTrajectory->processStart(true);
979 
980  if (pExp->getTimeData()[0] != mpModel->getInitialTime())
981  {
982  mpTrajectory->processStep(pExp->getTimeData()[0]);
983  }
984  }
985 
986  // We check after each simulation step whether the constraints are violated.
987  // Make sure the constraint values are up to date.
988  itRefresh = mExperimentConstraintRefreshes[i].begin();
989  endRefresh = mExperimentConstraintRefreshes[i].end();
990 
991  for (; itRefresh != endRefresh; ++itRefresh)
992  (**itRefresh)();
993 
994  ppConstraint = mExperimentConstraints[i];
995  ppConstraintEnd = ppConstraint + mExperimentConstraints.numCols();
996 
997  for (; ppConstraint != ppConstraintEnd; ++ppConstraint)
998  if (*ppConstraint)(*ppConstraint)->calculateConstraintViolation();
999 
1000  if (mStoreResults)
1001  mCalculateValue += pExp->sumOfSquaresStore(j, DependentValues);
1002  else
1003  mCalculateValue += pExp->sumOfSquares(j, Residuals);
1004 
1005  if (mStoreResults)
1006  {
1007  //additionally also store the the simulation result for the extended time series
1008  pExp->storeExtendedTimeSeriesData(pExp->getTimeData()[j]);
1009  }
1010  }
1011 
1013  {
1015  }
1016 
1017  break;
1018 
1019  default:
1020  break;
1021  }
1022 
1023  // restore independent data
1025 
1026  // restore experiment local values
1027  const C_FLOAT64 *pOriginal = mOriginalVariables.array();
1028  const C_FLOAT64 *pOriginalEnd = pOriginal + mOriginalVariables.size();
1029  bool RefreshNeeded = false;
1030 
1031  // set the global and experiment local fit item values.
1032  for (; pOriginal != pOriginalEnd; pOriginal++, pUndo++)
1033  if (*pUndo)
1034  {
1035  (**pUndo)(*pOriginal);
1036  RefreshNeeded = true;
1037  }
1038 
1039  if (RefreshNeeded)
1040  {
1041  // Update initial values which changed due to the fit item values.
1042  itRefresh = mExperimentInitialRefreshes[i].begin();
1043  endRefresh = mExperimentInitialRefreshes[i].end();
1044 
1045  while (itRefresh != endRefresh)
1046  (**itRefresh++)();
1047  }
1048  }
1049  }
1050 
1051  catch (CCopasiException &)
1052  {
1053  // We do not want to clog the message cue.
1055 
1056  mFailedCounter++;
1057 #ifdef XXX
1058  std::vector<COptItem * >::iterator it = mpOptItems->begin();
1059  std::vector<COptItem * >::iterator end = mpOptItems->end();
1060 
1061  for (; it != end; it++)
1062  std::cout << *(*it)->getObjectValue() << " ";
1063 
1064  std::cout << std::endl;
1065 #endif
1067 
1068  if (pExp)
1069  {
1071 
1072  // Update initial values which changed due to the fit item values.
1073  itRefresh = mExperimentInitialRefreshes[i].begin();
1074  endRefresh = mExperimentInitialRefreshes[i].end();
1075 
1076  while (itRefresh != endRefresh)
1077  (**itRefresh++)();
1078  }
1079  }
1080 
1081  catch (...)
1082  {
1083  mFailedCounter++;
1085 
1086  if (pExp)
1087  {
1089  // Update initial values which changed due to the fit item values.
1090  itRefresh = mExperimentInitialRefreshes[i].begin();
1091  endRefresh = mExperimentInitialRefreshes[i].end();
1092 
1093  while (itRefresh != endRefresh)
1094  (**itRefresh++)();
1095  }
1096  }
1097 
1098  if (isnan(mCalculateValue))
1100 
1102 
1103  return true;
1104 }
1105 
1106 bool CFitProblem::restore(const bool & updateModel)
1107 {
1108  bool success = true;
1109 
1110  if (mpTrajectory != NULL)
1111  {
1112  success &= mpTrajectory->restore();
1114  }
1115 
1116  if (mpSteadyState != NULL)
1117  success &= mpSteadyState->restore();
1118 
1119  if (mpTrajectoryProblem)
1121 
1122  success &= COptProblem::restore(updateModel);
1123 
1126 
1127  return success;
1128 }
1129 
1130 void CFitProblem::print(std::ostream * ostream) const
1131 {*ostream << *this;}
1132 
1133 void CFitProblem::printResult(std::ostream * ostream) const
1134 {
1135  std::ostream & os = *ostream;
1136 
1137  if (mSolutionVariables.size() == 0)
1138  {
1139  return;
1140  }
1141 
1142  os << "Objective Function Value:\t" << mSolutionValue << std::endl;
1143  os << "Standard Deviation:\t" << mSD << std::endl;
1144 
1145  CCopasiTimeVariable CPUTime = const_cast<CFitProblem *>(this)->mCPUTime.getElapsedTime();
1146 
1147  os << "Function Evaluations:\t" << mCounter << std::endl;
1148  os << "CPU Time [s]:\t"
1149  << CCopasiTimeVariable::LL2String(CPUTime.getSeconds(), 1) << "."
1150  << CCopasiTimeVariable::LL2String(CPUTime.getMilliSeconds(true), 3) << std::endl;
1151  os << "Evaluations/Second [1/s]:\t" << mCounter / (C_FLOAT64)(CPUTime.getMilliSeconds() / 1e3) << std::endl;
1152  os << std::endl;
1153 
1154  std::vector< COptItem * >::const_iterator itItem =
1155  mpOptItems->begin();
1156  std::vector< COptItem * >::const_iterator endItem =
1157  mpOptItems->end();
1158 
1159  CFitItem * pFitItem;
1160  CExperiment * pExperiment;
1161 
1162  size_t i, j;
1163 
1164  os << "\tParameter\tValue\tGradient\tStandard Deviation" << std::endl;
1165 
1166  for (i = 0; itItem != endItem; ++itItem, i++)
1167  {
1168  os << "\t" << (*itItem)->getObjectDisplayName();
1169  pFitItem = static_cast<CFitItem *>(*itItem);
1170 
1171  if (pFitItem->getExperimentCount() != 0)
1172  {
1173  os << " (";
1174 
1175  for (j = 0; j < pFitItem->getExperimentCount(); j++)
1176  {
1177  if (j) os << ", ";
1178 
1179  pExperiment =
1180  dynamic_cast< CExperiment * >(CCopasiRootContainer::getKeyFactory()->get(pFitItem->getExperiment(j)));
1181 
1182  if (pExperiment)
1183  os << pExperiment->getObjectName();
1184  }
1185 
1186  os << ")";
1187  }
1188 
1189  if (mHaveStatistics)
1190  {
1191  os << ":\t" << mSolutionVariables[i];
1192  os << "\t" << mGradient[i];
1193  os << "\t" << mParameterSD[i];
1194  }
1195  else
1196  {
1197  os << ":\t" << std::numeric_limits<C_FLOAT64>::quiet_NaN();
1198  os << "\t" << std::numeric_limits<C_FLOAT64>::quiet_NaN();
1199  os << "\t" << std::numeric_limits<C_FLOAT64>::quiet_NaN();
1200  }
1201 
1202  os << std::endl;
1203  }
1204 
1205  os << std::endl;
1206 
1207  size_t k, kmax = mpExperimentSet->getExperimentCount();
1208 
1209  for (k = 0; k < kmax; k++)
1210  {
1212  os << std::endl;
1213  }
1214 
1216  {
1217  os << "Fisher Information Matrix:" << std::endl;
1218  os << " " << mFisher << std::endl;
1219 
1220  os << "FIM Eigenvalues:" << std::endl;
1221  os << " " << mFisherEigenvalues << std::endl;
1222 
1223  os << "FIM Eigenvectors corresponding to Eigenvalues:" << std::endl;
1224  os << " " << mFisherEigenvectors << std::endl;
1225 
1226  os << "Fisher Information Matrix (scaled):" << std::endl;
1227  os << " " << mFisherScaled << std::endl;
1228 
1229  os << "FIM Eigenvalues (scaled):" << std::endl;
1230  os << " " << mFisherScaledEigenvalues << std::endl;
1231 
1232  os << "FIM Eigenvectors (scaled) corresponding to Eigenvalues:" << std::endl;
1233  os << " " << mFisherScaledEigenvectors << std::endl;
1234 
1235  os << "Correlation Matrix:" << std::endl;
1236  os << " " << mCorrelation << std::endl;
1237  }
1238 }
1239 
1240 std::ostream &operator<<(std::ostream &os, const CFitProblem & o)
1241 {
1242  os << "Problem Description:" << std::endl;
1243 
1244  os << "Subtask: " << std::endl;
1245 
1246  if (o.mpSteadyState)
1248 
1249  if (o.mpTrajectory)
1250  o.mpTrajectory->getDescription().print(&os);
1251 
1252  if (!o.mpTrajectory && !o.mpSteadyState)
1253  os << "No Subtask specified.";
1254 
1255  os << std::endl;
1256 
1257  os << "List of Fitting Items:" << std::endl;
1258 
1259  std::vector< COptItem * >::const_iterator itItem =
1260  o.mpOptItems->begin();
1261  std::vector< COptItem * >::const_iterator endItem =
1262  o.mpOptItems->end();
1263 
1264  for (; itItem != endItem; ++itItem)
1265  os << " " << *static_cast<CFitItem *>(*itItem) << std::endl;
1266 
1267  itItem = o.mpConstraintItems->begin();
1268  endItem = o.mpConstraintItems->end();
1269 
1270  for (; itItem != endItem; ++itItem)
1271  os << " " << *static_cast<CFitItem *>(*itItem) << std::endl;
1272 
1273  return os;
1274 }
1275 
1277 {
1279 }
1280 
1282 {return true;}
1283 
1284 bool CFitProblem::setResidualsRequired(const bool & required)
1285 {
1286  if (required)
1287  {
1289  size_t ii;
1290 
1291  for (ii = 0; ii < mResiduals.size(); ++ii)
1292  mResiduals[ii] = 0.0; //std::numeric_limits<C_FLOAT64>::quiet_NaN();
1293 
1294  //it would make sense to initialize it with NaN,
1295  //but some methods expect 0.0 as the residual for a missing data point
1296  }
1297  else
1298  mResiduals.resize(0);
1299 
1300  return true;
1301 }
1302 
1304 {
1305  return mResiduals;
1306 }
1307 
1309  const C_FLOAT64 & resolution)
1310 {
1311  // Set the current values to the solution values.
1312  size_t i, imax = mSolutionVariables.size();
1313  size_t l;
1314 
1315  mRMS = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1316  mSD = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1317 
1318  mCrossValidationRMS = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1319  mCrossValidationSD = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1320 
1321  mParameterSD.resize(imax);
1322  mParameterSD = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1323 
1324  mFisher = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1325  mFisherEigenvectors = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1326  mFisherEigenvalues = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1327  mFisherScaled = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1328  mFisherScaledEigenvectors = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1329  mFisherScaledEigenvalues = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1330  mGradient.resize(imax);
1331  mGradient = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1332 
1333  // Recalculate the best solution.
1334  for (i = 0; i < imax; i++)
1336 
1337  // For Output
1338  mStoreResults = true;
1339  calculate();
1340  mStoreResults = false;
1341 
1342  // The statistics need to be calculated for the result, i.e., now.
1344  size_t ValidDataCount = mpExperimentSet->getValidValueCount();
1345 
1346  if (ValidDataCount)
1347  mRMS = sqrt(mSolutionValue / ValidDataCount);
1348 
1349  if (ValidDataCount > imax)
1350  mSD = sqrt(mSolutionValue / (ValidDataCount - imax));
1351 
1353 
1355 
1356  size_t lmax = this->mCrossValidationDependentValues.size();
1357 
1358  if (lmax)
1360 
1361  if (lmax > imax)
1362  mCrossValidationSD = sqrt(mCrossValidationSolutionValue / (lmax - imax));
1363 
1364  mHaveStatistics = true;
1365 
1366  // Make sure the timer is accurate.
1367  (*mCPUTime.getRefresh())();
1368 
1369  if (mSolutionValue == mWorstValue)
1370  return false;
1371 
1373  {
1374  setResidualsRequired(true);
1375  size_t j, jmax = mResiduals.size();
1376  calculate();
1377 
1378  // Keep the results
1379  CVector< C_FLOAT64 > SolutionResiduals = mResiduals;
1380 
1381  CMatrix< C_FLOAT64 > DeltaResidualDeltaParameter;
1382 
1383  bool CalculateFIM = true;
1384 
1385  try
1386  {
1387  DeltaResidualDeltaParameter.resize(imax, jmax);
1388  }
1389 
1390  catch (CCopasiException & /*Exception*/)
1391  {
1392  CalculateFIM = false;
1393  }
1394 
1395  C_FLOAT64 * pDeltaResidualDeltaParameter = DeltaResidualDeltaParameter.array();
1396 
1397  C_FLOAT64 Current;
1398  C_FLOAT64 Delta;
1399 
1400  // Calculate the gradient
1401  for (i = 0; i < imax; i++)
1402  {
1403  Current = mSolutionVariables[i];
1404 
1405  if (fabs(Current) > resolution)
1406  {
1407  (*mUpdateMethods[i])(Current * (1.0 + factor));
1408  Delta = 1.0 / (Current * factor);
1409  }
1410  else
1411  {
1412  (*mUpdateMethods[i])(resolution);
1413  Delta = 1.0 / resolution;
1414  }
1415 
1416  calculate();
1417 
1418  mGradient[i] = (mCalculateValue - mSolutionValue) * Delta;
1419 
1420  C_FLOAT64 * pSolutionResidual = SolutionResiduals.array();
1421  C_FLOAT64 * pResidual = mResiduals.array();
1422 
1423  if (CalculateFIM)
1424  for (j = 0; j < jmax; j++, ++pDeltaResidualDeltaParameter, ++pSolutionResidual, ++pResidual)
1425  *pDeltaResidualDeltaParameter = (*pResidual - *pSolutionResidual) * Delta;
1426 
1427  // Restore the value
1428  (*mUpdateMethods[i])(Current);
1429  }
1430 
1431  if (!CalculateFIM)
1432  {
1433  // Make sure the timer is accurate.
1434  (*mCPUTime.getRefresh())();
1435 
1437  return false;
1438  }
1439 
1440  // Construct the fisher information matrix
1441  for (i = 0; i < imax; i++)
1442  for (l = 0; l <= i; l++)
1443  {
1444  C_FLOAT64 & tmp = mFisher(i, l);
1445 
1446  tmp = 0.0;
1447 
1448  C_FLOAT64 * pI = DeltaResidualDeltaParameter[i];
1449  C_FLOAT64 * pL = DeltaResidualDeltaParameter[l];
1450 
1451  for (j = 0; j < jmax; j++, ++pI, ++pL)
1452  tmp += *pI * *pL;
1453 
1454  tmp *= 2.0;
1455 
1456  // The Fisher matrix is symmetric.
1457  if (l != i)
1458  mFisher(l, i) = tmp;
1459  }
1460 
1461  // Construct the FIM Eigenvalue/-vector matrix
1462 
1463  /* int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a,
1464  integer *lda, doublereal *w, doublereal *work, integer *lwork,
1465  integer *info) */
1466 
1467  /* Purpose */
1468  /* ======= */
1469 
1470  /* DSYEV computes all eigenvalues and, optionally, eigenvectors of a */
1471  /* real symmetric matrix A. */
1472 
1473  /* Arguments */
1474  /* ========= */
1475 
1476  /* JOBZ (input) CHARACTER*1 */
1477  /* = 'N': Compute eigenvalues only; */
1478  /* = 'V': Compute eigenvalues and eigenvectors. */
1479 
1480  /* UPLO (input) CHARACTER*1 */
1481  /* = 'U': Upper triangle of A is stored; */
1482  /* = 'L': Lower triangle of A is stored. */
1483 
1484  /* N (input) INTEGER */
1485  /* The order of the matrix A. N >= 0. */
1486 
1487  /* A (input/output) DOUBLE PRECISION array, dimension (LDA, N) */
1488  /* On entry, the symmetric matrix A. If UPLO = 'U', the */
1489  /* leading N-by-N upper triangular part of A contains the */
1490  /* upper triangular part of the matrix A. If UPLO = 'L', */
1491  /* the leading N-by-N lower triangular part of A contains */
1492  /* the lower triangular part of the matrix A. */
1493  /* On exit, if JOBZ = 'V', then if INFO = 0, A contains the */
1494  /* orthonormal eigenvectors of the matrix A. */
1495  /* If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') */
1496  /* or the upper triangle (if UPLO='U') of A, including the */
1497  /* diagonal, is destroyed. */
1498 
1499  /* LDA (input) INTEGER */
1500  /* The leading dimension of the array A. LDA >= max(1,N). */
1501 
1502  /* W (output) DOUBLE PRECISION array, dimension (N) */
1503  /* If INFO = 0, the eigenvalues in ascending order. */
1504 
1505  /* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
1506  /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
1507 
1508  /* LWORK (input) INTEGER */
1509  /* The length of the array WORK. LWORK >= max(1,3*N-1). */
1510  /* For optimal efficiency, LWORK >= (NB+2)*N, */
1511  /* where NB is the blocksize for DSYTRD returned by ILAENV. */
1512 
1513  /* If LWORK = -1, then a workspace query is assumed; the routine */
1514  /* only calculates the optimal size of the WORK array, returns */
1515  /* this value as the first entry of the WORK array, and no error */
1516  /* message related to LWORK is issued by XERBLA. */
1517 
1518  /* INFO (output) INTEGER */
1519  /* = 0: successful exit */
1520  /* < 0: if INFO = -i, the i-th argument had an illegal value */
1521  /* > 0: if INFO = i, the algorithm failed to converge; i */
1522  /* off-diagonal elements of an intermediate tridiagonal */
1523  /* form did not converge to zero. */
1524 
1526 
1527  char JOBZ = 'V'; //also compute eigenvectors
1528  char UPLO = 'U'; //upper triangle
1529  C_INT N = (C_INT) imax;
1530  C_INT LDA = std::max< C_INT >(1, N);
1531  CVector<C_FLOAT64> WORK; //WORK space
1532  WORK.resize(1);
1533  C_INT LWORK = -1; //first query the memory need
1534  C_INT INFO = 0;
1535 
1536  dsyev_(&JOBZ,
1537  &UPLO,
1538  &N,
1540  &LDA,
1542  WORK.array(),
1543  &LWORK,
1544  &INFO);
1545 
1546  LWORK = (C_INT) WORK[0];
1547  WORK.resize(LWORK);
1548 
1549  //now do the real calculation
1550  dsyev_(&JOBZ,
1551  &UPLO,
1552  &N,
1554  &LDA,
1556  WORK.array(),
1557  &LWORK,
1558  &INFO);
1559 
1560  if (INFO != 0)
1561  {
1563 
1564  mFisherEigenvectors = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1565  mFisherEigenvalues = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1566  }
1567 
1568  for (i = 0; i < imax; i++)
1569  {
1570  for (j = 0; j < imax; j++)
1571  {
1573  }
1574  }
1575 
1577 
1578  //now do the real calculation
1579  dsyev_(&JOBZ,
1580  &UPLO,
1581  &N,
1583  &LDA,
1585  WORK.array(),
1586  &LWORK,
1587  &INFO);
1588 
1589  if (INFO != 0)
1590  {
1592 
1593  mFisherScaledEigenvectors = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1594  mFisherScaledEigenvalues = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1595  }
1596 
1598 
1599  // The Fisher Information matrix is a symmetric positive semidefinit matrix.
1600 
1601  /* int dpotrf_(char *uplo, integer *n, doublereal *a,
1602  * integer *lda, integer *info);
1603  *
1604  *
1605  * Purpose
1606  * =======
1607  *
1608  * DPOTRF computes the Cholesky factorization of a real symmetric
1609  * positive definite matrix A.
1610  *
1611  * The factorization has the form
1612  * A = U**T * U, if UPLO = 'U', or
1613  * A = L * L**T, if UPLO = 'L',
1614  * where U is an upper triangular matrix and L is lower triangular.
1615  *
1616  * This is the block version of the algorithm, calling Level 3 BLAS.
1617  *
1618  * Arguments
1619  * =========
1620  *
1621  * UPLO (input) CHARACTER*1
1622  * = 'U': Upper triangle of A is stored;
1623  * = 'L': Lower triangle of A is stored.
1624  *
1625  * N (input) INTEGER
1626  * The order of the matrix A. N >= 0.
1627  *
1628  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
1629  * On entry, the symmetric matrix A. If UPLO = 'U', the leading
1630  * N-by-N upper triangular part of A contains the upper
1631  * triangular part of the matrix A, and the strictly lower
1632  * triangular part of A is not referenced. If UPLO = 'L', the
1633  * leading N-by-N lower triangular part of A contains the lower
1634  * triangular part of the matrix A, and the strictly upper
1635  * triangular part of A is not referenced.
1636  *
1637  * On exit, if INFO = 0, the factor U or L from the Cholesky
1638  * factorization A = U**T*U or A = L*L**T.
1639  *
1640  * LDA (input) INTEGER
1641  * The leading dimension of the array A. LDA >= max(1,N).
1642  *
1643  * INFO (output) INTEGER
1644  * = 0: successful exit
1645  * < 0: if INFO = -i, the i-th argument had an illegal value
1646  * > 0: if INFO = i, the leading minor of order i is not
1647  * positive definite, and the factorization could not be
1648  * completed.
1649  *
1650  */
1651 
1652  char U = 'U';
1653 
1654  dpotrf_(&U, &N, mCorrelation.array(), &N, &INFO);
1655 
1656  if (INFO)
1657  {
1658  mCorrelation = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1659  mParameterSD = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1660 
1662 
1663  // Make sure the timer is accurate.
1664  (*mCPUTime.getRefresh())();
1665 
1666  return false;
1667  }
1668 
1669  /* int dpotri_(char *uplo, integer *n, doublereal *a,
1670  * integer *lda, integer *info);
1671  *
1672  *
1673  * Purpose
1674  * =======
1675  *
1676  * DPOTRI computes the inverse of a real symmetric positive definite
1677  * matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
1678  * computed by DPOTRF.
1679  *
1680  * Arguments
1681  * =========
1682  *
1683  * UPLO (input) CHARACTER*1
1684  * = 'U': Upper triangle of A is stored;
1685  * = 'L': Lower triangle of A is stored.
1686  *
1687  * N (input) INTEGER
1688  * The order of the matrix A. N >= 0.
1689  *
1690  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
1691  * On entry, the triangular factor U or L from the Cholesky
1692  * factorization A = U**T*U or A = L*L**T, as computed by
1693  * DPOTRF.
1694  * On exit, the upper or lower triangle of the (symmetric)
1695  * inverse of A, overwriting the input factor U or L.
1696  *
1697  * LDA (input) INTEGER
1698  * The leading dimension of the array A. LDA >= max(1,N).
1699  *
1700  * INFO (output) INTEGER
1701  * = 0: successful exit
1702  * < 0: if INFO = -i, the i-th argument had an illegal value
1703  * > 0: if INFO = i, the (i,i) element of the factor U or L is
1704  * zero, and the inverse could not be computed.
1705  *
1706  */
1707 
1708  dpotri_(&U, &N, mCorrelation.array(), &N, &INFO);
1709 
1710  if (INFO)
1711  {
1712  mCorrelation = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1713  mParameterSD = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1714 
1716 
1717  // Make sure the timer is accurate.
1718  (*mCPUTime.getRefresh())();
1719 
1720  return false;
1721  }
1722 
1723  // Assure that the inverse is completed.
1724 
1725  for (i = 0; i < imax; i++)
1726  for (l = 0; l < i; l++)
1727  mCorrelation(l, i) = mCorrelation(i, l);
1728 
1729  CVector< C_FLOAT64 > S(imax);
1730 
1731  // rescale the lower bound of the covariant matrix to have unit diagonal
1732  for (i = 0; i < imax; i++)
1733  {
1734  C_FLOAT64 & tmp = S[i];
1735 
1736  if (mCorrelation(i, i) > 0.0)
1737  {
1738  tmp = 1.0 / sqrt(mCorrelation(i, i));
1739  mParameterSD[i] = mSD / tmp;
1740  }
1741  else if (mCorrelation(i, i) < 0.0)
1742  {
1743  tmp = 1.0 / sqrt(- mCorrelation(i, i));
1744  mParameterSD[i] = mSD / tmp;
1745  }
1746  else
1747  {
1749  tmp = 1.0;
1750  mCorrelation(i, i) = 1.0;
1751  }
1752  }
1753 
1754  for (i = 0; i < imax; i++)
1755  for (l = 0; l < imax; l++)
1756  mCorrelation(i, l) *= S[i] * S[l];
1757 
1758  setResidualsRequired(false);
1759  mStoreResults = true;
1760  // This is necessary so that CExperiment::printResult shows the correct data.
1761  calculate();
1762 
1763  // Make sure the timer is accurate.
1764  (*mCPUTime.getRefresh())();
1765  }
1766 
1767  mStoreResults = false;
1768 
1769  return true;
1770 }
1771 
1773 {return mRMS;}
1774 
1776 {return mSD;}
1777 
1779 {return mParameterSD;}
1780 
1782 {return *mpFisherMatrix;}
1783 
1785 {return *mpFisherEigenvaluesMatrix;}
1786 
1788 {return *mpFisherEigenvectorsMatrix;}
1789 
1791 {return *mpFisherScaledMatrix;}
1792 
1795 
1798 
1800 {return *mpCorrelationMatrix;}
1801 
1803 {return *mpExperimentSet;}
1804 
1806 {return *mpCrossValidationSet;}
1807 
1809  const CVector< C_FLOAT64 > & variables)
1810 {
1811  bool Continue = COptProblem::setSolution(value, variables);
1812 
1813  if (Continue && mpCrossValidationSet->getExperimentCount() > 0)
1814  Continue = calculateCrossValidation();
1815 
1816  return Continue;
1817 }
1818 
1821 
1823 {return mCrossValidationRMS;}
1824 
1826 {return mCrossValidationSD;}
1827 
1829 {
1830  mCounter += 1;
1831  bool Continue = true;
1832 
1833  unsigned i, imax = mpCrossValidationSet->getExperimentCount();
1834  unsigned j;
1835  unsigned kmax;
1836  C_FLOAT64 CalculateValue = 0.0;
1837 
1838  CExperiment * pExp = NULL;
1839 
1840  C_FLOAT64 * Residuals = NULL;
1841  C_FLOAT64 * DependentValues = mCrossValidationDependentValues.array();
1842 
1845 
1846  C_FLOAT64 * pSolution = mSolutionVariables.array();
1847  C_FLOAT64 * pSolutionEnd = pSolution + mSolutionVariables.size();
1848 
1849  std::vector<COptItem *>::iterator itConstraint;
1850  std::vector<COptItem *>::iterator endConstraint = mpConstraintItems->end();
1851 
1852  std::vector< Refresh *>::const_iterator itRefresh;
1853  std::vector< Refresh *>::const_iterator endRefresh;
1854 
1855  // Reset the constraints memory
1856  for (itConstraint = mpConstraintItems->begin(); itConstraint != endConstraint; ++itConstraint)
1857  static_cast<CFitConstraint *>(*itConstraint)->resetConstraintViolation();
1858 
1860  CFitConstraint **ppConstraintEnd;
1861 
1862  try
1863  {
1864  for (i = 0; i < imax && Continue; i++) // For each CrossValidation
1865  {
1867 
1870 
1871  // set the global and CrossValidation local fit item values.
1872  for (; pSolution != pSolutionEnd; pSolution++, pUpdate++)
1873  if (*pUpdate)
1874  (**pUpdate)(*pSolution);
1875 
1876  kmax = pExp->getNumDataRows();
1877 
1878  switch (pExp->getExperimentType())
1879  {
1881 
1882  // set independent data
1883  for (j = 0; j < kmax && Continue; j++) // For each data row;
1884  {
1886 
1887  // We need to apply the parameter and independent
1888  // value updates as one unit.
1889  itRefresh = mCrossValidationInitialRefreshes[i].begin();
1890  endRefresh = mCrossValidationInitialRefreshes[i].end();
1891 
1892  for (; itRefresh != endRefresh; ++itRefresh)
1893  (**itRefresh)();
1894 
1895  Continue &= mpSteadyState->process(true);
1896 
1897  if (!Continue)
1898  {
1899  CalculateValue = mWorstValue;
1900  break;
1901  }
1902 
1903  // We check after each simulation whether the constraints are violated.
1904  // Make sure the constraint values are up to date.
1905  itRefresh = mCrossValidationConstraintRefreshes[i].begin();
1906  endRefresh = mCrossValidationConstraintRefreshes[i].end();
1907 
1908  for (; itRefresh != endRefresh; ++itRefresh)
1909  (**itRefresh)();
1910 
1911  ppConstraint = mCrossValidationConstraints[i];
1912  ppConstraintEnd = ppConstraint + mCrossValidationConstraints.numCols();
1913 
1914  for (; ppConstraint != ppConstraintEnd; ++ppConstraint)
1915  if (*ppConstraint)(*ppConstraint)->checkConstraint();
1916 
1917  if (mStoreResults)
1918  CalculateValue += pExp->sumOfSquaresStore(j, DependentValues);
1919  else
1920  CalculateValue += pExp->sumOfSquares(j, Residuals);
1921  }
1922 
1923  break;
1924 
1926 
1927  size_t numIntermediateSteps;
1928 
1929  if (mStoreResults)
1930  {
1931  //calculate a reasonable number of intermediate points
1932  numIntermediateSteps = 4; //TODO
1933  //resize the storage for the extended time series
1934  pExp->initExtendedTimeSeries(numIntermediateSteps * kmax - numIntermediateSteps + 1);
1935  }
1936 
1937  for (j = 0; j < kmax && Continue; j++) // For each data row;
1938  {
1939  if (j)
1940  {
1941  if (mStoreResults)
1942  {
1943  //do additional intermediate steps for nicer display
1944  C_FLOAT64 ttt;
1945  size_t ic;
1946 
1947  for (ic = 1; ic < numIntermediateSteps; ++ic)
1948  {
1949  ttt = pExp->getTimeData()[j - 1] + (pExp->getTimeData()[j] - pExp->getTimeData()[j - 1]) * (C_FLOAT64(ic) / numIntermediateSteps);
1950  mpTrajectory->processStep(ttt);
1951  //save the simulation results in the experiment
1952  pExp->storeExtendedTimeSeriesData(ttt);
1953  }
1954  }
1955 
1956  //do the regular step
1957  mpTrajectory->processStep(pExp->getTimeData()[j]);
1958  }
1959  else
1960  {
1961  // Set independent data. A time course only has one set of
1962  // independent data.
1964 
1965  // We need to apply the parameter and independent
1966  // value updates as one unit.
1967  itRefresh = mCrossValidationInitialRefreshes[i].begin();
1968  endRefresh = mCrossValidationInitialRefreshes[i].end();
1969 
1970  for (; itRefresh != endRefresh; ++itRefresh)
1971  (**itRefresh)();
1972 
1973  static_cast< CTrajectoryProblem * >(mpTrajectory->getProblem())->setStepNumber(1);
1974  mpTrajectory->processStart(true);
1975 
1976  if (pExp->getTimeData()[0] != mpModel->getInitialTime())
1977  {
1978  mpTrajectory->processStep(pExp->getTimeData()[0]);
1979  }
1980  }
1981 
1982  // We check after each simulation whether the constraints are violated.
1983  // Make sure the constraint values are up to date.
1984  itRefresh = mCrossValidationConstraintRefreshes[i].begin();
1985  endRefresh = mCrossValidationConstraintRefreshes[i].end();
1986 
1987  for (; itRefresh != endRefresh; ++itRefresh)
1988  (**itRefresh)();
1989 
1990  ppConstraint = mCrossValidationConstraints[i];
1991  ppConstraintEnd = ppConstraint + mCrossValidationConstraints.numCols();
1992 
1993  for (; ppConstraint != ppConstraintEnd; ++ppConstraint)
1994  if (*ppConstraint)(*ppConstraint)->checkConstraint();
1995 
1996  if (mStoreResults)
1997  CalculateValue += pExp->sumOfSquaresStore(j, DependentValues);
1998  else
1999  CalculateValue += pExp->sumOfSquares(j, Residuals);
2000 
2001  if (mStoreResults)
2002  {
2003  //additionally also store the the simulation result for the extended time series
2004  pExp->storeExtendedTimeSeriesData(pExp->getTimeData()[j]);
2005  }
2006  }
2007 
2008  break;
2009 
2010  default:
2011  break;
2012  }
2013 
2014  // restore independent data
2016 
2017  // restore experiment local values
2018  const C_FLOAT64 *pOriginal = mOriginalVariables.array();
2019  const C_FLOAT64 *pOriginalEnd = pOriginal + mOriginalVariables.size();
2020  bool RefreshNeeded = false;
2021 
2022  // set the global and experiment local fit item values.
2023  for (; pOriginal != pOriginalEnd; pOriginal++, pUndo++)
2024  if (*pUndo)
2025  {
2026  (**pUndo)(*pOriginal);
2027  RefreshNeeded = true;
2028  }
2029 
2030  if (RefreshNeeded)
2031  {
2032  // Update initial values which changed due to the fit item values.
2033  itRefresh = mCrossValidationInitialRefreshes[i].begin();
2034  endRefresh = mCrossValidationInitialRefreshes[i].end();
2035 
2036  while (itRefresh != endRefresh)
2037  (**itRefresh++)();
2038  }
2039  }
2040  }
2041 
2042  catch (CCopasiException &)
2043  {
2044  // We do not want to clog the message cue.
2046 
2047  mFailedCounter++;
2048  CalculateValue = mWorstValue;
2049 
2050  if (pExp)
2051  {
2053 
2054  // Update initial values which changed due to the fit item values.
2055  itRefresh = mCrossValidationInitialRefreshes[i].begin();
2056  endRefresh = mCrossValidationInitialRefreshes[i].end();
2057 
2058  while (itRefresh != endRefresh)
2059  (**itRefresh++)();
2060  }
2061  }
2062 
2063  catch (...)
2064  {
2065  mFailedCounter++;
2066  CalculateValue = mWorstValue;
2067 
2068  if (pExp)
2069  {
2071 
2072  // Update initial values which changed due to the fit item values.
2073  itRefresh = mCrossValidationInitialRefreshes[i].begin();
2074  endRefresh = mCrossValidationInitialRefreshes[i].end();
2075 
2076  while (itRefresh != endRefresh)
2077  (**itRefresh++)();
2078  }
2079  }
2080 
2081  if (isnan(CalculateValue))
2082  CalculateValue = mWorstValue;
2083 
2085  CalculateValue = mWorstValue;
2086 
2087  if (mpCallBack)
2088  Continue &= mpCallBack->progressItem(mhCounter);
2089 
2090  C_FLOAT64 CurrentObjective =
2093 
2094  if (CurrentObjective > mCrossValidationObjective)
2096  else
2097  {
2098  mThresholdCounter = 0;
2099  mCrossValidationObjective = CurrentObjective;
2100  mCrossValidationSolutionValue = CalculateValue;
2101  }
2102 
2104 
2105  return Continue;
2106 }
2107 
2109 {
2110  if (mpExperimentSet != NULL)
2111  {
2113  }
2114 
2115  if (mpCrossValidationSet != NULL)
2116  {
2118  }
2119 
2120  return;
2121 }
CMatrix< C_FLOAT64 > mFisherEigenvectors
Definition: CFitProblem.h:459
CCopasiDataModel * getObjectDataModel()
const std::string & getExperiment(const size_t &index) const
Definition: CFitItem.cpp:194
virtual Refresh * getRefresh() const
CCopasiContainer * getObjectAncestor(const std::string &type) const
CArrayAnnotation * mpFisherMatrix
Definition: CFitProblem.h:450
#define C_INT
Definition: copasi.h:115
virtual bool calculate()
CFitProblem(const CCopasiTask::Type &type=CCopasiTask::parameterFitting, const CCopasiContainer *pParent=NULL)
Definition: CFitProblem.cpp:42
static const CCopasiMessage & peekLastMessage()
C_FLOAT64 mRMS
Definition: CFitProblem.h:433
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpFisherScaledMatrixInterface
Definition: CFitProblem.h:467
CVector< C_FLOAT64 > mSolutionVariables
Definition: COptProblem.h:479
virtual void printResult(std::ostream *ostream) const
bool * mpCreateParameterSets
Definition: CFitProblem.h:491
#define pdelete(p)
Definition: copasi.h:215
C_FLOAT64 mSolutionValue
Definition: COptProblem.h:489
CVector< std::vector< Refresh * > > mCrossValidationInitialRefreshes
Definition: CFitProblem.h:371
CArrayAnnotation * mpFisherEigenvectorsMatrix
Definition: CFitProblem.h:461
std::string * mpParmSteadyStateCN
Definition: CFitProblem.h:296
C_FLOAT64 mCrossValidationObjective
Definition: CFitProblem.h:408
virtual bool initialize(const OutputFlag &of, COutputHandler *pOutputHandler, std::ostream *pOstream)
virtual bool process(const bool &useInitialValues)
unsigned C_INT32 mCounter
Definition: COptProblem.h:494
CCopasiProblem * getProblem()
CMatrix< C_FLOAT64 > mFisher
Definition: CFitProblem.h:448
bool setResidualsRequired(const bool &required)
virtual CCopasiObjectName getCN() const
const C_FLOAT64 & getCrossValidationSolutionValue() const
void updateInitialState()
const CCopasiObject * getObject() const
Definition: COptItem.cpp:128
const std::string & getObjectName() const
bool * mpParmCalculateStatistics
Definition: COptProblem.h:416
virtual size_t size() const
int dpotrf_(char *uplo, integer *n, doublereal *a, integer *lda, integer *info)
std::vector< COptItem * > * mpConstraintItems
Definition: COptProblem.h:436
const bool & isUpdateModel() const
std::vector< UpdateMethod * > mUpdateMethods
Definition: COptProblem.h:451
CMatrix< UpdateMethod * > mExperimentUpdateMethods
Definition: CFitProblem.h:323
CCopasiObject * get(const std::string &key)
CVector< C_FLOAT64 > mGradient
Definition: COptProblem.h:541
unsigned C_INT32 mFailedCounter
Definition: COptProblem.h:499
void setInitialState(const CState &state)
Definition: CModel.cpp:1774
bool processStep(const C_FLOAT64 &nextTime)
#define MCOptimization
#define fatalError()
virtual bool initialize()
C_FLOAT64 mCrossValidationSolutionValue
Definition: CFitProblem.h:392
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpFisherScaledEigenvectorsMatrixInterface
Definition: CFitProblem.h:478
CArrayAnnotation & getScaledFisherInformationEigenvectors() const
CMatrix< C_FLOAT64 > mCorrelation
Definition: CFitProblem.h:484
const CVector< C_FLOAT64 > & getVariableStdDeviations() const
const bool & getCreateParameterSets() const
virtual bool setModel(CModel *pModel)
Definition: CState.h:305
CTrajectoryProblem * mpTrajectoryProblem
Definition: CFitProblem.h:418
std::string * mpParmTimeCourseCN
Definition: CFitProblem.h:301
bool calculateCrossValidation()
virtual bool calculateStatistics(const C_FLOAT64 &factor=1.0e-003, const C_FLOAT64 &resolution=1.0e-009)
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
CState * mpInitialState
Definition: CFitProblem.h:423
const std::set< const CCopasiObject * > & getUptoDateObjects() const
Definition: CModel.cpp:1178
#define C_INVALID_INDEX
Definition: copasi.h:222
CVector< C_FLOAT64 > mCrossValidationDependentValues
Definition: CFitProblem.h:387
const CVector< C_FLOAT64 > & getResiduals() const
const CExperimentSet & getExperiementSet() const
CMatrix< CFitConstraint * > mExperimentConstraints
Definition: CFitProblem.h:339
CExperiment * getExperiment(const size_t &index)
virtual size_t getIndex(const std::string &name) const
CArrayAnnotation * mpFisherEigenvaluesMatrix
Definition: CFitProblem.h:457
std::ostream & operator<<(std::ostream &os, const CFitProblem &o)
void initExtendedTimeSeries(size_t s)
void processStart(const bool &useInitialValues)
CTrajectoryTask * mpTrajectory
Definition: CFitProblem.h:318
CVector< C_FLOAT64 > mExperimentDependentValues
Definition: CFitProblem.h:350
CArrayAnnotation * mpCorrelationMatrix
Definition: CFitProblem.h:486
bool calculateStatistics()
CArrayAnnotation & getCorrelations() const
const std::set< const CCopasiObject * > & getIndependentObjects() const
const C_FLOAT64 & getCrossValidationRMS() const
virtual bool elevateChildren()
virtual bool progressItem(const size_t &handle)
#define MCCopasiMessage
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
unsigned C_INT32 mThresholdCounter
Definition: CFitProblem.h:413
CArrayAnnotation * mpFisherScaledEigenvaluesMatrix
Definition: CFitProblem.h:475
const C_FLOAT64 & getLocalValue() const
Definition: CFitItem.cpp:175
virtual bool restore(const bool &updateModel)
bool mHaveStatistics
Definition: COptProblem.h:536
CMatrix< C_FLOAT64 > mFisherEigenvalues
Definition: CFitProblem.h:455
virtual ~CFitProblem()
unsigned C_INT32 mConstraintCounter
Definition: COptProblem.h:504
CRegisteredObjectName * pCN
const CCopasiTimeVariable & getElapsedTime() const
void setAnnotationString(size_t d, size_t i, const std::string s)
unsigned C_INT32 mFailedConstraintCounter
Definition: COptProblem.h:509
size_t mhCounter
Definition: COptProblem.h:524
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpFisherMatrixInterface
Definition: CFitProblem.h:449
virtual bool restore(const bool &updateModel)
const size_t & getValidValueCount() const
CArrayAnnotation & getFisherInformation() const
CTSSATask * pTask
virtual bool setCallBack(CProcessReport *pCallBack)
size_t getNumDataRows() const
bool removeParameter(const std::string &name)
CVector< std::vector< Refresh * > > mExperimentInitialRefreshes
Definition: CFitProblem.h:334
C_FLOAT64 mCalculateValue
Definition: COptProblem.h:474
C_FLOAT64 mSD
Definition: CFitProblem.h:438
int dpotri_(char *uplo, integer *n, doublereal *a, integer *lda, integer *info)
static std::string LL2String(const C_INT64 &value, const C_INT32 &digits=0)
Definition: CopasiTime.cpp:246
const CDescription & getDescription() const
void initObjects()
std::vector< COptItem * > * mpOptItems
Definition: COptProblem.h:431
CVector< C_FLOAT64 > mParameterSD
Definition: CFitProblem.h:443
void storeExtendedTimeSeriesData(C_FLOAT64 time)
virtual bool add(const CType &src)
const C_FLOAT64 & getWeight() const
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpFisherEigenvaluesMatrixInterface
Definition: CFitProblem.h:456
CArrayAnnotation & getScaledFisherInformation() const
virtual bool setModel(CModel *pModel)
bool restoreModelIndependentData()
bool mStoreResults
Definition: COptProblem.h:530
#define MCFitting
CSteadyStateTask * mpSteadyState
Definition: CFitProblem.h:312
CMatrix< C_FLOAT64 > mFisherScaledEigenvalues
Definition: CFitProblem.h:473
CArrayAnnotation & getFisherInformationEigenvalues() const
std::string * mpParmSubtaskCN
Definition: COptProblem.h:390
CArrayAnnotation * mpFisherScaledEigenvectorsMatrix
Definition: CFitProblem.h:479
const Value & getValue() const
CArrayAnnotation & getFisherInformationEigenvectors() const
const C_FLOAT64 & getInitialTime() const
Definition: CModel.cpp:1184
const CCrossValidationSet & getCrossValidationSet() const
CCopasiVectorN< CCopasiTask > * getTaskList()
bool hasDataForTaskType(const CCopasiTask::Type &type) const
bool setValue(const std::string &name, const CType &value)
bool compile(const std::vector< CCopasiContainer * > listOfContainer=CCopasiContainer::EmptyList)
CCopasiTimer mCPUTime
Definition: COptProblem.h:514
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpFisherScaledEigenvaluesMatrixInterface
Definition: CFitProblem.h:474
CCopasiParameterGroup * assertGroup(const std::string &name)
std::string * mpParmObjectiveExpression
Definition: COptProblem.h:401
CArrayAnnotation & getScaledFisherInformationEigenvalues() const
const CCopasiVectorN< CModelParameterSet > & getModelParameterSets() const
Definition: CModel.cpp:1066
C_FLOAT64 sumOfSquares(const size_t &index, C_FLOAT64 *&residuals) const
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpCorrelationMatrixInterface
Definition: CFitProblem.h:485
CCopasiParameter * getParameter(const std::string &name)
CMatrix< UpdateMethod * > mExperimentUndoMethods
Definition: CFitProblem.h:328
C_FLOAT64 sumOfSquaresStore(const size_t &index, C_FLOAT64 *&dependentValues)
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
CMatrix< C_FLOAT64 > mFisherScaledEigenvectors
Definition: CFitProblem.h:477
CMatrix< CFitConstraint * > mCrossValidationConstraints
Definition: CFitProblem.h:376
virtual bool initialize()
void fixBuild55()
virtual bool elevateChildren()
bool * mpParmMaximize
Definition: COptProblem.h:406
size_t keyToIndex(const std::string &key) const
virtual bool setCallBack(CProcessReport *pCallBack)
C_FLOAT64 mCrossValidationRMS
Definition: CFitProblem.h:397
virtual bool createObjectiveFunction()
size_t size() const
Definition: CVector.h:100
const C_FLOAT64 & getStdDeviation() const
std::set< const CObjectInterface * > ObjectSet
int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a, integer *lda, doublereal *w, doublereal *work, integer *lwork, integer *info)
virtual const DataObjectSet & getDirectDependencies(const DataObjectSet &context=DataObjectSet()) const
const CCopasiTask::Type & getExperimentType() const
virtual std::string getObjectDisplayName(bool regular=true, bool richtext=false) const
const CVector< C_FLOAT64 > & getTimeData() const
size_t getCrossValidationCount() const
Definition: CFitItem.cpp:257
static CKeyFactory * getKeyFactory()
C_INT64 getSeconds(const bool &bounded=false) const
Definition: CopasiTime.cpp:126
virtual bool checkFunctionalConstraints()
Header file of class CArrayAnnotation.
#define C_FLOAT64
Definition: copasi.h:92
const CState & getInitialState() const
Definition: CModel.cpp:1768
CCrossValidationSet * mpCrossValidationSet
Definition: CFitProblem.h:355
virtual bool initialize(const OutputFlag &of, COutputHandler *pOutputHandler, std::ostream *pOstream)
virtual void print(std::ostream *ostream) const
CType * array()
Definition: CVector.h:139
bool updateBounds(std::vector< COptItem * >::iterator it)
Definition: CFitItem.cpp:279
size_t getDataPointCount() const
const unsigned C_INT32 & getThreshold() const
virtual void print(std::ostream *ostream) const
static CCopasiMessage getLastMessage()
const C_FLOAT64 & getCrossValidationSD() const
std::string UTCTimeStamp()
Definition: utility.cpp:64
static std::vector< Refresh * > buildUpdateSequence(const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())
C_FLOAT64 mWorstValue
Definition: COptProblem.h:385
CMatrix< C_FLOAT64 > mFisherScaled
Definition: CFitProblem.h:466
void initializeParameter()
void createParameterSetsForExperiment(CExperiment *pExp)
CVector< std::vector< Refresh * > > mCrossValidationConstraintRefreshes
Definition: CFitProblem.h:382
bool updateInitialValues()
Definition: CModel.cpp:1461
CVector< C_FLOAT64 > mOriginalVariables
Definition: COptProblem.h:484
const CCopasiParameter::Value & getValue(const std::string &name) const
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
Definition: CModel.h:50
virtual bool restore()
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpFisherEigenvectorsMatrixInterface
Definition: CFitProblem.h:460
CModel * mpModel
void setCreateParameterSets(const bool &create)
CCopasiParameterGroup * getGroup(const std::string &name)
CExperimentSet * mpExperimentSet
Definition: CFitProblem.h:306
CArrayAnnotation * mpFisherScaledMatrix
Definition: CFitProblem.h:468
virtual void printResult(std::ostream *ostream) const
virtual size_t numCols() const
Definition: CMatrix.h:144
std::vector< Refresh * > buildInitialRefreshSequence(std::set< const CCopasiObject * > &changedObjects)
Definition: CModel.cpp:4164
bool updateModelWithIndependentData(const size_t &index)
CProcessReport * mpCallBack
std::string getObjectDisplayName() const
Definition: COptItem.cpp:134
bool setObjectName(const std::string &name)
void setUpdateModel(const bool &updateModel)
const C_FLOAT64 & getRMS() const
C_INT64 getMilliSeconds(const bool &bounded=false) const
Definition: CopasiTime.cpp:118
CCopasiObject * ObjectFromName(const std::vector< CCopasiContainer * > &listOfContainer, const CCopasiObjectName &CN) const
CMatrix< UpdateMethod * > mCrossValidationUndoMethods
Definition: CFitProblem.h:365
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
virtual bool restore()
std::string getExperiments() const
Definition: CFitItem.cpp:210
bool mTrajectoryUpdate
Definition: CFitProblem.h:496
virtual CType * array()
Definition: CMatrix.h:337
size_t getExperimentCount() const
const std::string & getCrossValidation(const size_t &index) const
Definition: CFitItem.cpp:244
CMatrix< UpdateMethod * > mCrossValidationUpdateMethods
Definition: CFitProblem.h:360
C_FLOAT64 mCrossValidationSD
Definition: CFitProblem.h:402
size_t getExperimentCount() const
Definition: CFitItem.cpp:207
CCopasiContainer * getObjectParent() const
CVector< C_FLOAT64 > mResiduals
Definition: CFitProblem.h:428
CVector< std::vector< Refresh * > > mExperimentConstraintRefreshes
Definition: CFitProblem.h:345
static CCopasiMessage::Type getHighestSeverity()