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 
25 #include "CopasiDataModel/CCopasiDataModel.h"
26 #include "report/CCopasiRootContainer.h"
27 #include "model/CModel.h"
28 #include "model/CState.h"
29 #include "report/CCopasiObjectReference.h"
30 #include "report/CKeyFactory.h"
31 #include "steadystate/CSteadyStateTask.h"
32 #include "trajectory/CTrajectoryTask.h"
33 #include "trajectory/CTrajectoryProblem.h"
34 #include "utilities/CProcessReport.h"
35 #include "utilities/CCopasiException.h"
36 #include "utilities/CAnnotatedMatrix.h"
37 
38 #include "lapack/blaswrap.h" //use blas
39 #include "lapack/lapackwrap.h" //use CLAPACK
40 
41 // Default constructor
42 CFitProblem::CFitProblem(const CCopasiTask::Type & type,
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();
95  initializeParameter();
96 }
97 
98 // copy constructor
99 CFitProblem::CFitProblem(const CFitProblem& src,
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();
151  initializeParameter();
152 }
153 
154 // Destructor
155 CFitProblem::~CFitProblem()
156 {
157  pdelete(mpTrajectoryProblem);
158  pdelete(mpInitialState);
159  pdelete(mpFisherMatrixInterface);
160  pdelete(mpFisherMatrix);
161  pdelete(mpFisherEigenvaluesMatrixInterface);
162  pdelete(mpFisherEigenvaluesMatrix);
163  pdelete(mpFisherEigenvectorsMatrixInterface);
164  pdelete(mpFisherEigenvectorsMatrix);
165  pdelete(mpFisherScaledMatrixInterface);
166  pdelete(mpFisherScaledMatrix);
167  pdelete(mpFisherScaledEigenvaluesMatrixInterface);
168  pdelete(mpFisherScaledEigenvaluesMatrix);
169  pdelete(mpFisherScaledEigenvectorsMatrixInterface);
170  pdelete(mpFisherScaledEigenvectorsMatrix);
171  pdelete(mpCorrelationMatrixInterface);
172  pdelete(mpCorrelationMatrix);
173 }
174 
175 void CFitProblem::initObjects()
176 {
177  addObjectReference("Validation Solution", mCrossValidationSolutionValue, CCopasiObject::ValueDbl);
178  addObjectReference("Validation Objective", mCrossValidationObjective, CCopasiObject::ValueDbl);
179 
180  mpFisherMatrixInterface = new CCopasiMatrixInterface< CMatrix< C_FLOAT64 > >(&mFisher);
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 
187  mpFisherEigenvaluesMatrixInterface = new CCopasiMatrixInterface< CMatrix< C_FLOAT64 > >(&mFisherEigenvalues);
188  mpFisherEigenvaluesMatrix = new CArrayAnnotation("FIM Eigenvalues", this, mpFisherEigenvaluesMatrixInterface, false);
189  mpFisherEigenvaluesMatrix->setDescription("FIM Eigenvalues");
190  mpFisherEigenvaluesMatrix->setDimensionDescription(0, "Eigenvalues");
191  mpFisherEigenvaluesMatrix->setDimensionDescription(1, "Result");
192  mpFisherEigenvaluesMatrix->setMode(CArrayAnnotation::NUMBERS);
193 
194  mpFisherEigenvectorsMatrixInterface = new CCopasiMatrixInterface< CMatrix< C_FLOAT64 > >(&mFisherEigenvectors);
195  mpFisherEigenvectorsMatrix = new CArrayAnnotation("FIM Eigenvectors", this, mpFisherEigenvectorsMatrixInterface, false);
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 
202  mpFisherScaledMatrixInterface = new CCopasiMatrixInterface< CMatrix< C_FLOAT64 > >(&mFisherScaled);
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 
209  mpFisherScaledEigenvaluesMatrixInterface = new CCopasiMatrixInterface< CMatrix< C_FLOAT64 > >(&mFisherScaledEigenvalues);
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 
216  mpFisherScaledEigenvectorsMatrixInterface = new CCopasiMatrixInterface< CMatrix< C_FLOAT64 > >(&mFisherScaledEigenvectors);
217  mpFisherScaledEigenvectorsMatrix = new CArrayAnnotation("FIM Eigenvectors (scaled)", this, mpFisherScaledEigenvectorsMatrixInterface, false);
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 
224  mpCorrelationMatrixInterface = new CCopasiMatrixInterface< CMatrix< C_FLOAT64 > >(&mCorrelation);
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 
232 void CFitProblem::initializeParameter()
233 {
234  removeParameter("Subtask");
235  mpParmSubtaskCN = NULL;
236  removeParameter("ObjectiveExpression");
237  mpParmObjectiveExpression = NULL;
238  *mpParmMaximize = false;
239 
240  mpParmSteadyStateCN =
241  assertParameter("Steady-State", CCopasiParameter::CN, CCopasiObjectName(""))->getValue().pCN;
242  mpParmTimeCourseCN =
243  assertParameter("Time-Course", CCopasiParameter::CN, CCopasiObjectName(""))->getValue().pCN;
244 
245  mpCreateParameterSets =
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 
258 const bool & CFitProblem::getCreateParameterSets() const
259 {return *mpCreateParameterSets;}
260 
261 bool CFitProblem::elevateChildren()
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.
268  mpParmSteadyStateCN =
269  assertParameter("Steady-State", CCopasiParameter::CN, CCopasiObjectName(""))->getValue().pCN;
270  mpParmTimeCourseCN =
271  assertParameter("Time-Course", CCopasiParameter::CN, CCopasiObjectName(""))->getValue().pCN;
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 
319  mpExperimentSet =
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"));
357  mpCrossValidationSet =
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 
414 bool CFitProblem::setModel(CModel * pModel)
415 {return COptProblem::setModel(pModel);}
416 
417 bool CFitProblem::setCallBack(CProcessReport * pCallBack)
418 {return COptProblem::setCallBack(pCallBack);}
419 
420 bool CFitProblem::initialize()
421 {
422  bool success = true;
423 
424  mHaveStatistics = false;
425  mStoreResults = false;
426 
427  if (!COptProblem::initialize())
428  {
429  while (CCopasiMessage::peekLastMessage().getNumber() == MCOptimization + 5 ||
430  CCopasiMessage::peekLastMessage().getNumber() == MCOptimization + 7)
431  CCopasiMessage::getLastMessage();
432 
433  if (CCopasiMessage::getHighestSeverity() > CCopasiMessage::WARNING &&
434  CCopasiMessage::peekLastMessage().getNumber() != MCCopasiMessage + 1)
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.
445  if (mpExperimentSet->hasDataForTaskType(CCopasiTask::steadyState))
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 
458  *mpParmSteadyStateCN = mpSteadyState->getCN();
459  mpSteadyState->initialize(CCopasiTask::NO_OUTPUT, NULL, NULL);
460  }
461  else
462  {
463  mpSteadyState = NULL;
464  }
465 
466  pdelete(mpTrajectoryProblem);
467 
468  // We only need to initialize the trajectory task if time course data is present.
469  if (mpExperimentSet->hasDataForTaskType(CCopasiTask::timeCourse))
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 
482  *mpParmTimeCourseCN = mpTrajectory->getCN();
483 
484  // do not update initial values when running fit
485  mTrajectoryUpdate = mpTrajectory->isUpdateModel();
486  mpTrajectory->setUpdateModel(false);
487 
488  mpTrajectory->initialize(CCopasiTask::NO_OUTPUT, NULL, NULL);
489 
490  mpTrajectoryProblem =
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  {
528  ObjectSet[i] = mpExperimentSet->getExperiment(i)->getIndependentObjects();
529  }
530 
531  // Build a matrix of experiment and experiment local items.
532  mExperimentUpdateMethods.resize(mpExperimentSet->getExperimentCount(),
533  mpOptItems->size());
534  mExperimentUpdateMethods = NULL;
535 
536  // Build a matrix of experiment and experiment local undo items.
537  mExperimentUndoMethods.resize(mpExperimentSet->getExperimentCount(),
538  mpOptItems->size());
539  mExperimentUndoMethods = NULL;
540 
541  mExperimentInitialRefreshes.resize(mpExperimentSet->getExperimentCount());
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);
555  mpFisherMatrix->resize();
556  mFisherEigenvalues.resize(imax, 1);
557  mpFisherEigenvaluesMatrix->resize();
558  mFisherEigenvectors.resize(imax, imax);
559  mpFisherEigenvectorsMatrix->resize();
560  mFisherScaled.resize(imax, imax);
561  mpFisherScaledMatrix->resize();
562  mFisherScaledEigenvalues.resize(imax, 1);
563  mpFisherScaledEigenvaluesMatrix->resize();
564  mFisherScaledEigenvectors.resize(imax, imax);
565  mpFisherScaledEigenvectorsMatrix->resize();
566  mCorrelation.resize(imax, imax);
567  mpCorrelationMatrix->resize();
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);
606  mpFisherEigenvectorsMatrix->setAnnotationString(1, j, Annotation);
607  mpFisherScaledMatrix->setAnnotationString(0, j, Annotation);
608  mpFisherScaledMatrix->setAnnotationString(1, j, Annotation);
609  mpFisherScaledEigenvectorsMatrix->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);
618  mExperimentInitialRefreshes[i] = mpModel->buildInitialRefreshSequence(ObjectSet[i]);
619  }
620 
621  // Build a matrix of experiment and constraint items;
622  mExperimentConstraints.resize(mpExperimentSet->getExperimentCount(),
623  mpConstraintItems->size());
624  mExperimentConstraints = NULL;
625  mExperimentConstraintRefreshes.resize(mpExperimentSet->getExperimentCount());
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++)
665  mExperimentConstraintRefreshes[i] = CCopasiObject::buildUpdateSequence(ObjectSet[i], mpModel->getUptoDateObjects());
666 
667  mExperimentDependentValues.resize(mpExperimentSet->getDataPointCount());
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  {
677  ObjectSet[i] = mpCrossValidationSet->getExperiment(i)->getIndependentObjects();
678  }
679 
680  // Build a matrix of cross validation experiments and local items.
681  mCrossValidationUpdateMethods.resize(mpCrossValidationSet->getExperimentCount(),
682  mpOptItems->size());
683  mCrossValidationUpdateMethods = NULL;
684 
685  // Build a matrix of experiment and experiment local undo items.
686  mCrossValidationUndoMethods.resize(mpCrossValidationSet->getExperimentCount(),
687  mpOptItems->size());
688  mCrossValidationUndoMethods = NULL;
689 
690  mCrossValidationInitialRefreshes.resize(mpCrossValidationSet->getExperimentCount());
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  {
714  if ((Index = mpCrossValidationSet->keyToIndex(pItem->getCrossValidation(i))) == C_INVALID_INDEX)
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  {
730  mCrossValidationInitialRefreshes[i] = mpModel->buildInitialRefreshSequence(ObjectSet[i]);
731  }
732 
733  // Build a matrix of cross validation experiments and constraint items;
734  mCrossValidationConstraints.resize(mpCrossValidationSet->getExperimentCount(),
735  mpConstraintItems->size());
736  mCrossValidationConstraints = NULL;
737  mCrossValidationConstraintRefreshes.resize(mpCrossValidationSet->getExperimentCount());
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++)
772  mCrossValidationConstraintRefreshes[i] = CCopasiObject::buildUpdateSequence(ObjectSet[i], mpModel->getUptoDateObjects());
773 
774  mCrossValidationDependentValues.resize(mpCrossValidationSet->getDataPointCount());
775 
776  mCrossValidationObjective = mWorstValue;
777  mThresholdCounter = 0;
778 
779  setResidualsRequired(false);
780 
781  pdelete(mpInitialState);
782  mpInitialState = new CState(mpModel->getInitialState());
783 
784  return true;
785 }
786 
787 bool CFitProblem::checkFunctionalConstraints()
788 {
789  std::vector< COptItem * >::const_iterator it = mpConstraintItems->begin();
790  std::vector< COptItem * >::const_iterator end = mpConstraintItems->end();
791 
792  mConstraintCounter++;
793 
794  for (; it != end; ++it)
795  if (static_cast<CFitConstraint *>(*it)->getConstraintViolation() > 0.0)
796  {
797  mFailedConstraintCounter++;
798  return false;
799  }
800 
801  return true;
802 }
803 
804 /**
805  * Utility function creating a parameter set for each experiment
806  */
807 void createParameterSetsForExperiment(CExperiment* pExp)
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 
827 bool CFitProblem::calculate()
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();
841  UpdateMethod ** pUpdate = mExperimentUpdateMethods.array();
842  UpdateMethod ** pUndo = mExperimentUndoMethods.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.
866  mpModel->setInitialState(*mpInitialState);
867  mpModel->updateInitialValues();
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  {
878  case CCopasiTask::steadyState:
879 
880  // set independent data
881  for (j = 0; j < kmax && Continue; j++) // For each data row;
882  {
883  pExp->updateModelWithIndependentData(j);
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++;
898  mCalculateValue = mWorstValue;
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 
922  if (mStoreResults && *mpCreateParameterSets)
923  {
924  createParameterSetsForExperiment(pExp);
925  }
926 
927  break;
928 
929  case CCopasiTask::timeCourse:
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);
954  mpTrajectory->processStep(ttt);
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.
967  pExp->updateModelWithIndependentData(0);
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 
1012  if (mStoreResults && *mpCreateParameterSets)
1013  {
1014  createParameterSetsForExperiment(pExp);
1015  }
1016 
1017  break;
1018 
1019  default:
1020  break;
1021  }
1022 
1023  // restore independent data
1024  pExp->restoreModelIndependentData();
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.
1054  CCopasiMessage::getLastMessage();
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
1066  mCalculateValue = mWorstValue;
1067 
1068  if (pExp)
1069  {
1070  pExp->restoreModelIndependentData();
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++;
1084  mCalculateValue = mWorstValue;
1085 
1086  if (pExp)
1087  {
1088  pExp->restoreModelIndependentData();
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))
1099  mCalculateValue = mWorstValue;
1100 
1101  if (mpCallBack) return mpCallBack->progressItem(mhCounter);
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();
1113  mpTrajectory->setUpdateModel(mTrajectoryUpdate);
1114  }
1115 
1116  if (mpSteadyState != NULL)
1117  success &= mpSteadyState->restore();
1118 
1119  if (mpTrajectoryProblem)
1120  *mpTrajectory->getProblem() = *mpTrajectoryProblem;
1121 
1122  success &= COptProblem::restore(updateModel);
1123 
1124  pdelete(mpTrajectoryProblem);
1125  pdelete(mpInitialState);
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  {
1211  mpExperimentSet->getExperiment(k)->printResult(ostream);
1212  os << std::endl;
1213  }
1214 
1215  if (*mpParmCalculateStatistics)
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)
1247  o.mpSteadyState->getDescription().print(&os);
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 
1276 void CFitProblem::updateInitialState()
1277 {
1278  *mpInitialState = mpModel->getInitialState();
1279 }
1280 
1281 bool CFitProblem::createObjectiveFunction()
1282 {return true;}
1283 
1284 bool CFitProblem::setResidualsRequired(const bool & required)
1285 {
1286  if (required)
1287  {
1288  mResiduals.resize(mpExperimentSet->getDataPointCount(), false);
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 
1303 const CVector< C_FLOAT64 > & CFitProblem::getResiduals() const
1304 {
1305  return mResiduals;
1306 }
1307 
1308 bool CFitProblem::calculateStatistics(const C_FLOAT64 & factor,
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++)
1335  (*mUpdateMethods[i])(mSolutionVariables[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.
1343  mpExperimentSet->calculateStatistics();
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 
1352  calculateCrossValidation();
1353 
1354  mpCrossValidationSet->calculateStatistics();
1355 
1356  size_t lmax = this->mCrossValidationDependentValues.size();
1357 
1358  if (lmax)
1359  mCrossValidationRMS = sqrt(mCrossValidationSolutionValue / 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 
1372  if (*mpParmCalculateStatistics)
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 
1436  CCopasiMessage(CCopasiMessage::WARNING, MCFitting + 13);
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 
1525  mFisherEigenvectors = mFisher;
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,
1539  mFisherEigenvectors.array(),
1540  &LDA,
1541  mFisherEigenvalues.array(),
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,
1553  mFisherEigenvectors.array(),
1554  &LDA,
1555  mFisherEigenvalues.array(),
1556  WORK.array(),
1557  &LWORK,
1558  &INFO);
1559 
1560  if (INFO != 0)
1561  {
1562  CCopasiMessage(CCopasiMessage::WARNING, MCFitting + 14);
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  {
1572  mFisherScaled[i][j] = mFisher[i][j] * mSolutionVariables[i] * mSolutionVariables[j];
1573  }
1574  }
1575 
1576  mFisherScaledEigenvectors = mFisherScaled;
1577 
1578  //now do the real calculation
1579  dsyev_(&JOBZ,
1580  &UPLO,
1581  &N,
1582  mFisherScaledEigenvectors.array(),
1583  &LDA,
1584  mFisherScaledEigenvalues.array(),
1585  WORK.array(),
1586  &LWORK,
1587  &INFO);
1588 
1589  if (INFO != 0)
1590  {
1591  CCopasiMessage(CCopasiMessage::WARNING, MCFitting + 15);
1592 
1593  mFisherScaledEigenvectors = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1594  mFisherScaledEigenvalues = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1595  }
1596 
1597  mCorrelation = mFisher;
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 
1661  CCopasiMessage(CCopasiMessage::WARNING, MCFitting + 12);
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 
1715  CCopasiMessage(CCopasiMessage::WARNING, MCFitting + 1, INFO);
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  {
1748  mParameterSD[i] = mWorstValue;
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 
1772 const C_FLOAT64 & CFitProblem::getRMS() const
1773 {return mRMS;}
1774 
1775 const C_FLOAT64 & CFitProblem::getStdDeviation() const
1776 {return mSD;}
1777 
1778 const CVector< C_FLOAT64 > & CFitProblem::getVariableStdDeviations() const
1779 {return mParameterSD;}
1780 
1781 CArrayAnnotation & CFitProblem::getFisherInformation() const
1782 {return *mpFisherMatrix;}
1783 
1784 CArrayAnnotation & CFitProblem::getFisherInformationEigenvalues() const
1785 {return *mpFisherEigenvaluesMatrix;}
1786 
1787 CArrayAnnotation & CFitProblem::getFisherInformationEigenvectors() const
1788 {return *mpFisherEigenvectorsMatrix;}
1789 
1790 CArrayAnnotation & CFitProblem::getScaledFisherInformation() const
1791 {return *mpFisherScaledMatrix;}
1792 
1793 CArrayAnnotation & CFitProblem::getScaledFisherInformationEigenvalues() const
1794 {return *mpFisherScaledEigenvaluesMatrix;}
1795 
1796 CArrayAnnotation & CFitProblem::getScaledFisherInformationEigenvectors() const
1797 {return *mpFisherScaledEigenvectorsMatrix;}
1798 
1799 CArrayAnnotation & CFitProblem::getCorrelations() const
1800 {return *mpCorrelationMatrix;}
1801 
1802 const CExperimentSet & CFitProblem::getExperiementSet() const
1803 {return *mpExperimentSet;}
1804 
1805 const CCrossValidationSet & CFitProblem::getCrossValidationSet() const
1806 {return *mpCrossValidationSet;}
1807 
1808 bool CFitProblem::setSolution(const C_FLOAT64 & value,
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 
1819 const C_FLOAT64 & CFitProblem::getCrossValidationSolutionValue() const
1820 {return mCrossValidationSolutionValue;}
1821 
1822 const C_FLOAT64 & CFitProblem::getCrossValidationRMS() const
1823 {return mCrossValidationRMS;}
1824 
1825 const C_FLOAT64 & CFitProblem::getCrossValidationSD() const
1826 {return mCrossValidationSD;}
1827 
1828 bool CFitProblem::calculateCrossValidation()
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 
1843  UpdateMethod ** pUpdate = mCrossValidationUpdateMethods.array();
1844  UpdateMethod ** pUndo = mCrossValidationUndoMethods.array();
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 
1859  CFitConstraint **ppConstraint = mCrossValidationConstraints.array();
1860  CFitConstraint **ppConstraintEnd;
1861 
1862  try
1863  {
1864  for (i = 0; i < imax && Continue; i++) // For each CrossValidation
1865  {
1866  pExp = mpCrossValidationSet->getExperiment(i);
1867 
1868  mpModel->setInitialState(*mpInitialState);
1869  mpModel->updateInitialValues();
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  {
1880  case CCopasiTask::steadyState:
1881 
1882  // set independent data
1883  for (j = 0; j < kmax && Continue; j++) // For each data row;
1884  {
1885  pExp->updateModelWithIndependentData(j);
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 
1925  case CCopasiTask::timeCourse:
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.
1963  pExp->updateModelWithIndependentData(0);
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
2015  pExp->restoreModelIndependentData();
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.
2045  CCopasiMessage::getLastMessage();
2046 
2047  mFailedCounter++;
2048  CalculateValue = mWorstValue;
2049 
2050  if (pExp)
2051  {
2052  pExp->restoreModelIndependentData();
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  {
2070  pExp->restoreModelIndependentData();
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 
2084  if (!checkFunctionalConstraints())
2085  CalculateValue = mWorstValue;
2086 
2087  if (mpCallBack)
2088  Continue &= mpCallBack->progressItem(mhCounter);
2089 
2090  C_FLOAT64 CurrentObjective =
2091  (1.0 - mpCrossValidationSet->getWeight()) * mSolutionValue
2092  + mpCrossValidationSet->getWeight() * CalculateValue * mpCrossValidationSet->getDataPointCount() / mpExperimentSet->getDataPointCount();
2093 
2094  if (CurrentObjective > mCrossValidationObjective)
2095  mThresholdCounter++;
2096  else
2097  {
2098  mThresholdCounter = 0;
2099  mCrossValidationObjective = CurrentObjective;
2100  mCrossValidationSolutionValue = CalculateValue;
2101  }
2102 
2103  Continue &= (mThresholdCounter < mpCrossValidationSet->getThreshold());
2104 
2105  return Continue;
2106 }
2107 
2108 void CFitProblem::fixBuild55()
2109 {
2110  if (mpExperimentSet != NULL)
2111  {
2112  mpExperimentSet->fixBuild55();
2113  }
2114 
2115  if (mpCrossValidationSet != NULL)
2116  {
2117  mpCrossValidationSet->fixBuild55();
2118  }
2119 
2120  return;
2121 }
CFitProblem::mFisherEigenvectors
CMatrix< C_FLOAT64 > mFisherEigenvectors
Definition: CFitProblem.h:459
CCopasiObject::getObjectDataModel
CCopasiDataModel * getObjectDataModel()
Definition: CCopasiObject.cpp:666
CFitItem::getExperiment
const std::string & getExperiment(const size_t &index) const
Definition: CFitItem.cpp:194
CCopasiObject::getRefresh
virtual Refresh * getRefresh() const
Definition: CCopasiObject.cpp:633
CCopasiObject::getObjectAncestor
CCopasiContainer * getObjectAncestor(const std::string &type) const
Definition: CCopasiObject.cpp:279
CFitProblem::mpFisherMatrix
CArrayAnnotation * mpFisherMatrix
Definition: CFitProblem.h:450
C_INT
#define C_INT
Definition: copasi.h:115
CFitProblem::calculate
virtual bool calculate()
Definition: CFitProblem.cpp:827
CFitProblem::CFitProblem
CFitProblem(const CCopasiTask::Type &type=CCopasiTask::parameterFitting, const CCopasiContainer *pParent=NULL)
Definition: CFitProblem.cpp:42
CCopasiMessage::peekLastMessage
static const CCopasiMessage & peekLastMessage()
Definition: CCopasiMessage.cpp:60
CFitProblem::mRMS
C_FLOAT64 mRMS
Definition: CFitProblem.h:433
CFitProblem::mpFisherScaledMatrixInterface
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpFisherScaledMatrixInterface
Definition: CFitProblem.h:467
COptProblem::mSolutionVariables
CVector< C_FLOAT64 > mSolutionVariables
Definition: COptProblem.h:479
CFitProblem::printResult
virtual void printResult(std::ostream *ostream) const
Definition: CFitProblem.cpp:1133
CFitProblem::mpCreateParameterSets
bool * mpCreateParameterSets
Definition: CFitProblem.h:491
pdelete
#define pdelete(p)
Definition: copasi.h:215
COptProblem::mSolutionValue
C_FLOAT64 mSolutionValue
Definition: COptProblem.h:489
CFitProblem::mCrossValidationInitialRefreshes
CVector< std::vector< Refresh * > > mCrossValidationInitialRefreshes
Definition: CFitProblem.h:371
CFitProblem::mpFisherEigenvectorsMatrix
CArrayAnnotation * mpFisherEigenvectorsMatrix
Definition: CFitProblem.h:461
CFitProblem::mpParmSteadyStateCN
std::string * mpParmSteadyStateCN
Definition: CFitProblem.h:296
CFitProblem::mCrossValidationObjective
C_FLOAT64 mCrossValidationObjective
Definition: CFitProblem.h:408
CTrajectoryTask::initialize
virtual bool initialize(const OutputFlag &of, COutputHandler *pOutputHandler, std::ostream *pOstream)
Definition: CTrajectoryTask.cpp:161
CSteadyStateTask::process
virtual bool process(const bool &useInitialValues)
Definition: CSteadyStateTask.cpp:301
COptProblem::mCounter
unsigned C_INT32 mCounter
Definition: COptProblem.h:494
CCopasiTask::getProblem
CCopasiProblem * getProblem()
Definition: CCopasiTask.cpp:296
CFitProblem::mFisher
CMatrix< C_FLOAT64 > mFisher
Definition: CFitProblem.h:448
CFitProblem::setResidualsRequired
bool setResidualsRequired(const bool &required)
Definition: CFitProblem.cpp:1284
CCopasiObject::getCN
virtual CCopasiObjectName getCN() const
Definition: CCopasiObject.cpp:103
CFitProblem::getCrossValidationSolutionValue
const C_FLOAT64 & getCrossValidationSolutionValue() const
Definition: CFitProblem.cpp:1819
CFitProblem::updateInitialState
void updateInitialState()
Definition: CFitProblem.cpp:1276
COptItem::getObject
const CCopasiObject * getObject() const
Definition: COptItem.cpp:128
CCopasiObject::getObjectName
const std::string & getObjectName() const
Definition: CCopasiObject.cpp:258
COptProblem::mpParmCalculateStatistics
bool * mpParmCalculateStatistics
Definition: COptProblem.h:416
CCopasiVector::size
virtual size_t size() const
Definition: CCopasiVector.h:386
dpotrf_
int dpotrf_(char *uplo, integer *n, doublereal *a, integer *lda, integer *info)
COptProblem::mpConstraintItems
std::vector< COptItem * > * mpConstraintItems
Definition: COptProblem.h:436
CCopasiTask::isUpdateModel
const bool & isUpdateModel() const
Definition: CCopasiTask.cpp:188
COptProblem::mUpdateMethods
std::vector< UpdateMethod * > mUpdateMethods
Definition: COptProblem.h:451
CFitProblem::mExperimentUpdateMethods
CMatrix< UpdateMethod * > mExperimentUpdateMethods
Definition: CFitProblem.h:323
CKeyFactory::get
CCopasiObject * get(const std::string &key)
Definition: CKeyFactory.cpp:235
COptProblem::mGradient
CVector< C_FLOAT64 > mGradient
Definition: COptProblem.h:541
COptProblem::mFailedCounter
unsigned C_INT32 mFailedCounter
Definition: COptProblem.h:499
CModel::setInitialState
void setInitialState(const CState &state)
Definition: CModel.cpp:1774
CTrajectoryTask::processStep
bool processStep(const C_FLOAT64 &nextTime)
Definition: CTrajectoryTask.cpp:364
MCOptimization
#define MCOptimization
Definition: CCopasiMessage.h:45
fatalError
#define fatalError()
Definition: CCopasiMessage.h:74
COptProblem::initialize
virtual bool initialize()
Definition: COptProblem.cpp:328
CFitProblem::mCrossValidationSolutionValue
C_FLOAT64 mCrossValidationSolutionValue
Definition: CFitProblem.h:392
CFitProblem::mpFisherScaledEigenvectorsMatrixInterface
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpFisherScaledEigenvectorsMatrixInterface
Definition: CFitProblem.h:478
CFitProblem::getScaledFisherInformationEigenvectors
CArrayAnnotation & getScaledFisherInformationEigenvectors() const
Definition: CFitProblem.cpp:1796
CFitProblem::mCorrelation
CMatrix< C_FLOAT64 > mCorrelation
Definition: CFitProblem.h:484
CFitProblem::getVariableStdDeviations
const CVector< C_FLOAT64 > & getVariableStdDeviations() const
Definition: CFitProblem.cpp:1778
CFitProblem::getCreateParameterSets
const bool & getCreateParameterSets() const
Definition: CFitProblem.cpp:258
COptProblem::setModel
virtual bool setModel(CModel *pModel)
Definition: COptProblem.cpp:263
CState
Definition: CState.h:305
CFitProblem::mpTrajectoryProblem
CTrajectoryProblem * mpTrajectoryProblem
Definition: CFitProblem.h:418
CFitProblem::mpParmTimeCourseCN
std::string * mpParmTimeCourseCN
Definition: CFitProblem.h:301
CFitProblem::calculateCrossValidation
bool calculateCrossValidation()
Definition: CFitProblem.cpp:1828
CFitProblem::calculateStatistics
virtual bool calculateStatistics(const C_FLOAT64 &factor=1.0e-003, const C_FLOAT64 &resolution=1.0e-009)
Definition: CFitProblem.cpp:1308
CVector::resize
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
CFitProblem::mpInitialState
CState * mpInitialState
Definition: CFitProblem.h:423
CModel::getUptoDateObjects
const std::set< const CCopasiObject * > & getUptoDateObjects() const
Definition: CModel.cpp:1178
C_INVALID_INDEX
#define C_INVALID_INDEX
Definition: copasi.h:222
CFitProblem::mCrossValidationDependentValues
CVector< C_FLOAT64 > mCrossValidationDependentValues
Definition: CFitProblem.h:387
CFitProblem::getResiduals
const CVector< C_FLOAT64 > & getResiduals() const
Definition: CFitProblem.cpp:1303
CFitProblem::getExperiementSet
const CExperimentSet & getExperiementSet() const
Definition: CFitProblem.cpp:1802
CFitProblem::mExperimentConstraints
CMatrix< CFitConstraint * > mExperimentConstraints
Definition: CFitProblem.h:339
CExperimentSet::getExperiment
CExperiment * getExperiment(const size_t &index)
Definition: CExperimentSet.cpp:352
CCopasiVectorN::getIndex
virtual size_t getIndex(const std::string &name) const
Definition: CCopasiVector.h:774
CFitProblem::mpFisherEigenvaluesMatrix
CArrayAnnotation * mpFisherEigenvaluesMatrix
Definition: CFitProblem.h:457
operator<<
std::ostream & operator<<(std::ostream &os, const CFitProblem &o)
Definition: CFitProblem.cpp:1240
CExperiment::initExtendedTimeSeries
void initExtendedTimeSeries(size_t s)
Definition: CExperiment.cpp:569
CTrajectoryTask::processStart
void processStart(const bool &useInitialValues)
Definition: CTrajectoryTask.cpp:351
CFitProblem::mpTrajectory
CTrajectoryTask * mpTrajectory
Definition: CFitProblem.h:318
CFitProblem::mExperimentDependentValues
CVector< C_FLOAT64 > mExperimentDependentValues
Definition: CFitProblem.h:350
CFitProblem::mpCorrelationMatrix
CArrayAnnotation * mpCorrelationMatrix
Definition: CFitProblem.h:486
CExperimentSet::calculateStatistics
bool calculateStatistics()
Definition: CExperimentSet.cpp:177
CFitProblem::getCorrelations
CArrayAnnotation & getCorrelations() const
Definition: CFitProblem.cpp:1799
CExperiment::getIndependentObjects
const std::set< const CCopasiObject * > & getIndependentObjects() const
Definition: CExperiment.cpp:1676
CFitProblem::getCrossValidationRMS
const C_FLOAT64 & getCrossValidationRMS() const
Definition: CFitProblem.cpp:1822
CFitProblem::elevateChildren
virtual bool elevateChildren()
Definition: CFitProblem.cpp:261
CProcessReport::progressItem
virtual bool progressItem(const size_t &handle)
Definition: CProcessReport.cpp:160
MCCopasiMessage
#define MCCopasiMessage
Definition: CCopasiMessage.h:43
CFitProblem::setSolution
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
Definition: CFitProblem.cpp:1808
CFitProblem::mThresholdCounter
unsigned C_INT32 mThresholdCounter
Definition: CFitProblem.h:413
CFitProblem::mpFisherScaledEigenvaluesMatrix
CArrayAnnotation * mpFisherScaledEigenvaluesMatrix
Definition: CFitProblem.h:475
CFitItem::getLocalValue
const C_FLOAT64 & getLocalValue() const
Definition: CFitItem.cpp:175
CFitProblem::restore
virtual bool restore(const bool &updateModel)
Definition: CFitProblem.cpp:1106
COptProblem::mHaveStatistics
bool mHaveStatistics
Definition: COptProblem.h:536
CFitProblem::mFisherEigenvalues
CMatrix< C_FLOAT64 > mFisherEigenvalues
Definition: CFitProblem.h:455
CFitProblem::~CFitProblem
virtual ~CFitProblem()
Definition: CFitProblem.cpp:155
COptProblem::mConstraintCounter
unsigned C_INT32 mConstraintCounter
Definition: COptProblem.h:504
CCopasiParameter::Value::pCN
CRegisteredObjectName * pCN
Definition: CCopasiParameter.h:64
CCopasiTimer::getElapsedTime
const CCopasiTimeVariable & getElapsedTime() const
Definition: CCopasiTimer.cpp:111
CArrayAnnotation::setAnnotationString
void setAnnotationString(size_t d, size_t i, const std::string s)
Definition: CAnnotatedMatrix.cpp:106
COptProblem::mFailedConstraintCounter
unsigned C_INT32 mFailedConstraintCounter
Definition: COptProblem.h:509
COptProblem::mhCounter
size_t mhCounter
Definition: COptProblem.h:524
CFitProblem::mpFisherMatrixInterface
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpFisherMatrixInterface
Definition: CFitProblem.h:449
COptProblem::restore
virtual bool restore(const bool &updateModel)
Definition: COptProblem.cpp:478
CExperimentSet::getValidValueCount
const size_t & getValidValueCount() const
Definition: CExperimentSet.cpp:461
CFitProblem::getFisherInformation
CArrayAnnotation & getFisherInformation() const
Definition: CFitProblem.cpp:1781
pTask
CTSSATask * pTask
Definition: CQTSSAResultWidget.cpp:31
COptProblem::setCallBack
virtual bool setCallBack(CProcessReport *pCallBack)
Definition: COptProblem.cpp:275
CExperiment::getNumDataRows
size_t getNumDataRows() const
Definition: CExperiment.cpp:1347
CCopasiParameterGroup::removeParameter
bool removeParameter(const std::string &name)
Definition: CCopasiParameterGroup.cpp:309
CFitProblem::mExperimentInitialRefreshes
CVector< std::vector< Refresh * > > mExperimentInitialRefreshes
Definition: CFitProblem.h:334
COptProblem::mCalculateValue
C_FLOAT64 mCalculateValue
Definition: COptProblem.h:474
CFitProblem::mSD
C_FLOAT64 mSD
Definition: CFitProblem.h:438
dpotri_
int dpotri_(char *uplo, integer *n, doublereal *a, integer *lda, integer *info)
CCopasiTimeVariable::LL2String
static std::string LL2String(const C_INT64 &value, const C_INT32 &digits=0)
Definition: CopasiTime.cpp:246
CCopasiTask::getDescription
const CDescription & getDescription() const
Definition: CCopasiTask.cpp:318
CFitProblem::initObjects
void initObjects()
Definition: CFitProblem.cpp:175
COptProblem::mpOptItems
std::vector< COptItem * > * mpOptItems
Definition: COptProblem.h:431
CFitProblem::mParameterSD
CVector< C_FLOAT64 > mParameterSD
Definition: CFitProblem.h:443
CExperiment::storeExtendedTimeSeriesData
void storeExtendedTimeSeriesData(C_FLOAT64 time)
Definition: CExperiment.cpp:576
CCopasiVectorN::add
virtual bool add(const CType &src)
Definition: CCopasiVector.h:621
CCrossValidationSet::getWeight
const C_FLOAT64 & getWeight() const
Definition: CExperimentSet.cpp:512
CFitProblem::mpFisherEigenvaluesMatrixInterface
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpFisherEigenvaluesMatrixInterface
Definition: CFitProblem.h:456
CFitProblem::getScaledFisherInformation
CArrayAnnotation & getScaledFisherInformation() const
Definition: CFitProblem.cpp:1790
CFitProblem::setModel
virtual bool setModel(CModel *pModel)
Definition: CFitProblem.cpp:414
CExperiment::restoreModelIndependentData
bool restoreModelIndependentData()
Definition: CExperiment.cpp:1227
COptProblem::mStoreResults
bool mStoreResults
Definition: COptProblem.h:530
MCFitting
#define MCFitting
Definition: CCopasiMessage.h:54
CFitProblem::mpSteadyState
CSteadyStateTask * mpSteadyState
Definition: CFitProblem.h:312
CFitProblem::mFisherScaledEigenvalues
CMatrix< C_FLOAT64 > mFisherScaledEigenvalues
Definition: CFitProblem.h:473
CFitProblem::getFisherInformationEigenvalues
CArrayAnnotation & getFisherInformationEigenvalues() const
Definition: CFitProblem.cpp:1784
COptProblem::mpParmSubtaskCN
std::string * mpParmSubtaskCN
Definition: COptProblem.h:390
CFitProblem::mpFisherScaledEigenvectorsMatrix
CArrayAnnotation * mpFisherScaledEigenvectorsMatrix
Definition: CFitProblem.h:479
CCopasiParameter::getValue
const Value & getValue() const
Definition: CCopasiParameter.cpp:247
CFitProblem::getFisherInformationEigenvectors
CArrayAnnotation & getFisherInformationEigenvectors() const
Definition: CFitProblem.cpp:1787
CModel::getInitialTime
const C_FLOAT64 & getInitialTime() const
Definition: CModel.cpp:1184
CFitProblem::getCrossValidationSet
const CCrossValidationSet & getCrossValidationSet() const
Definition: CFitProblem.cpp:1805
CCopasiDataModel::getTaskList
CCopasiVectorN< CCopasiTask > * getTaskList()
Definition: CCopasiDataModel.cpp:1415
CExperimentSet::hasDataForTaskType
bool hasDataForTaskType(const CCopasiTask::Type &type) const
Definition: CExperimentSet.cpp:364
CCopasiParameterGroup::setValue
bool setValue(const std::string &name, const CType &value)
Definition: CCopasiParameterGroup.h:361
CExperimentSet::compile
bool compile(const std::vector< CCopasiContainer * > listOfContainer=CCopasiContainer::EmptyList)
Definition: CExperimentSet.cpp:96
COptProblem::mCPUTime
CCopasiTimer mCPUTime
Definition: COptProblem.h:514
COptProblem::setSolution
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
Definition: COptProblem.cpp:682
CFitProblem::mpFisherScaledEigenvaluesMatrixInterface
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpFisherScaledEigenvaluesMatrixInterface
Definition: CFitProblem.h:474
CCopasiParameterGroup::assertGroup
CCopasiParameterGroup * assertGroup(const std::string &name)
Definition: CCopasiParameterGroup.cpp:297
COptProblem::mpParmObjectiveExpression
std::string * mpParmObjectiveExpression
Definition: COptProblem.h:401
CFitProblem::getScaledFisherInformationEigenvalues
CArrayAnnotation & getScaledFisherInformationEigenvalues() const
Definition: CFitProblem.cpp:1793
CModel::getModelParameterSets
const CCopasiVectorN< CModelParameterSet > & getModelParameterSets() const
Definition: CModel.cpp:1066
CExperiment::sumOfSquares
C_FLOAT64 sumOfSquares(const size_t &index, C_FLOAT64 *&residuals) const
Definition: CExperiment.cpp:426
CFitProblem::mpCorrelationMatrixInterface
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpCorrelationMatrixInterface
Definition: CFitProblem.h:485
CCopasiParameterGroup::getParameter
CCopasiParameter * getParameter(const std::string &name)
Definition: CCopasiParameterGroup.cpp:343
CFitProblem::mExperimentUndoMethods
CMatrix< UpdateMethod * > mExperimentUndoMethods
Definition: CFitProblem.h:328
CExperiment::sumOfSquaresStore
C_FLOAT64 sumOfSquaresStore(const size_t &index, C_FLOAT64 *&dependentValues)
Definition: CExperiment.cpp:511
CMatrix::resize
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
CFitProblem::mFisherScaledEigenvectors
CMatrix< C_FLOAT64 > mFisherScaledEigenvectors
Definition: CFitProblem.h:477
CFitProblem::mCrossValidationConstraints
CMatrix< CFitConstraint * > mCrossValidationConstraints
Definition: CFitProblem.h:376
CFitProblem::initialize
virtual bool initialize()
Definition: CFitProblem.cpp:420
CFitProblem::fixBuild55
void fixBuild55()
Definition: CFitProblem.cpp:2108
COptProblem::elevateChildren
virtual bool elevateChildren()
Definition: COptProblem.cpp:167
COptProblem::mpParmMaximize
bool * mpParmMaximize
Definition: COptProblem.h:406
CExperimentSet::keyToIndex
size_t keyToIndex(const std::string &key) const
Definition: CExperimentSet.cpp:389
CFitProblem::setCallBack
virtual bool setCallBack(CProcessReport *pCallBack)
Definition: CFitProblem.cpp:417
CFitProblem::mCrossValidationRMS
C_FLOAT64 mCrossValidationRMS
Definition: CFitProblem.h:397
CFitProblem::createObjectiveFunction
virtual bool createObjectiveFunction()
Definition: CFitProblem.cpp:1281
CVectorCore::size
size_t size() const
Definition: CVector.h:100
CFitProblem::getStdDeviation
const C_FLOAT64 & getStdDeviation() const
Definition: CFitProblem.cpp:1775
CObjectInterface::ObjectSet
std::set< const CObjectInterface * > ObjectSet
Definition: CCopasiObject.h:157
dsyev_
int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a, integer *lda, doublereal *w, doublereal *work, integer *lwork, integer *info)
CCopasiObject::getDirectDependencies
virtual const DataObjectSet & getDirectDependencies(const DataObjectSet &context=DataObjectSet()) const
Definition: CCopasiObject.cpp:304
CExperiment::getExperimentType
const CCopasiTask::Type & getExperimentType() const
Definition: CExperiment.cpp:1237
CCopasiParameter::getObjectDisplayName
virtual std::string getObjectDisplayName(bool regular=true, bool richtext=false) const
Definition: CCopasiParameter.cpp:560
CExperiment::getTimeData
const CVector< C_FLOAT64 > & getTimeData() const
Definition: CExperiment.cpp:1270
CFitItem::getCrossValidationCount
size_t getCrossValidationCount() const
Definition: CFitItem.cpp:257
CCopasiRootContainer::getKeyFactory
static CKeyFactory * getKeyFactory()
Definition: CCopasiRootContainer.cpp:163
CCopasiTimeVariable::getSeconds
C_INT64 getSeconds(const bool &bounded=false) const
Definition: CopasiTime.cpp:126
CFitProblem::checkFunctionalConstraints
virtual bool checkFunctionalConstraints()
Definition: CFitProblem.cpp:787
CAnnotatedMatrix.h
Header file of class CArrayAnnotation.
C_FLOAT64
#define C_FLOAT64
Definition: copasi.h:92
CModel::getInitialState
const CState & getInitialState() const
Definition: CModel.cpp:1768
CFitProblem::mpCrossValidationSet
CCrossValidationSet * mpCrossValidationSet
Definition: CFitProblem.h:355
CSteadyStateTask::initialize
virtual bool initialize(const OutputFlag &of, COutputHandler *pOutputHandler, std::ostream *pOstream)
Definition: CSteadyStateTask.cpp:220
CFitProblem::print
virtual void print(std::ostream *ostream) const
Definition: CFitProblem.cpp:1130
CVectorCore::array
CType * array()
Definition: CVector.h:139
CFitItem::updateBounds
bool updateBounds(std::vector< COptItem * >::iterator it)
Definition: CFitItem.cpp:279
CExperimentSet::getDataPointCount
size_t getDataPointCount() const
Definition: CExperimentSet.cpp:449
CCrossValidationSet::getThreshold
const unsigned C_INT32 & getThreshold() const
Definition: CExperimentSet.cpp:518
CCopasiTask::CDescription::print
virtual void print(std::ostream *ostream) const
Definition: CCopasiTask.cpp:408
CCopasiMessage::getLastMessage
static CCopasiMessage getLastMessage()
Definition: CCopasiMessage.cpp:81
CFitProblem::getCrossValidationSD
const C_FLOAT64 & getCrossValidationSD() const
Definition: CFitProblem.cpp:1825
UTCTimeStamp
std::string UTCTimeStamp()
Definition: utility.cpp:64
CCopasiObject::buildUpdateSequence
static std::vector< Refresh * > buildUpdateSequence(const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())
Definition: CCopasiObject.cpp:438
COptProblem::mWorstValue
C_FLOAT64 mWorstValue
Definition: COptProblem.h:385
CFitProblem::mFisherScaled
CMatrix< C_FLOAT64 > mFisherScaled
Definition: CFitProblem.h:466
CFitProblem::initializeParameter
void initializeParameter()
Definition: CFitProblem.cpp:232
createParameterSetsForExperiment
void createParameterSetsForExperiment(CExperiment *pExp)
Definition: CFitProblem.cpp:807
CFitProblem::mCrossValidationConstraintRefreshes
CVector< std::vector< Refresh * > > mCrossValidationConstraintRefreshes
Definition: CFitProblem.h:382
CModel::updateInitialValues
bool updateInitialValues()
Definition: CModel.cpp:1461
COptProblem::mOriginalVariables
CVector< C_FLOAT64 > mOriginalVariables
Definition: COptProblem.h:484
CCopasiParameterGroup::getValue
const CCopasiParameter::Value & getValue(const std::string &name) const
Definition: CCopasiParameterGroup.cpp:395
CCopasiParameterGroup::assertParameter
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
Definition: CCopasiParameterGroup.h:190
CModel
Definition: CModel.h:50
CTrajectoryTask::restore
virtual bool restore()
Definition: CTrajectoryTask.cpp:497
CFitProblem::mpFisherEigenvectorsMatrixInterface
CCopasiMatrixInterface< CMatrix< C_FLOAT64 > > * mpFisherEigenvectorsMatrixInterface
Definition: CFitProblem.h:460
CCopasiProblem::mpModel
CModel * mpModel
Definition: CCopasiProblem.h:53
CFitProblem::setCreateParameterSets
void setCreateParameterSets(const bool &create)
Definition: CFitProblem.cpp:255
CCopasiParameterGroup::getGroup
CCopasiParameterGroup * getGroup(const std::string &name)
Definition: CCopasiParameterGroup.cpp:383
CFitProblem::mpExperimentSet
CExperimentSet * mpExperimentSet
Definition: CFitProblem.h:306
CFitProblem::mpFisherScaledMatrix
CArrayAnnotation * mpFisherScaledMatrix
Definition: CFitProblem.h:468
CExperiment::printResult
virtual void printResult(std::ostream *ostream) const
Definition: CExperiment.cpp:1409
CMatrix::numCols
virtual size_t numCols() const
Definition: CMatrix.h:144
CModel::buildInitialRefreshSequence
std::vector< Refresh * > buildInitialRefreshSequence(std::set< const CCopasiObject * > &changedObjects)
Definition: CModel.cpp:4164
CExperiment::updateModelWithIndependentData
bool updateModelWithIndependentData(const size_t &index)
Definition: CExperiment.cpp:1217
CCopasiProblem::mpCallBack
CProcessReport * mpCallBack
Definition: CCopasiProblem.h:58
COptItem::getObjectDisplayName
std::string getObjectDisplayName() const
Definition: COptItem.cpp:134
CCopasiObject::setObjectName
bool setObjectName(const std::string &name)
Definition: CCopasiObject.cpp:158
CCopasiTask::setUpdateModel
void setUpdateModel(const bool &updateModel)
Definition: CCopasiTask.cpp:186
CFitProblem::getRMS
const C_FLOAT64 & getRMS() const
Definition: CFitProblem.cpp:1772
CCopasiTimeVariable::getMilliSeconds
C_INT64 getMilliSeconds(const bool &bounded=false) const
Definition: CopasiTime.cpp:118
CCopasiDataModel::ObjectFromName
CCopasiObject * ObjectFromName(const std::vector< CCopasiContainer * > &listOfContainer, const CCopasiObjectName &CN) const
Definition: CCopasiDataModel.cpp:1940
CFitProblem::mCrossValidationUndoMethods
CMatrix< UpdateMethod * > mCrossValidationUndoMethods
Definition: CFitProblem.h:365
CCopasiContainer::addObjectReference
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
Definition: CCopasiContainer.h:89
CSteadyStateTask::restore
virtual bool restore()
Definition: CSteadyStateTask.cpp:422
CFitItem::getExperiments
std::string getExperiments() const
Definition: CFitItem.cpp:210
CFitProblem::mTrajectoryUpdate
bool mTrajectoryUpdate
Definition: CFitProblem.h:496
CMatrix::array
virtual CType * array()
Definition: CMatrix.h:337
CExperimentSet::getExperimentCount
size_t getExperimentCount() const
Definition: CExperimentSet.cpp:324
CFitItem::getCrossValidation
const std::string & getCrossValidation(const size_t &index) const
Definition: CFitItem.cpp:244
CFitProblem::mCrossValidationUpdateMethods
CMatrix< UpdateMethod * > mCrossValidationUpdateMethods
Definition: CFitProblem.h:360
CFitProblem::mCrossValidationSD
C_FLOAT64 mCrossValidationSD
Definition: CFitProblem.h:402
CFitItem::getExperimentCount
size_t getExperimentCount() const
Definition: CFitItem.cpp:207
CCopasiObject::getObjectParent
CCopasiContainer * getObjectParent() const
Definition: CCopasiObject.cpp:276
CFitProblem::mResiduals
CVector< C_FLOAT64 > mResiduals
Definition: CFitProblem.h:428
CFitProblem::mExperimentConstraintRefreshes
CVector< std::vector< Refresh * > > mExperimentConstraintRefreshes
Definition: CFitProblem.h:345
CCopasiMessage::getHighestSeverity
static CCopasiMessage::Type getHighestSeverity()
Definition: CCopasiMessage.cpp:119
Generated on Thu Jul 2 2015 14:50:13 for COPASI API by   doxygen 1.8.6