48 mSteadyStateResolution(1.0e-9),
50 mpSteadyStateTask(NULL),
52 mReducedStoichiometry(),
53 mElasticityDependencies()
66 mSteadyStateResolution(src.mSteadyStateResolution),
68 mpSteadyStateTask(NULL),
69 mLinkZero(src.mLinkZero),
70 mReducedStoichiometry(src.mReducedStoichiometry),
71 mElasticityDependencies(src.mElasticityDependencies)
88 tmp =
new CArrayAnnotation(
"Unscaled concentration control coefficients",
this,
91 tmp->
setDescription(
"Unscaled concentration control coefficients");
112 tmp =
new CArrayAnnotation(
"Scaled concentration control coefficients",
this,
144 if ((pParm =
getParameter(
"MCA.ModulationFactor")) != NULL)
221 for (; itReac != endReac; ++itReac)
224 (*itReac)->getParticleFluxReference()->getAllDependencies(ReactionDependencies,
DataObjectSet());
228 for (; itSpecies != endSpecies; ++itSpecies, ++pElasticityDependency)
230 *pElasticityDependency = (ReactionDependencies.find((*itSpecies)->getValueReference()) != ReactionDependencies.end());
259 for (
size_t j = 0; itSpecies != endSpecies; ++itSpecies, ++j)
261 Store = (*itSpecies)->getValue();
279 InvDelta = 1.0 / (X1 - X2);
282 (*itSpecies)->setValue(X1);
289 (*itSpecies)->setValue(X2);
300 for (; pElasticity < pElasticityEnd; pElasticity += numMetabs, ++pY1, ++pY2)
301 * pElasticity = (*pY1 - *pY2) * InvDelta;
304 (*itSpecies)->setValue(Store);
351 dgemm_(&TRANSA, &TRANSB, &M, &N, &K, &Alpha, aux1.
array(), &LDA,
361 if (info != 0)
return false;
368 lwork = (
C_INT) work[0];
374 if (info != 0)
return false;
382 K = (
C_INT) aux2.numCols();
385 LDB = (
C_INT) aux2.numCols();
392 aux2.array(), &LDB, &Beta, aux1.
array(), &LDC);
444 size_t numSpeciesReaction =
459 bool * pElasticityDependency;
462 for (col = 0; itSpecies != endSpecies; ++itSpecies, col++)
464 C_FLOAT64 VolumeInv = 1.0 / (*itSpecies)->getCompartment()->getValue();
465 C_FLOAT64 Number = (*itSpecies)->getValue();
467 for (itReaction = reacs.
begin(),
471 itReaction != endReaction;
473 pUnscaled += numSpeciesReaction,
474 pScaled += numSpeciesReaction,
475 pElasticityDependency += numSpeciesReaction)
477 if (!*pElasticityDependency)
481 else if (fabs((*itReaction)->getFlux() * VolumeInv) >= res)
483 *pScaled = *pUnscaled * Number / (*itReaction)->getParticleFlux();
487 *pScaled = (((*itReaction)->getFlux() < 0.0) ? - std::numeric_limits<C_FLOAT64>::infinity() : std::numeric_limits<C_FLOAT64>::infinity());
496 mScaledConcCC = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
497 mScaledFluxCC = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
504 itSpecies = metabs.
begin();
508 for (; itSpecies != endSpecies; ++itSpecies)
511 (*(*itSpecies)->getConcentrationReference()->getRefresh())();
512 C_FLOAT64 alt = fabs((*itSpecies)->getConcentration());
514 for (itReaction = reacs.
begin(); itReaction != endReaction; ++itReaction, ++pUnscaled, ++pScaled)
518 *pUnscaled * (*itReaction)->getParticleFlux() / (*itSpecies)->getValue();
520 *pScaled = *pUnscaled * std::numeric_limits<C_FLOAT64>::infinity();
525 itReaction = reacs.
begin();
529 for (; itReaction != endReaction; ++itReaction)
531 const CCompartment * pCompartment = (*itReaction)->getLargestCompartment();
534 C_FLOAT64 Scale = (*itReaction)->getFlux();
542 for (itReactionCol = reacs.
begin(); pSum != pSumEnd; ++itReactionCol, ++pSum)
544 tmp += *pSum * (*itReactionCol)->getFlux();
550 if (fabs(tmp) < Resolution && eq >= Resolution)
552 Scale = std::numeric_limits< C_FLOAT64 >::infinity();
557 for (itReactionCol = reacs.
begin(); itReactionCol != endReaction; ++itReactionCol, ++pUnscaled, ++pScaled)
562 *pScaled = (itReactionCol != itReaction) ? 0.0 : 1.0;
564 else if (itReactionCol == itReaction)
566 *pScaled = *pUnscaled;
568 else if (fabs(Scale) >= res)
570 *pScaled = *pUnscaled * (*itReactionCol)->getFlux() / Scale;
574 const CCompartment * pCompartmentCol = (*itReactionCol)->getLargestCompartment();
575 C_FLOAT64 ScaleCol = (pCompartmentCol == NULL) ?
576 fabs((*itReactionCol)->getFlux()) :
577 fabs((*itReactionCol)->getFlux() / pCompartmentCol->
getValue());
579 if (fabs(ScaleCol) <= res)
581 *pScaled = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
585 *pScaled = (((*itReaction)->getFlux() < 0.0) ?
586 - std::numeric_limits<C_FLOAT64>::infinity() :
587 std::numeric_limits<C_FLOAT64>::infinity());
615 for (; pScaled != pScaledEnd; pScaledRowEnd +=
mScaledConcCC.
numCols(), ++pSum, ++pMax, ++pError)
617 for (; pScaled != pScaledRowEnd; ++pScaled)
619 success &= !isnan(*pScaled);
622 *pMax =
std::max(*pMax, fabs(*pScaled));
626 success &= *pError < resolution;
641 pError = Error.array();
643 for (; pScaled != pScaledEnd; pScaledRowEnd +=
mScaledConcCC.
numCols(), ++pSum, ++pMax, ++pError)
645 for (; pScaled != pScaledRowEnd; ++pScaled)
647 success &= !isnan(*pScaled);
650 *pMax =
std::max(*pMax, fabs(*pScaled));
654 success &= *pError < resolution;
676 bool SummationTheoremsOK =
false;
704 if (!SummationTheoremsOK)
749 if ((Fail = configBuffer.
getVariable(
"SSMCAUnscaled",
"C_INT16",
835 for (; it != end; ++it)
int dgetrf_(integer *m, integer *n, doublereal *a, integer *lda, integer *ipiv, integer *info)
CArrayAnnotation * mScaledElasticitiesAnn
bool build(const CMatrix< C_FLOAT64 > &matrix, size_t maxRank=C_INVALID_INDEX)
const CSteadyStateMethod::ReturnCode & getResult() const
CMatrix< C_FLOAT64 > mUnscaledConcCC
bool CalculateMCA(C_FLOAT64 res)
CSteadyStateMethod::ReturnCode mSSStatus
void setSteadyStateTask(CSteadyStateTask *pSteadyStateTask)
CMatrix< C_FLOAT64 > mScaledElasticities
int dlaset_(char *uplo, integer *m, integer *n, doublereal *alpha, doublereal *beta, doublereal *a, integer *lda)
virtual size_t size() const
CSteadyStateTask * mpSteadyStateTask
CArrayAnnotation * mScaledConcCCAnn
const CVector< C_FLOAT64 > & getParticleFlux() const
void updateSimulatedValues(const bool &updateMoieties)
virtual size_t numRows() const
CMatrix< C_FLOAT64 > mUnscaledFluxCC
bool leftMultiply(const CMatrix< C_FLOAT64 > &M, CMatrix< C_FLOAT64 > &P) const
void resize(size_t size, const bool ©=false)
void setDescription(const std::string &s)
CArrayAnnotation * mScaledFluxCCAnn
virtual void resizeAllMatrices()
const size_t & getNumIndependent() const
virtual bool elevateChildren()
const CMatrix< C_FLOAT64 > & getRedStoi() const
bool scaleMCA(const bool &status, C_FLOAT64 res)
bool rightMultiply(const C_FLOAT64 &alpha, const CMatrix< C_FLOAT64 > &M, CMatrix< C_FLOAT64 > &P) const
size_t getTotSteps() const
void setSteadyStateResolution(C_FLOAT64 factor)
const CModel * getModel() const
virtual bool isValidProblem(const CCopasiProblem *pProblem)
virtual bool isValidProblem(const CCopasiProblem *pProblem)
bool checkSummationTheorems(const C_FLOAT64 &resolution)
bool removeParameter(const std::string &name)
std::vector< CType * >::const_iterator const_iterator
CMatrix< C_FLOAT64 > mScaledConcCC
bool doRowPivot(CMatrix< C_FLOAT64 > &matrix) const
size_t getNumDependentReactionMetabs() const
bool createLinkMatrix(const bool &useSmallbone=false)
const Value & getValue() const
void setFactor(C_FLOAT64 factor)
bool setValue(const std::string &name, const CType &value)
void setCopasiVector(size_t d, const CCopasiContainer *v)
CModelEntity ** beginIndependent()
CCopasiParameter * getParameter(const std::string &name)
CMatrix< C_FLOAT64 > mReducedStoichiometry
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)
CMatrix< bool > mElasticityDependencies
virtual void resize(size_t rows, size_t cols, const bool ©=false)
size_t getNumIndependentReactionMetabs() const
bool doColumnPivot(CMatrix< C_FLOAT64 > &matrix) const
CArrayAnnotation * mUnscaledElasticitiesAnn
void calculateUnscaledElasticities(C_FLOAT64 res)
void setMode(size_t d, Mode m)
CCopasiVectorNS< CCompartment > & getCompartments()
void setDimensionDescription(size_t d, const std::string &s)
CMatrix< C_FLOAT64 > mUnscaledElasticities
const C_FLOAT64 & getValue() const
void setModel(CModel *model)
CArrayAnnotation * mUnscaledConcCCAnn
const CStateTemplate & getStateTemplate() const
CMCAMethod(const CCopasiContainer *pParent=NULL)
C_INT32 load(CReadConfig &configBuffer)
virtual size_t size() const
CCopasiVectorNS< CReaction > & getReactions()
bool undoColumnPivot(CMatrix< C_FLOAT64 > &matrix) const
const CLinkMatrix & getL0() const
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
static CMCAMethod * createMethod(CCopasiMethod::SubType subType=CCopasiMethod::mcaMethodReder)
C_INT32 getVariable(const std::string &name, const std::string &type, void *pout, CReadConfig::Mode mode=CReadConfig::NEXT)
const CCopasiVector< CMetab > & getMetabolitesX() const
const CMatrix< C_FLOAT64 > & getStoi() const
int dgetri_(integer *n, doublereal *a, integer *lda, integer *ipiv, doublereal *work, integer *lwork, integer *info)
virtual size_t numCols() const
CArrayAnnotation * mUnscaledFluxCCAnn
CMatrix< C_FLOAT64 > mScaledFluxCC
bool undoRowPivot(CMatrix< C_FLOAT64 > &matrix) const
void initializeParameter()
std::set< const CCopasiObject * > DataObjectSet
C_FLOAT64 mSteadyStateResolution
const CMatrix< C_FLOAT64 > & getJacobian() const
CModel * getModel() const
CModelEntity ** endIndependent()
bool calculateUnscaledConcentrationCC()
bool calculateUnscaledFluxCC(const bool &status)