COPASI API  4.16.103
CMCAMethod.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) 2004 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include <cmath>
16 #include <limits>
17 
18 #include "copasi.h"
19 #include "CMCAMethod.h"
20 #include "CMCAProblem.h"
21 #include "CSteadyStateTask.h"
22 #include "model/CModel.h"
23 #include "utilities/CReadConfig.h"
24 #include "utilities/utility.h"
25 #include "utilities/CLinkMatrix.h"
26 
27 #include "lapack/blaswrap.h"
28 #include "lapack/lapackwrap.h"
29 
30 //TODO: put all matrix resizing and annotations creation in one place, so
31 // that it has to be done only once if several MCA are calculated (e.g. in a scan)
32 
33 // static
35 {
36  return new CMCAMethod();
37 }
38 
39 /**
40  * Default constructor
41  */
43  CCopasiMethod(CCopasiTask::mca, CCopasiMethod::mcaMethodReder, pParent),
44  mpModel(NULL),
45  mpUseReeder(NULL),
46  mpUseSmallbone(NULL),
47  mFactor(1.0e-9),
48  mSteadyStateResolution(1.0e-9),
49  mSSStatus(CSteadyStateMethod::notFound),
50  mpSteadyStateTask(NULL),
51  mLinkZero(),
52  mReducedStoichiometry(),
53  mElasticityDependencies()
54 {
56  initObjects();
57 }
58 
60  const CCopasiContainer * pParent):
61  CCopasiMethod(src, pParent),
62  mpModel(NULL),
63  mpUseReeder(NULL),
64  mpUseSmallbone(NULL),
65  mFactor(src.mFactor),
66  mSteadyStateResolution(src.mSteadyStateResolution),
67  mSSStatus(CSteadyStateMethod::notFound),
68  mpSteadyStateTask(NULL),
69  mLinkZero(src.mLinkZero),
70  mReducedStoichiometry(src.mReducedStoichiometry),
71  mElasticityDependencies(src.mElasticityDependencies)
72 {
74  initObjects();
75 }
76 
78 {
80  tmp = new CArrayAnnotation("Unscaled elasticities", this,
83  tmp->setDescription("Unscaled elasticity matrix");
84  tmp->setDimensionDescription(0, "Reactions (reduced system)");
85  tmp->setDimensionDescription(1, "Species (reduced system)");
87 
88  tmp = new CArrayAnnotation("Unscaled concentration control coefficients", this,
91  tmp->setDescription("Unscaled concentration control coefficients");
92  tmp->setDimensionDescription(0, "Species (reduced system)");
93  tmp->setDimensionDescription(1, "Reactions (reduced system)");
94  mUnscaledConcCCAnn = tmp;
95 
96  tmp = new CArrayAnnotation("Unscaled flux control coefficients", this,
99  tmp->setDescription("Unscaled flux control coefficients");
100  tmp->setDimensionDescription(0, "Reactions (reduced system)");
101  tmp->setDimensionDescription(1, "Reactions (reduced system)");
102  mUnscaledFluxCCAnn = tmp;
103 
104  tmp = new CArrayAnnotation("Scaled elasticities", this,
107  tmp->setDescription("Scaled elasticity matrix");
108  tmp->setDimensionDescription(0, "Reactions (reduced system)");
109  tmp->setDimensionDescription(1, "Species (reduced system)");
111 
112  tmp = new CArrayAnnotation("Scaled concentration control coefficients", this,
115  tmp->setDescription("Scaled concentration control coefficients");
116  tmp->setDimensionDescription(0, "Species (reduced system)");
117  tmp->setDimensionDescription(1, "Reactions (reduced system)");
118  mScaledConcCCAnn = tmp;
119 
120  tmp = new CArrayAnnotation("Scaled flux control coefficients", this,
123  tmp->setDescription("Scaled flux control coefficients");
124  tmp->setDimensionDescription(0, "Reactions (reduced system)");
125  tmp->setDimensionDescription(1, "Reactions (reduced system)");
126  mScaledFluxCCAnn = tmp;
127 }
128 
129 /**
130  * Deconstructor
131  */
133 {
135  //delSsipvt();
136 }
137 
139 {
140  CCopasiParameter *pParm;
141 
142  assertParameter("Modulation Factor", CCopasiParameter::UDOUBLE, 1.0e-009);
143 
144  if ((pParm = getParameter("MCA.ModulationFactor")) != NULL)
145  {
146  setValue("Modulation Factor", *pParm->getValue().pUDOUBLE);
147  removeParameter("MCA.ModulationFactor");
148  }
149 
152 }
153 
155 {
157  return true;
158 }
159 
161 {
162  if (!mpModel) return;
163 
169 
171  mpModel->getTotSteps());
175 
180 
185 
186  // Reactions are columns, species are rows
191 
192  // Reactions are columns and rows
197 
199 }
200 
201 //this calculates the elasticities as d(particle flux)/d(particle number)
202 //which is the same as d(flux of substance)/d(amount of substance)
204 {
205  assert(mpModel);
206 
207  // We need the number of metabolites determined by reactions.
208  size_t numMetabs =
210 
212  CCopasiVector< CMetab >::iterator endSpecies = itSpecies + numMetabs;
213 
214  // Calculate the dependencies of the elasticities this is helpful for scaling to determine
215  // whether 0/0 is 0 or NaN
216 
217  bool * pElasticityDependency = mElasticityDependencies.array();
220 
221  for (; itReac != endReac; ++itReac)
222  {
223  CCopasiObject::DataObjectSet ReactionDependencies;
224  (*itReac)->getParticleFluxReference()->getAllDependencies(ReactionDependencies, DataObjectSet());
225 
226  itSpecies = mpModel->getMetabolitesX().begin();
227 
228  for (; itSpecies != endSpecies; ++itSpecies, ++pElasticityDependency)
229  {
230  *pElasticityDependency = (ReactionDependencies.find((*itSpecies)->getValueReference()) != ReactionDependencies.end());
231  }
232  }
233 
234  const CVector< C_FLOAT64 > & ParticleFlux = mpModel->getParticleFlux();
235 
236  // mUnscaledElasticities.resize(numReacs, numMetabs);
237  C_FLOAT64 * pElasticity;
238 
240 
241  C_FLOAT64 Store, InvDelta;
242 
243  C_FLOAT64 X1, X2;
244 
245  // Arrays to store function value
246  size_t numReacs = mpModel->getReactions().size();
247 
248  CVector< C_FLOAT64 > Y1(numReacs);
249 
250  C_FLOAT64 * pY1;
251 
252  CVector< C_FLOAT64 > Y2(numReacs);
253 
254  C_FLOAT64 * pY2;
255 
256  // calculate elasticities
257  itSpecies = mpModel->getMetabolitesX().begin();
258 
259  for (size_t j = 0; itSpecies != endSpecies; ++itSpecies, ++j)
260  {
261  Store = (*itSpecies)->getValue();
262 
263  // We only need to make sure that we do not have an underflow problem
264  if (fabs(Store) < 100 * std::numeric_limits< C_FLOAT64 >::min())
265  {
266  X1 = 0.0;
267 
268  if (Store < 0.0)
270  else
272  }
273  else
274  {
275  X1 = Store * (1.0 + mFactor);
276  X2 = Store * (1.0 - mFactor);
277  }
278 
279  InvDelta = 1.0 / (X1 - X2);
280 
281  // let's take X+dx
282  (*itSpecies)->setValue(X1);
283  mpModel->updateSimulatedValues(false); // TODO test if true or false should be used.
284 
285  // get the fluxes
286  Y1 = ParticleFlux;
287 
288  // now X-dx
289  (*itSpecies)->setValue(X2);
290  mpModel->updateSimulatedValues(false); // TODO test if true or false should be used.
291 
292  // get the fluxes
293  Y2 = ParticleFlux;
294 
295  // set column j of Dxv
296  pElasticity = mUnscaledElasticities.array() + j;
297  pY1 = Y1.array();
298  pY2 = Y2.array();
299 
300  for (; pElasticity < pElasticityEnd; pElasticity += numMetabs, ++pY1, ++pY2)
301  * pElasticity = (*pY1 - *pY2) * InvDelta;
302 
303  // restore the value of the species
304  (*itSpecies)->setValue(Store);
305  }
306 
307  // make sure the fluxes are correct afterwards (needed for scaling of the MCA results)
309 }
310 
312 {
313  assert(mpModel);
314 
315  // TODO CRITICAL We must not use the reduced stoichiometry matrix
316  // Instead use N * mLinkZero
317 
318  // Calculate RedStoi * mUnscaledElasticities;
319  // Note the columns of mUnscaledElasticities must be reordered
320 
322 
323  // Initialize the unscaled concentration control coefficients to 0.0
324  mUnscaledConcCC = 0.0;
325 
326  // aux1 := mUnscaledElasticities * L
327  CMatrix<C_FLOAT64> aux1;
329 
330  // We can now undo the column pivoting
332 
333  assert(mReducedStoichiometry.numCols() == aux1.numRows());
334 
335  // aux2 := RedStoi * aux1
336  // DGEMM (TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
337  // C := alpha A B + beta C
339 
340  char TRANSA = 'N';
341  char TRANSB = 'N';
342  C_INT M = (C_INT) aux2.numCols(); /* LDA, LDC */
343  C_INT N = (C_INT) aux2.numRows();
345  C_FLOAT64 Alpha = 1.0;
346  C_INT LDA = (C_INT) aux1.numCols();
348  C_FLOAT64 Beta = 0.0;
349  C_INT LDC = (C_INT) aux2.numCols();
350 
351  dgemm_(&TRANSA, &TRANSB, &M, &N, &K, &Alpha, aux1.array(), &LDA,
352  mReducedStoichiometry.array(), &LDB, &Beta, aux2.array(), &LDC);
353 
354  // Invert aux2
355  C_INT info;
356  CVector<C_INT> Ipiv(M);
357 
358  // LU decomposition of aux2 (for inversion)
359  dgetrf_(&M, &M, aux2.array(), &M, Ipiv.array(), &info);
360 
361  if (info != 0) return false;
362 
363  C_INT lwork = -1; // Instruct dgetri_ to determine work array size.
364  CVector< C_FLOAT64 > work(1);
365 
366  dgetri_(&M, aux2.array(), &M, Ipiv.array(), work.array(), &lwork, &info);
367 
368  lwork = (C_INT) work[0];
369  work.resize(lwork);
370 
371  // now invert aux2 (result in aux2)
372  dgetri_(&M, aux2.array(), &M, Ipiv.array(), work.array(), &lwork, &info);
373 
374  if (info != 0) return false;
375 
376  // aux1 := -1.0 * aux2 * RedStoi
377 
378  aux1.resize(aux2.numRows(), mReducedStoichiometry.numCols());
379 
380  M = (C_INT) aux1.numCols();
381  N = (C_INT) aux1.numRows();
382  K = (C_INT) aux2.numCols();
383  Alpha = -1.0;
385  LDB = (C_INT) aux2.numCols();
386  Beta = 0.0;
387  LDC = (C_INT) aux1.numCols();
388 
389  // DGEMM (TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
390  // C := alpha A B + beta C
391  dgemm_(&TRANSA, &TRANSB, &M, &N, &K, &Alpha, mReducedStoichiometry.array(), &LDA,
392  aux2.array(), &LDB, &Beta, aux1.array(), &LDC);
393 
394  // mUnscaledConcCC := L * aux1
396 
397  // We need to swap the rows since they are with respect to reordered stoichiometry .
399 
400  return true;
401 }
402 
403 bool CMCAMethod::calculateUnscaledFluxCC(const bool & status)
404 {
405  assert(mpModel);
406  //size_t i, j, k;
407 
408  // mUnscaledFluxCC := I + mUnscaledElasticities * mUnscaledConcCC
409 
410  char UPLO = 'A';
412  C_FLOAT64 Alpha = 0.0;
413  C_FLOAT64 Beta = 1.0;
414 
415  // Initialize mUnscaledFluxCC to the identity matrix;
416  dlaset_(&UPLO, &M, &M, &Alpha, &Beta, mUnscaledFluxCC.array(), &M);
417 
418  if (status)
419  {
420  char TRANSA = 'N';
421  char TRANSB = 'N';
422  M = (C_INT) mUnscaledFluxCC.numCols(); /* LDA, LDC */
425 
429 
430  Alpha = 1.0;
431  Beta = 1.0;
432 
433  dgemm_(&TRANSA, &TRANSB, &M, &N, &K, &Alpha, mUnscaledConcCC.array(), &LDA,
434  mUnscaledElasticities.array(), &LDB, &Beta, mUnscaledFluxCC.array(), &LDC);
435  }
436 
437  return status;
438 }
439 
440 bool CMCAMethod::scaleMCA(const bool & status, C_FLOAT64 res)
441 {
442  assert(mpModel);
443  // The number of metabs determined by reaction.
444  size_t numSpeciesReaction =
447  CCopasiVector< CMetab >::const_iterator itSpecies = metabs.begin();
448  CCopasiVector< CMetab >::const_iterator endSpecies = itSpecies + numSpeciesReaction;
449 
452  CCopasiVector< CReaction >::const_iterator endReaction = reacs.end();
453 
454  size_t col;
455 
456  // Scale Elasticities
457  C_FLOAT64 * pUnscaled;
458  C_FLOAT64 * pScaled;
459  bool * pElasticityDependency;
460 
461  // Reactions are rows, species are columns
462  for (col = 0; itSpecies != endSpecies; ++itSpecies, col++)
463  {
464  C_FLOAT64 VolumeInv = 1.0 / (*itSpecies)->getCompartment()->getValue();
465  C_FLOAT64 Number = (*itSpecies)->getValue();
466 
467  for (itReaction = reacs.begin(),
468  pUnscaled = mUnscaledElasticities.array() + col,
469  pScaled = mScaledElasticities.array() + col,
470  pElasticityDependency = mElasticityDependencies.array() + col;
471  itReaction != endReaction;
472  ++itReaction,
473  pUnscaled += numSpeciesReaction,
474  pScaled += numSpeciesReaction,
475  pElasticityDependency += numSpeciesReaction)
476  {
477  if (!*pElasticityDependency)
478  {
479  *pScaled = 0.0;
480  }
481  else if (fabs((*itReaction)->getFlux() * VolumeInv) >= res)
482  {
483  *pScaled = *pUnscaled * Number / (*itReaction)->getParticleFlux();
484  }
485  else
486  {
487  *pScaled = (((*itReaction)->getFlux() < 0.0) ? - std::numeric_limits<C_FLOAT64>::infinity() : std::numeric_limits<C_FLOAT64>::infinity());
488  }
489  }
490  }
491 
492  //if we are not in a steady state we cannot calculate CCs
494  !status)
495  {
496  mScaledConcCC = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
497  mScaledFluxCC = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
498 
499  return false;
500  }
501 
502  // Scale ConcCC
503  // Reactions are columns, species are rows
504  itSpecies = metabs.begin();
505  pUnscaled = mUnscaledConcCC.array();
506  pScaled = mScaledConcCC.array();
507 
508  for (; itSpecies != endSpecies; ++itSpecies)
509  {
510  // In rare occasions the concentration might not be updated
511  (*(*itSpecies)->getConcentrationReference()->getRefresh())();
512  C_FLOAT64 alt = fabs((*itSpecies)->getConcentration());
513 
514  for (itReaction = reacs.begin(); itReaction != endReaction; ++itReaction, ++pUnscaled, ++pScaled)
515  {
516  if (alt >= res)
517  *pScaled =
518  *pUnscaled * (*itReaction)->getParticleFlux() / (*itSpecies)->getValue();
519  else
520  *pScaled = *pUnscaled * std::numeric_limits<C_FLOAT64>::infinity();
521  }
522  }
523 
525  itReaction = reacs.begin();
526  pUnscaled = mUnscaledFluxCC.array();
527  pScaled = mScaledFluxCC.array();
528 
529  for (; itReaction != endReaction; ++itReaction)
530  {
531  const CCompartment * pCompartment = (*itReaction)->getLargestCompartment();
532  C_FLOAT64 Resolution = res * pCompartment->getValue();
533 
534  C_FLOAT64 Scale = (*itReaction)->getFlux();
535  C_FLOAT64 tmp = 0.0;
536  C_FLOAT64 eq = 0.0;
537 
538  // We use the summation theorem to verify the scaling
539  C_FLOAT64 *pSum = pUnscaled;
540  C_FLOAT64 *pSumEnd = pSum + mScaledFluxCC.numCols();
541 
542  for (itReactionCol = reacs.begin(); pSum != pSumEnd; ++itReactionCol, ++pSum)
543  {
544  tmp += *pSum * (*itReactionCol)->getFlux();
545  eq += fabs(*pSum);
546  }
547 
548  eq /= mScaledFluxCC.numCols();
549 
550  if (fabs(tmp) < Resolution && eq >= Resolution)
551  {
552  Scale = std::numeric_limits< C_FLOAT64 >::infinity();
553  }
554 
555  tmp = 0.0;
556 
557  for (itReactionCol = reacs.begin(); itReactionCol != endReaction; ++itReactionCol, ++pUnscaled, ++pScaled)
558  {
559  // In the diagonal the scaling factors cancel.
560  if (eq < res)
561  {
562  *pScaled = (itReactionCol != itReaction) ? 0.0 : 1.0;
563  }
564  else if (itReactionCol == itReaction)
565  {
566  *pScaled = *pUnscaled;
567  }
568  else if (fabs(Scale) >= res)
569  {
570  *pScaled = *pUnscaled * (*itReactionCol)->getFlux() / Scale;
571  }
572  else
573  {
574  const CCompartment * pCompartmentCol = (*itReactionCol)->getLargestCompartment();
575  C_FLOAT64 ScaleCol = (pCompartmentCol == NULL) ?
576  fabs((*itReactionCol)->getFlux()) :
577  fabs((*itReactionCol)->getFlux() / pCompartmentCol->getValue());
578 
579  if (fabs(ScaleCol) <= res)
580  {
581  *pScaled = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
582  }
583  else
584  {
585  *pScaled = (((*itReaction)->getFlux() < 0.0) ?
586  - std::numeric_limits<C_FLOAT64>::infinity() :
587  std::numeric_limits<C_FLOAT64>::infinity());
588  }
589  }
590 
591  tmp += *pScaled;
592  }
593  }
594 
595  return true;
596 }
597 
599 {
600  bool success = true;
601 
602  C_FLOAT64 * pScaled = mScaledConcCC.array();
603  C_FLOAT64 * pScaledRowEnd = pScaled + mScaledConcCC.numCols();
604  C_FLOAT64 * pScaledEnd = pScaled + mScaledConcCC.size();
605 
609  Sum = 0.0;
610  Max = 0.0;
611  C_FLOAT64 * pSum = Sum.array();
612  C_FLOAT64 * pMax = Max.array();
613  C_FLOAT64 * pError = Error.array();
614 
615  for (; pScaled != pScaledEnd; pScaledRowEnd += mScaledConcCC.numCols(), ++pSum, ++pMax, ++pError)
616  {
617  for (; pScaled != pScaledRowEnd; ++pScaled)
618  {
619  success &= !isnan(*pScaled);
620 
621  *pSum += *pScaled;
622  *pMax = std::max(*pMax, fabs(*pScaled));
623  }
624 
625  *pError = (*pMax > 100.0 * std::numeric_limits< C_FLOAT64 >::min()) ? fabs(*pSum) / *pMax : 0.0;
626  success &= *pError < resolution;
627  }
628 
629  pScaled = mScaledFluxCC.array();
630  pScaledRowEnd = pScaled + mScaledFluxCC.numCols();
631  pScaledEnd = pScaled + mScaledFluxCC.size();
632 
634  Max.resize(mScaledFluxCC.numRows());
635  Error.resize(mScaledFluxCC.numRows());
636  Sum = 0.0;
637  Max = 0.0;
638 
639  pSum = Sum.array();
640  pMax = Max.array();
641  pError = Error.array();
642 
643  for (; pScaled != pScaledEnd; pScaledRowEnd += mScaledConcCC.numCols(), ++pSum, ++pMax, ++pError)
644  {
645  for (; pScaled != pScaledRowEnd; ++pScaled)
646  {
647  success &= !isnan(*pScaled);
648 
649  *pSum += *pScaled;
650  *pMax = std::max(*pMax, fabs(*pScaled));
651  }
652 
653  *pError = (*pMax > 100.0 * std::numeric_limits< C_FLOAT64 >::min()) ? fabs(1.0 - *pSum) / *pMax : 0.0;
654  success &= *pError < resolution;
655  }
656 
657  return success;
658 }
659 
660 /**
661  * Set the Model
662  */
664 {
665  mpModel = model;
666 }
667 
668 /**
669  * the steady state MCA entry point
670  * @param ss_solution refer to steady-state solution
671  * @param refer to the resolution
672  */
674 {
675  bool success = true;
676  bool SummationTheoremsOK = false;
677 
678  assert(mpModel);
679 
681 
683  {
684  if (*mpUseReeder)
685  {
688  success &= calculateUnscaledFluxCC(success);
689  success &= scaleMCA(success, res);
690  SummationTheoremsOK = checkSummationTheorems(res);
691  }
692 
693  if (*mpUseSmallbone && !SummationTheoremsOK)
694  {
695  success = true;
696 
697  createLinkMatrix(true);
699  success &= calculateUnscaledFluxCC(success);
700  success &= scaleMCA(success, res);
701  SummationTheoremsOK = checkSummationTheorems(res);
702  }
703 
704  if (!SummationTheoremsOK)
705  {
707  }
708  }
709  else
710  {
711  mUnscaledConcCC = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
712  mUnscaledFluxCC = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
713  }
714 
715  return success;
716 }
717 
718 bool CMCAMethod::createLinkMatrix(const bool & useSmallbone)
719 {
720  if (mpModel == NULL ||
721  mpSteadyStateTask == NULL)
722  {
723  return false;
724  }
725 
726  if (useSmallbone)
727  {
732  }
733  else
734  {
735  mLinkZero = mpModel->getL0();
737  }
738 
739  return true;
740 }
741 
742 /**
743  * Read SSMCAUnscaled from configuration file
744  */
746 {
747  C_INT32 Fail = 0;
748 
749  if ((Fail = configBuffer.getVariable("SSMCAUnscaled", "C_INT16",
750  (void *) & mSSReder,
752  return Fail;
753 
754  return Fail;
755 }
756 
758 {
759  // check if current state is a steady state
760  // if not, calculate TimeMCA only
761  if (1 /*mIsSteadyState*/)
762  {
764  }
765  else
766  {
767  // CalculateTimeMCA(mSteadyStateResolution);
768  }
769 
770  return true;
771 }
772 
774 {
775  mpSteadyStateTask = pSteadyStateTask;
776 
777  if (mpSteadyStateTask != NULL)
778  {
780  }
781  else
782  {
784  }
785 }
786 
788 {
789  this->mFactor = factor;
790 }
791 
793 {
794  this->mSteadyStateResolution = resolution;
795 }
796 
798 {
799  return this->mpModel;
800 }
801 
802 //virtual
804 {
805  if (!CCopasiMethod::isValidProblem(pProblem)) return false;
806 
807  const CMCAProblem * pP = dynamic_cast<const CMCAProblem *>(pProblem);
808 
809  if (!pP)
810  {
811  //not a TrajectoryProblem
812  CCopasiMessage(CCopasiMessage::ERROR, "Problem is not an MCA problem.");
813  return false;
814  }
815 
816  CModel * pModel = pP->getModel();
817 
818  if (pModel == NULL)
819  return false;
820 
821  // Check if the model contains an ODE.
822  size_t NumODE =
824 
825  if (pModel->getNumIndependentReactionMetabs() < NumODE)
826  {
827  CCopasiMessage(CCopasiMessage::ERROR, "MCA is not applicable for a system with explicit ODEs.");
828  return false;
829  }
830 
831  // Check if the model has a compartment with an assignment
834 
835  for (; it != end; ++it)
836  if ((*it)->getStatus() == CModelEntity::ASSIGNMENT)
837  {
838  CCopasiMessage(CCopasiMessage::ERROR, "MCA is not applicable for a system with changing volumes.");
839  return false;
840  }
841 
842  if (!*mpUseReeder && !*mpUseSmallbone)
843  {
844  CCopasiMessage(CCopasiMessage::ERROR, "At least one of the algorithm Reeder or Smallbone must be selected.");
845  return false;
846  }
847 
848  return true;
849 }
int dgetrf_(integer *m, integer *n, doublereal *a, integer *lda, integer *ipiv, integer *info)
CArrayAnnotation * mScaledElasticitiesAnn
Definition: CMCAMethod.h:191
#define C_INT
Definition: copasi.h:115
bool build(const CMatrix< C_FLOAT64 > &matrix, size_t maxRank=C_INVALID_INDEX)
Definition: CLinkMatrix.cpp:44
const CSteadyStateMethod::ReturnCode & getResult() const
CMatrix< C_FLOAT64 > mUnscaledConcCC
Definition: CMCAMethod.h:184
bool CalculateMCA(C_FLOAT64 res)
Definition: CMCAMethod.cpp:673
CSteadyStateMethod::ReturnCode mSSStatus
Definition: CMCAMethod.h:217
#define K
void setSteadyStateTask(CSteadyStateTask *pSteadyStateTask)
Definition: CMCAMethod.cpp:773
CMatrix< C_FLOAT64 > mScaledElasticities
Definition: CMCAMethod.h:190
int dlaset_(char *uplo, integer *m, integer *n, doublereal *alpha, doublereal *beta, doublereal *a, integer *lda)
virtual size_t size() const
CSteadyStateTask * mpSteadyStateTask
Definition: CMCAMethod.h:219
CLinkMatrix mLinkZero
Definition: CMCAMethod.h:221
CArrayAnnotation * mScaledConcCCAnn
Definition: CMCAMethod.h:194
#define MCMCA
const CVector< C_FLOAT64 > & getParticleFlux() const
Definition: CModel.cpp:1045
iterator begin()
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
virtual size_t numRows() const
Definition: CMatrix.h:138
CMatrix< C_FLOAT64 > mUnscaledFluxCC
Definition: CMCAMethod.h:187
bool leftMultiply(const CMatrix< C_FLOAT64 > &M, CMatrix< C_FLOAT64 > &P) const
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
void setDescription(const std::string &s)
void initObjects()
Definition: CMCAMethod.cpp:77
CArrayAnnotation * mScaledFluxCCAnn
Definition: CMCAMethod.h:197
#define C_INT32
Definition: copasi.h:90
virtual void resizeAllMatrices()
Definition: CMCAMethod.cpp:160
virtual ~CMCAMethod()
Definition: CMCAMethod.cpp:132
const size_t & getNumIndependent() const
virtual bool elevateChildren()
Definition: CMCAMethod.cpp:154
const CMatrix< C_FLOAT64 > & getRedStoi() const
Definition: CModel.cpp:1154
bool scaleMCA(const bool &status, C_FLOAT64 res)
Definition: CMCAMethod.cpp:440
bool rightMultiply(const C_FLOAT64 &alpha, const CMatrix< C_FLOAT64 > &M, CMatrix< C_FLOAT64 > &P) const
size_t getTotSteps() const
Definition: CModel.cpp:1136
void setSteadyStateResolution(C_FLOAT64 factor)
Definition: CMCAMethod.cpp:792
const CModel * getModel() const
Definition: CMCAMethod.cpp:797
virtual bool isValidProblem(const CCopasiProblem *pProblem)
Definition: CMCAMethod.cpp:803
virtual bool isValidProblem(const CCopasiProblem *pProblem)
bool checkSummationTheorems(const C_FLOAT64 &resolution)
Definition: CMCAMethod.cpp:598
bool removeParameter(const std::string &name)
#define DESTRUCTOR_TRACE
Definition: copasi.h:206
iterator end()
std::vector< CType * >::const_iterator const_iterator
Definition: CCopasiVector.h:57
CMatrix< C_FLOAT64 > mScaledConcCC
Definition: CMCAMethod.h:193
bool doRowPivot(CMatrix< C_FLOAT64 > &matrix) const
size_t getNumDependentReactionMetabs() const
Definition: CModel.cpp:1133
bool createLinkMatrix(const bool &useSmallbone=false)
Definition: CMCAMethod.cpp:718
virtual bool process()
Definition: CMCAMethod.cpp:757
const Value & getValue() const
void setFactor(C_FLOAT64 factor)
Definition: CMCAMethod.cpp:787
bool setValue(const std::string &name, const CType &value)
void setCopasiVector(size_t d, const CCopasiContainer *v)
CModelEntity ** beginIndependent()
Definition: CState.cpp:208
CCopasiParameter * getParameter(const std::string &name)
CMatrix< C_FLOAT64 > mReducedStoichiometry
Definition: CMCAMethod.h:223
int dgemm_(char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc)
bool * mpUseSmallbone
Definition: CMCAMethod.h:176
CMatrix< bool > mElasticityDependencies
Definition: CMCAMethod.h:225
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
size_t getNumIndependentReactionMetabs() const
Definition: CModel.cpp:1130
bool doColumnPivot(CMatrix< C_FLOAT64 > &matrix) const
CArrayAnnotation * mUnscaledElasticitiesAnn
Definition: CMCAMethod.h:182
void calculateUnscaledElasticities(C_FLOAT64 res)
Definition: CMCAMethod.cpp:203
void setMode(size_t d, Mode m)
CCopasiVectorNS< CCompartment > & getCompartments()
Definition: CModel.cpp:1145
void setDimensionDescription(size_t d, const std::string &s)
CMatrix< C_FLOAT64 > mUnscaledElasticities
Definition: CMCAMethod.h:181
const C_FLOAT64 & getValue() const
void setModel(CModel *model)
Definition: CMCAMethod.cpp:663
#define C_FLOAT64
Definition: copasi.h:92
CArrayAnnotation * mUnscaledConcCCAnn
Definition: CMCAMethod.h:185
CType * array()
Definition: CVector.h:139
const CStateTemplate & getStateTemplate() const
Definition: CModel.cpp:1172
CMCAMethod(const CCopasiContainer *pParent=NULL)
Definition: CMCAMethod.cpp:42
C_INT32 load(CReadConfig &configBuffer)
Definition: CMCAMethod.cpp:745
virtual size_t size() const
Definition: CMatrix.h:132
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
bool undoColumnPivot(CMatrix< C_FLOAT64 > &matrix) const
const CLinkMatrix & getL0() const
Definition: CModel.cpp:1169
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
static CMCAMethod * createMethod(CCopasiMethod::SubType subType=CCopasiMethod::mcaMethodReder)
Definition: CMCAMethod.cpp:34
Definition: CModel.h:50
C_INT32 getVariable(const std::string &name, const std::string &type, void *pout, CReadConfig::Mode mode=CReadConfig::NEXT)
Definition: CReadConfig.cpp:81
CModel * mpModel
Definition: CMCAMethod.h:172
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
const CMatrix< C_FLOAT64 > & getStoi() const
Definition: CModel.cpp:1160
int dgetri_(integer *n, doublereal *a, integer *lda, integer *ipiv, doublereal *work, integer *lwork, integer *info)
virtual size_t numCols() const
Definition: CMatrix.h:144
CArrayAnnotation * mUnscaledFluxCCAnn
Definition: CMCAMethod.h:188
CMatrix< C_FLOAT64 > mScaledFluxCC
Definition: CMCAMethod.h:196
bool undoRowPivot(CMatrix< C_FLOAT64 > &matrix) const
void initializeParameter()
Definition: CMCAMethod.cpp:138
std::set< const CCopasiObject * > DataObjectSet
bool * mpUseReeder
Definition: CMCAMethod.h:174
C_FLOAT64 mSteadyStateResolution
Definition: CMCAMethod.h:215
const CMatrix< C_FLOAT64 > & getJacobian() const
virtual CType * array()
Definition: CMatrix.h:337
CModel * getModel() const
CModelEntity ** endIndependent()
Definition: CState.cpp:209
#define min(a, b)
Definition: f2c.h:175
bool calculateUnscaledConcentrationCC()
Definition: CMCAMethod.cpp:311
C_INT16 mSSReder
Definition: CMCAMethod.h:202
C_FLOAT64 mFactor
Definition: CMCAMethod.h:213
bool calculateUnscaledFluxCC(const bool &status)
Definition: CMCAMethod.cpp:403
#define max(a, b)
Definition: f2c.h:176