#include <CMCAMethod.h>

Inheritance diagram for CMCAMethod:

Collaboration diagram for CMCAMethod:

Public Member Functions
bool	CalculateMCA (C_FLOAT64 res)

bool	calculateUnscaledConcentrationCC ()

void	calculateUnscaledElasticities (C_FLOAT64 res)

bool	calculateUnscaledFluxCC (const bool &status)

bool	checkSummationTheorems (const C_FLOAT64 &resolution)

	CMCAMethod (const CCopasiContainer *pParent=NULL)

	CMCAMethod (const CMCAMethod &src, const CCopasiContainer *pParent=NULL)

virtual bool	elevateChildren ()

const CModel *	getModel () const

const CMatrix< C_FLOAT64 > &	getScaledConcentrationCC () const

const CArrayAnnotation *	getScaledConcentrationCCAnn () const

const CMatrix< C_FLOAT64 > &	getScaledElasticities () const

const CArrayAnnotation *	getScaledElasticitiesAnn () const

const CMatrix< C_FLOAT64 > &	getScaledFluxCC () const

const CArrayAnnotation *	getScaledFluxCCAnn () const

const CSteadyStateMethod::ReturnCode &	getSteadyStateStatus () const

const CMatrix< C_FLOAT64 > &	getUnscaledConcentrationCC () const

const CArrayAnnotation *	getUnscaledConcentrationCCAnn () const

const CMatrix< C_FLOAT64 > &	getUnscaledElasticities () const

const CArrayAnnotation *	getUnscaledElasticitiesAnn () const

const CMatrix< C_FLOAT64 > &	getUnscaledFluxCC () const

const CArrayAnnotation *	getUnscaledFluxCCAnn () const

virtual bool	isValidProblem (const CCopasiProblem *pProblem)

C_INT32	load (CReadConfig &configBuffer)

virtual bool	process ()

virtual void	resizeAllMatrices ()

bool	scaleMCA (const bool &status, C_FLOAT64 res)

void	setFactor (C_FLOAT64 factor)

void	setModel (CModel *model)

void	setSteadyStateResolution (C_FLOAT64 factor)

void	setSteadyStateTask (CSteadyStateTask *pSteadyStateTask)

virtual	~CMCAMethod ()

Public Member Functions inherited from CCopasiMethod
	CCopasiMethod (const CCopasiMethod &src, const CCopasiContainer *pParent=NULL)

const CCopasiMethod::SubType &	getSubType () const

const CCopasiTask::Type &	getType () const

virtual void	load (CReadConfig &configBuffer, CReadConfig::Mode mode=CReadConfig::SEARCH)

virtual void	print (std::ostream *ostream) const

virtual void	printResult (std::ostream *ostream) const

virtual bool	setCallBack (CProcessReport *pCallBack)

virtual	~CCopasiMethod ()

Public Member Functions inherited from CCopasiParameterGroup
bool	addGroup (const std::string &name)

bool	addParameter (const CCopasiParameter &parameter)

bool	addParameter (const std::string &name, const CCopasiParameter::Type type)

template<class CType >
bool	addParameter (const std::string &name, const CCopasiParameter::Type type, const CType &value)

void	addParameter (CCopasiParameter *pParameter)

CCopasiParameterGroup *	assertGroup (const std::string &name)

template<class CType >
CCopasiParameter *	assertParameter (const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)

index_iterator	beginIndex () const

name_iterator	beginName () const

	CCopasiParameterGroup (const CCopasiParameterGroup &src, const CCopasiContainer *pParent=NULL)

	CCopasiParameterGroup (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &objectType="ParameterGroup")

void	clear ()

index_iterator	endIndex () const

name_iterator	endName () const

CCopasiParameterGroup *	getGroup (const std::string &name)

const CCopasiParameterGroup *	getGroup (const std::string &name) const

CCopasiParameterGroup *	getGroup (const size_t &index)

const CCopasiParameterGroup *	getGroup (const size_t &index) const

size_t	getIndex (const std::string &name) const

std::string	getKey (const std::string &name) const

std::string	getKey (const size_t &index) const

virtual const std::string &	getName (const size_t &index) const

virtual const CObjectInterface *	getObject (const CCopasiObjectName &cn) const

CCopasiParameter *	getParameter (const std::string &name)

const CCopasiParameter *	getParameter (const std::string &name) const

CCopasiParameter *	getParameter (const size_t &index)

const CCopasiParameter *	getParameter (const size_t &index) const

CCopasiParameter::Type	getType (const std::string &name) const

CCopasiParameter::Type	getType (const size_t &index) const

std::string	getUniqueParameterName (const CCopasiParameter *pParameter) const

const CCopasiParameter::Value &	getValue (const std::string &name) const

const CCopasiParameter::Value &	getValue (const size_t &index) const

CCopasiParameter::Value &	getValue (const std::string &name)

CCopasiParameter::Value &	getValue (const size_t &index)

CCopasiParameterGroup &	operator= (const CCopasiParameterGroup &rhs)

bool	removeParameter (const std::string &name)

bool	removeParameter (const size_t &index)

template<class CType >
bool	setValue (const std::string &name, const CType &value)

template<class CType >
bool	setValue (const size_t &index, const CType &value)

size_t	size () const

bool	swap (const size_t &iFrom, const size_t &iTo)

bool	swap (index_iterator &from, index_iterator &to)

virtual	~CCopasiParameterGroup ()

Public Member Functions inherited from CCopasiParameter
	CCopasiParameter (const CCopasiParameter &src, const CCopasiContainer *pParent=NULL)

	CCopasiParameter (const std::string &name, const Type &type, const void pValue=NULL, const CCopasiContainer pParent=NULL, const std::string &objectType="Parameter")

virtual CCopasiObjectName	getCN () const

virtual const std::string &	getKey () const

virtual std::string	getObjectDisplayName (bool regular=true, bool richtext=false) const

const CCopasiParameter::Type &	getType () const

const Value &	getValue () const

Value &	getValue ()

virtual void *	getValuePointer () const

CCopasiObject *	getValueReference () const

bool	isValidValue (const C_FLOAT64 &value) const

bool	isValidValue (const C_INT32 &value) const

bool	isValidValue (const unsigned C_INT32 &value) const

bool	isValidValue (const bool &value) const

bool	isValidValue (const std::string &value) const

bool	isValidValue (const CCopasiObjectName &value) const

bool	isValidValue (const std::vector< CCopasiParameter * > &value) const

CCopasiParameter &	operator= (const CCopasiParameter &rhs)

template<class CType >
bool	setValue (const CType &value)

bool	setValue (const std::vector< CCopasiParameter * > &value)

virtual	~CCopasiParameter ()

Public Member Functions inherited from CCopasiContainer
virtual bool	add (CCopasiObject *pObject, const bool &adopt=true)

	CCopasiContainer (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &type="CN", const unsigned C_INT32 &flag=CCopasiObject::Container)

	CCopasiContainer (const CCopasiContainer &src, const CCopasiContainer *pParent=NULL)

virtual std::string	getChildObjectUnits (const CCopasiObject *pObject) const

virtual const objectMap &	getObjects () const

virtual std::string	getUnits () const

virtual const CCopasiObject *	getValueObject () const

virtual bool	remove (CCopasiObject *pObject)

virtual	~CCopasiContainer ()

Public Member Functions inherited from CCopasiObject
void	addDirectDependency (const CCopasiObject *pObject)

	CCopasiObject (const CCopasiObject &src, const CCopasiContainer *pParent=NULL)

void	clearDirectDependencies ()

void	clearRefresh ()

bool	dependsOn (DataObjectSet candidates, const DataObjectSet &context=DataObjectSet()) const

void	getAllDependencies (DataObjectSet &dependencies, const DataObjectSet &context) const

virtual const DataObjectSet &	getDirectDependencies (const DataObjectSet &context=DataObjectSet()) const

CCopasiContainer *	getObjectAncestor (const std::string &type) const

CCopasiDataModel *	getObjectDataModel ()

const CCopasiDataModel *	getObjectDataModel () const

const std::string &	getObjectName () const

CCopasiContainer *	getObjectParent () const

const std::string &	getObjectType () const

virtual const CObjectInterface::ObjectSet &	getPrerequisites () const

virtual Refresh *	getRefresh () const

UpdateMethod *	getUpdateMethod () const

bool	hasCircularDependencies (DataObjectSet &candidates, DataObjectSet &verified, const DataObjectSet &context) const

bool	hasUpdateMethod () const

bool	isArray () const

bool	isContainer () const

bool	isDataModel () const

bool	isMatrix () const

bool	isNameVector () const

bool	isNonUniqueName () const

virtual bool	isPrerequisiteForContext (const CObjectInterface *pObject, const CMath::SimulationContextFlag &context, const CObjectInterface::ObjectSet &changedObjects) const

bool	isReference () const

bool	isRoot () const

bool	isSeparator () const

bool	isStaticString () const

bool	isValueBool () const

bool	isValueDbl () const

bool	isValueInt () const

bool	isValueInt64 () const

bool	isValueString () const

bool	isVector () const

virtual bool	mustBeDeleted (const DataObjectSet &deletedObjects) const

void	removeDirectDependency (const CCopasiObject *pObject)

void	setDirectDependencies (const DataObjectSet &directDependencies)

bool	setObjectName (const std::string &name)

virtual bool	setObjectParent (const CCopasiContainer *pParent)

void	setObjectValue (const C_FLOAT64 &value)

void	setObjectValue (const C_INT32 &value)

void	setObjectValue (const bool &value)

template<class CType >
void	setRefresh (CType pType, void(CType::method)(void))

template<class CType >
void	setUpdateMethod (CType pType, void(CType::method)(const C_FLOAT64 &))

template<class CType >
void	setUpdateMethod (CType pType, void(CType::method)(const C_INT32 &))

template<class CType >
void	setUpdateMethod (CType pType, void(CType::method)(const bool &))

virtual	~CCopasiObject ()

Public Member Functions inherited from CObjectInterface
	CObjectInterface ()

virtual	~CObjectInterface ()

Static Public Member Functions
static CMCAMethod *	createMethod (CCopasiMethod::SubType subType=CCopasiMethod::mcaMethodReder)

Static Public Member Functions inherited from CCopasiObject
static std::vector< Refresh * >	buildUpdateSequence (const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())

static void	setRenameHandler (CRenameHandler *rh)

Private Member Functions
bool	createLinkMatrix (const bool &useSmallbone=false)

void	initializeParameter ()

void	initObjects ()

Private Attributes
CMatrix< bool >	mElasticityDependencies

C_FLOAT64	mFactor

CLinkMatrix	mLinkZero

CModel *	mpModel

CSteadyStateTask *	mpSteadyStateTask

bool *	mpUseReeder

bool *	mpUseSmallbone

CMatrix< C_FLOAT64 >	mReducedStoichiometry

CMatrix< C_FLOAT64 >	mScaledConcCC

CArrayAnnotation *	mScaledConcCCAnn

CMatrix< C_FLOAT64 >	mScaledElasticities

CArrayAnnotation *	mScaledElasticitiesAnn

CMatrix< C_FLOAT64 >	mScaledFluxCC

CArrayAnnotation *	mScaledFluxCCAnn

C_INT16	mSSReder

CSteadyStateMethod::ReturnCode	mSSStatus

std::vector< C_FLOAT64 >	mSsx

C_FLOAT64	mSteadyStateResolution

CMatrix< C_FLOAT64 >	mUnscaledConcCC

CArrayAnnotation *	mUnscaledConcCCAnn

CMatrix< C_FLOAT64 >	mUnscaledElasticities

CArrayAnnotation *	mUnscaledElasticitiesAnn

CMatrix< C_FLOAT64 >	mUnscaledFluxCC

CArrayAnnotation *	mUnscaledFluxCCAnn

Additional Inherited Members
Public Types inherited from CCopasiMethod
enum	SubType { unset = 0, RandomSearch, RandomSearchMaster, SimulatedAnnealing, CoranaWalk, DifferentialEvolution, ScatterSearch, GeneticAlgorithm, EvolutionaryProgram, SteepestDescent, HybridGASA, GeneticAlgorithmSR, HookeJeeves, LevenbergMarquardt, NelderMead, SRES, Statistics, ParticleSwarm, Praxis, TruncatedNewton, Newton, deterministic, LSODAR, directMethod, stochastic, tauLeap, adaptiveSA, hybrid, hybridLSODA, hybridODE45, DsaLsodar, tssILDM, tssILDMModified, tssCSP, mcaMethodReder, scanMethod, lyapWolf, sensMethod, EFMAlgorithm, EFMBitPatternTreeAlgorithm, EFMBitPatternAlgorithm, Householder, crossSectionMethod, linearNoiseApproximation }

Public Types inherited from CCopasiParameterGroup
typedef parameterGroup::iterator	index_iterator

typedef CCopasiContainer::objectMap::iterator	name_iterator

typedef std::vector < CCopasiParameter * >	parameterGroup

Public Types inherited from CCopasiParameter
enum	Type { DOUBLE = 0, UDOUBLE, INT, UINT, BOOL, GROUP, STRING, CN, KEY, FILE, EXPRESSION, INVALID }

Public Types inherited from CCopasiContainer
typedef std::multimap < std::string, CCopasiObject * >	objectMap

Public Types inherited from CCopasiObject
typedef std::set< const CCopasiObject * >	DataObjectSet

typedef std::vector< Refresh * >	DataUpdateSequence

Public Types inherited from CObjectInterface
typedef std::set< const CObjectInterface * >	ObjectSet

typedef std::vector < CObjectInterface * >	UpdateSequence

Static Public Attributes inherited from CCopasiMethod
static const std::string	SubTypeName []

static const char *	XMLSubType []

Static Public Attributes inherited from CCopasiParameter
static const std::string	TypeName []

static const char *	XMLType []

Static Public Attributes inherited from CCopasiContainer
static const std::vector < CCopasiContainer * >	EmptyList

Protected Types inherited from CCopasiObject
enum	Flag { Container = 0x1, Vector = 0x2, Matrix = 0x4, NameVector = 0x8, Reference = 0x10, ValueBool = 0x20, ValueInt = 0x40, ValueInt64 = 0x80, ValueDbl = 0x100, NonUniqueName = 0x200, StaticString = 0x400, ValueString = 0x800, Separator = 0x1000, ModelEntity = 0x2000, Array = 0x4000, DataModel = 0x8000, Root = 0x10000, Gui = 0x20000 }

Protected Member Functions inherited from CCopasiMethod
	CCopasiMethod (const CCopasiTask::Type &taskType, const SubType &subType, const CCopasiContainer *pParent=NULL)

Protected Member Functions inherited from CCopasiParameterGroup
	CCopasiParameterGroup ()

Protected Member Functions inherited from CCopasiContainer
template<class CType >
CCopasiObject *	addMatrixReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)

template<class CType >
CCopasiObject *	addObjectReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)

template<class CType >
CCopasiObject *	addVectorReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)

void	initObjects ()

Protected Member Functions inherited from CCopasiObject
	CCopasiObject ()

	CCopasiObject (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &type="CN", const unsigned C_INT32 &flag=0)

Protected Attributes inherited from CCopasiMethod
CProcessReport *	mpCallBack

Protected Attributes inherited from CCopasiParameter
std::string	mKey

CCopasiObject *	mpValueReference

size_t	mSize

Value	mValue

Protected Attributes inherited from CCopasiContainer
objectMap	mObjects

Static Protected Attributes inherited from CCopasiObject
static CRenameHandler *	smpRenameHandler = NULL

Detailed Description

Definition at line 27 of file CMCAMethod.h.

Constructor & Destructor Documentation

CMCAMethod::CMCAMethod ( const CCopasiContainer * pParent = NULL )

Default constructor

Parameters

const CCopasiContainer * pParent (Default: NULL)

Default constructor

Definition at line 42 of file CMCAMethod.cpp.

References initializeParameter(), and initObjects().

Referenced by createMethod().

                                                      :
   CCopasiMethod(CCopasiTask::mca, CCopasiMethod::mcaMethodReder, pParent),
   mpModel(NULL),
   mpUseReeder(NULL),
   mpUseSmallbone(NULL),
   mFactor(1.0e-9),
   mSteadyStateResolution(1.0e-9),
   mSSStatus(CSteadyStateMethod::notFound),
   mpSteadyStateTask(NULL),
   mLinkZero(),
   mReducedStoichiometry(),
   mElasticityDependencies()
 {
   initializeParameter();
   initObjects();
 }

CMCAMethod::CMCAMethod	(	const CMCAMethod &	src,
		const CCopasiContainer *	pParent = `NULL`
	)

Copy constructor

Parameters

const	CMCAMethod & src
const	CCopasiContainer * pParent (Default: NULL)

Definition at line 59 of file CMCAMethod.cpp.

References initializeParameter(), and initObjects().

                                                         :
   CCopasiMethod(src, pParent),
   mpModel(NULL),
   mpUseReeder(NULL),
   mpUseSmallbone(NULL),
   mFactor(src.mFactor),
   mSteadyStateResolution(src.mSteadyStateResolution),
   mSSStatus(CSteadyStateMethod::notFound),
   mpSteadyStateTask(NULL),
   mLinkZero(src.mLinkZero),
   mReducedStoichiometry(src.mReducedStoichiometry),
   mElasticityDependencies(src.mElasticityDependencies)
 {
   initializeParameter();
   initObjects();
 }

CMCAMethod::~CMCAMethod ( )

virtual

Destructor

Deconstructor

Definition at line 132 of file CMCAMethod.cpp.

References DESTRUCTOR_TRACE.

 {
   DESTRUCTOR_TRACE;
   //delSsipvt();
 }

Member Function Documentation

bool CMCAMethod::CalculateMCA ( C_FLOAT64 res )

the steady state MCA entry point

Parameters

ss_solution	refer to steady-state solution
refer	to the resolution

Definition at line 673 of file CMCAMethod.cpp.

References calculateUnscaledConcentrationCC(), calculateUnscaledElasticities(), calculateUnscaledFluxCC(), checkSummationTheorems(), createLinkMatrix(), CSteadyStateMethod::found, MCMCA, mpModel, mpUseReeder, mpUseSmallbone, mSSStatus, mUnscaledConcCC, mUnscaledFluxCC, scaleMCA(), and CCopasiMessage::WARNING.

Referenced by process().

 {
   bool success = true;
   bool SummationTheoremsOK = false;
 
   assert(mpModel);
 
   calculateUnscaledElasticities(res);
 
   if (mSSStatus == CSteadyStateMethod::found)
     {
       if (*mpUseReeder)
         {
           createLinkMatrix();
           success &= calculateUnscaledConcentrationCC();
           success &= calculateUnscaledFluxCC(success);
           success &= scaleMCA(success, res);
           SummationTheoremsOK = checkSummationTheorems(res);
         }
 
       if (*mpUseSmallbone && !SummationTheoremsOK)
         {
           success = true;
 
           createLinkMatrix(true);
           success &= calculateUnscaledConcentrationCC();
           success &= calculateUnscaledFluxCC(success);
           success &= scaleMCA(success, res);
           SummationTheoremsOK = checkSummationTheorems(res);
         }
 
       if (!SummationTheoremsOK)
         {
           CCopasiMessage(CCopasiMessage::WARNING, MCMCA + 1);
         }
     }
   else
     {
       mUnscaledConcCC = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
       mUnscaledFluxCC = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
     }
 
   return success;
 }

bool CMCAMethod::calculateUnscaledConcentrationCC ( )

Definition at line 311 of file CMCAMethod.cpp.

References CVectorCore< CType >::array(), CMatrix< CType >::array(), C_FLOAT64, C_INT, dgemm_(), dgetrf_(), dgetri_(), CLinkMatrix::doColumnPivot(), K, CLinkMatrix::leftMultiply(), mLinkZero, mpModel, mReducedStoichiometry, mUnscaledConcCC, mUnscaledElasticities, CMatrix< CType >::numCols(), CMatrix< CType >::numRows(), CMatrix< CType >::resize(), CVector< CType >::resize(), CLinkMatrix::rightMultiply(), CLinkMatrix::undoColumnPivot(), and CLinkMatrix::undoRowPivot().

Referenced by CalculateMCA().

 {
   assert(mpModel);
 
   // TODO CRITICAL We must not use the reduced stoichiometry matrix
   // Instead use N * mLinkZero
 
   // Calculate RedStoi * mUnscaledElasticities;
   // Note the columns of mUnscaledElasticities must be reordered
 
   mLinkZero.doColumnPivot(mUnscaledElasticities);
 
   // Initialize the unscaled concentration control coefficients to 0.0
   mUnscaledConcCC = 0.0;
 
   // aux1 := mUnscaledElasticities * L
   CMatrix<C_FLOAT64> aux1;
   mLinkZero.rightMultiply(1.0, mUnscaledElasticities, aux1);
 
   // We can now undo the column pivoting
   mLinkZero.undoColumnPivot(mUnscaledElasticities);
 
   assert(mReducedStoichiometry.numCols() == aux1.numRows());
 
   // aux2 := RedStoi * aux1
   // DGEMM (TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
   // C := alpha A B + beta C
   CMatrix<C_FLOAT64> aux2(mReducedStoichiometry.numRows(), aux1.numCols());
 
   char TRANSA = 'N';
   char TRANSB = 'N';
   C_INT M = (C_INT) aux2.numCols(); /* LDA, LDC */
   C_INT N = (C_INT) aux2.numRows();
   C_INT K = (C_INT) mReducedStoichiometry.numCols();
   C_FLOAT64 Alpha = 1.0;
   C_INT LDA = (C_INT) aux1.numCols();
   C_INT LDB = (C_INT) mReducedStoichiometry.numCols();
   C_FLOAT64 Beta = 0.0;
   C_INT LDC = (C_INT) aux2.numCols();
 
   dgemm_(&TRANSA, &TRANSB, &M, &N, &K, &Alpha, aux1.array(), &LDA,
          mReducedStoichiometry.array(), &LDB, &Beta, aux2.array(), &LDC);
 
   // Invert aux2
   C_INT info;
   CVector<C_INT> Ipiv(M);
 
   // LU decomposition of aux2 (for inversion)
   dgetrf_(&M, &M, aux2.array(), &M, Ipiv.array(), &info);
 
   if (info != 0) return false;
 
   C_INT lwork = -1; // Instruct dgetri_ to determine work array size.
   CVector< C_FLOAT64 > work(1);
 
   dgetri_(&M, aux2.array(), &M, Ipiv.array(), work.array(), &lwork, &info);
 
   lwork = (C_INT) work[0];
   work.resize(lwork);
 
   // now invert aux2 (result in aux2)
   dgetri_(&M, aux2.array(), &M, Ipiv.array(), work.array(), &lwork, &info);
 
   if (info != 0) return false;
 
   // aux1 := -1.0 * aux2 * RedStoi
 
   aux1.resize(aux2.numRows(), mReducedStoichiometry.numCols());
 
   M = (C_INT) aux1.numCols();
   N = (C_INT) aux1.numRows();
   K = (C_INT) aux2.numCols();
   Alpha = -1.0;
   LDA = (C_INT) mReducedStoichiometry.numCols();
   LDB = (C_INT) aux2.numCols();
   Beta = 0.0;
   LDC = (C_INT) aux1.numCols();
 
   // DGEMM (TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
   // C := alpha A B + beta C
   dgemm_(&TRANSA, &TRANSB, &M, &N, &K, &Alpha, mReducedStoichiometry.array(), &LDA,
          aux2.array(), &LDB, &Beta, aux1.array(), &LDC);
 
   // mUnscaledConcCC := L * aux1
   mLinkZero.leftMultiply(aux1, mUnscaledConcCC);
 
   // We need to swap the rows since they are with respect to reordered stoichiometry .
   mLinkZero.undoRowPivot(mUnscaledConcCC);
 
   return true;
 }

void CMCAMethod::calculateUnscaledElasticities ( C_FLOAT64 res )

Definition at line 203 of file CMCAMethod.cpp.

References CVectorCore< CType >::array(), CMatrix< CType >::array(), CCopasiVector< T >::begin(), C_FLOAT64, CCopasiVector< T >::end(), CModel::getMetabolitesX(), CModel::getNumDependentReactionMetabs(), CModel::getNumIndependentReactionMetabs(), CModel::getParticleFlux(), CModel::getReactions(), mElasticityDependencies, mFactor, min, mpModel, mUnscaledElasticities, CMatrix< CType >::size(), CCopasiVector< T >::size(), and CModel::updateSimulatedValues().

Referenced by CalculateMCA().

 {
   assert(mpModel);
 
   // We need the number of metabolites determined by reactions.
   size_t numMetabs =
     mpModel->getNumIndependentReactionMetabs() + mpModel->getNumDependentReactionMetabs();
 
   CCopasiVector< CMetab >::iterator itSpecies = mpModel->getMetabolitesX().begin();
   CCopasiVector< CMetab >::iterator endSpecies = itSpecies + numMetabs;
 
   // Calculate the dependencies of the elasticities this is helpful for scaling to determine
   // whether 0/0 is 0 or NaN
 
   bool * pElasticityDependency = mElasticityDependencies.array();
   CCopasiVector< CReaction >::const_iterator itReac = mpModel->getReactions().begin();
   CCopasiVector< CReaction >::const_iterator endReac = mpModel->getReactions().end();
 
   for (; itReac != endReac; ++itReac)
     {
       CCopasiObject::DataObjectSet ReactionDependencies;
       (*itReac)->getParticleFluxReference()->getAllDependencies(ReactionDependencies, DataObjectSet());
 
       itSpecies = mpModel->getMetabolitesX().begin();
 
       for (; itSpecies != endSpecies; ++itSpecies, ++pElasticityDependency)
         {
           *pElasticityDependency = (ReactionDependencies.find((*itSpecies)->getValueReference()) != ReactionDependencies.end());
         }
     }
 
   const CVector< C_FLOAT64 > & ParticleFlux = mpModel->getParticleFlux();
 
   //   mUnscaledElasticities.resize(numReacs, numMetabs);
   C_FLOAT64 * pElasticity;
 
   C_FLOAT64 * pElasticityEnd = mUnscaledElasticities.array() + mUnscaledElasticities.size();
 
   C_FLOAT64 Store, InvDelta;
 
   C_FLOAT64 X1, X2;
 
   // Arrays to store function value
   size_t numReacs = mpModel->getReactions().size();
 
   CVector< C_FLOAT64 > Y1(numReacs);
 
   C_FLOAT64 * pY1;
 
   CVector< C_FLOAT64 > Y2(numReacs);
 
   C_FLOAT64 * pY2;
 
   // calculate elasticities
   itSpecies = mpModel->getMetabolitesX().begin();
 
   for (size_t j = 0; itSpecies != endSpecies; ++itSpecies, ++j)
     {
       Store = (*itSpecies)->getValue();
 
       // We only need to make sure that we do not have an underflow problem
       if (fabs(Store) < 100 * std::numeric_limits< C_FLOAT64 >::min())
         {
           X1 = 0.0;
 
           if (Store < 0.0)
             X2 = -200.0 * std::numeric_limits< C_FLOAT64 >::min();
           else
             X2 = 200.0 * std::numeric_limits< C_FLOAT64 >::min();;
         }
       else
         {
           X1 = Store * (1.0 + mFactor);
           X2 = Store * (1.0 - mFactor);
         }
 
       InvDelta = 1.0 / (X1 - X2);
 
       // let's take X+dx
       (*itSpecies)->setValue(X1);
       mpModel->updateSimulatedValues(false); // TODO test if true or false should be used.
 
       // get the fluxes
       Y1 = ParticleFlux;
 
       // now X-dx
       (*itSpecies)->setValue(X2);
       mpModel->updateSimulatedValues(false); // TODO test if true or false should be used.
 
       // get the fluxes
       Y2 = ParticleFlux;
 
       // set column j of Dxv
       pElasticity = mUnscaledElasticities.array() + j;
       pY1 = Y1.array();
       pY2 = Y2.array();
 
       for (; pElasticity < pElasticityEnd; pElasticity += numMetabs, ++pY1, ++pY2)
         * pElasticity = (*pY1 - *pY2) * InvDelta;
 
       // restore the value of the species
       (*itSpecies)->setValue(Store);
     }
 
   // make sure the fluxes are correct afterwards (needed for scaling of the MCA results)
   mpModel->updateSimulatedValues(false);
 }

bool CMCAMethod::calculateUnscaledFluxCC ( const bool & status )

Definition at line 403 of file CMCAMethod.cpp.

References CMatrix< CType >::array(), C_FLOAT64, C_INT, dgemm_(), dlaset_(), K, mpModel, mUnscaledConcCC, mUnscaledElasticities, mUnscaledFluxCC, CMatrix< CType >::numCols(), and CMatrix< CType >::numRows().

Referenced by CalculateMCA().

 {
   assert(mpModel);
   //size_t i, j, k;
 
   // mUnscaledFluxCC := I + mUnscaledElasticities * mUnscaledConcCC
 
   char UPLO = 'A';
   C_INT M = mUnscaledFluxCC.numRows();
   C_FLOAT64 Alpha = 0.0;
   C_FLOAT64 Beta = 1.0;
 
   // Initialize mUnscaledFluxCC to the identity matrix;
   dlaset_(&UPLO, &M, &M, &Alpha, &Beta, mUnscaledFluxCC.array(), &M);
 
   if (status)
     {
       char TRANSA = 'N';
       char TRANSB = 'N';
       M = (C_INT) mUnscaledFluxCC.numCols(); /* LDA, LDC */
       C_INT N = (C_INT) mUnscaledFluxCC.numRows();
       C_INT K = (C_INT) mUnscaledElasticities.numCols();
 
       C_INT LDA = (C_INT) mUnscaledConcCC.numCols();
       C_INT LDB = (C_INT) mUnscaledElasticities.numCols();
       C_INT LDC = (C_INT) mUnscaledFluxCC.numCols();
 
       Alpha = 1.0;
       Beta = 1.0;
 
       dgemm_(&TRANSA, &TRANSB, &M, &N, &K, &Alpha, mUnscaledConcCC.array(), &LDA,
              mUnscaledElasticities.array(), &LDB, &Beta, mUnscaledFluxCC.array(), &LDC);
     }
 
   return status;
 }

bool CMCAMethod::checkSummationTheorems ( const C_FLOAT64 & resolution )

Check whether the summation theorems hold.

Definition at line 598 of file CMCAMethod.cpp.

References CVectorCore< CType >::array(), CMatrix< CType >::array(), C_FLOAT64, max, min, mScaledConcCC, mScaledFluxCC, CMatrix< CType >::numCols(), CMatrix< CType >::numRows(), CVector< CType >::resize(), and CMatrix< CType >::size().

Referenced by CalculateMCA().

 {
   bool success = true;
 
   C_FLOAT64 * pScaled = mScaledConcCC.array();
   C_FLOAT64 * pScaledRowEnd = pScaled + mScaledConcCC.numCols();
   C_FLOAT64 * pScaledEnd = pScaled + mScaledConcCC.size();
 
   CVector< C_FLOAT64 > Sum(mScaledConcCC.numRows());
   CVector< C_FLOAT64 > Max(mScaledConcCC.numRows());
   CVector< C_FLOAT64 > Error(mScaledConcCC.numRows());
   Sum = 0.0;
   Max = 0.0;
   C_FLOAT64 * pSum = Sum.array();
   C_FLOAT64 * pMax = Max.array();
   C_FLOAT64 * pError = Error.array();
 
   for (; pScaled != pScaledEnd; pScaledRowEnd += mScaledConcCC.numCols(), ++pSum, ++pMax, ++pError)
     {
       for (; pScaled != pScaledRowEnd; ++pScaled)
         {
           success &= !isnan(*pScaled);
 
           *pSum += *pScaled;
           *pMax = std::max(*pMax, fabs(*pScaled));
         }
 
       *pError = (*pMax > 100.0 * std::numeric_limits< C_FLOAT64 >::min()) ? fabs(*pSum) / *pMax : 0.0;
       success &= *pError < resolution;
     }
 
   pScaled = mScaledFluxCC.array();
   pScaledRowEnd = pScaled + mScaledFluxCC.numCols();
   pScaledEnd = pScaled + mScaledFluxCC.size();
 
   Sum.resize(mScaledFluxCC.numRows());
   Max.resize(mScaledFluxCC.numRows());
   Error.resize(mScaledFluxCC.numRows());
   Sum = 0.0;
   Max = 0.0;
 
   pSum = Sum.array();
   pMax = Max.array();
   pError = Error.array();
 
   for (; pScaled != pScaledEnd; pScaledRowEnd += mScaledConcCC.numCols(), ++pSum, ++pMax, ++pError)
     {
       for (; pScaled != pScaledRowEnd; ++pScaled)
         {
           success &= !isnan(*pScaled);
 
           *pSum += *pScaled;
           *pMax = std::max(*pMax, fabs(*pScaled));
         }
 
       *pError = (*pMax > 100.0 * std::numeric_limits< C_FLOAT64 >::min()) ? fabs(1.0 - *pSum) / *pMax : 0.0;
       success &= *pError < resolution;
     }
 
   return success;
 }

bool CMCAMethod::createLinkMatrix ( const bool & useSmallbone = false )

private

Definition at line 718 of file CMCAMethod.cpp.

References CLinkMatrix::build(), CLinkMatrix::doRowPivot(), CSteadyStateTask::getJacobian(), CModel::getL0(), CLinkMatrix::getNumIndependent(), CModel::getNumIndependentReactionMetabs(), CModel::getRedStoi(), CModel::getStoi(), mLinkZero, mpModel, mpSteadyStateTask, mReducedStoichiometry, CMatrix< CType >::numCols(), and CMatrix< CType >::resize().

Referenced by CalculateMCA().

 {
   if (mpModel == NULL ||
       mpSteadyStateTask == NULL)
     {
       return false;
     }
 
   if (useSmallbone)
     {
       mLinkZero.build(mpSteadyStateTask->getJacobian(), mpModel->getNumIndependentReactionMetabs());
       mReducedStoichiometry = mpModel->getStoi();
       mLinkZero.doRowPivot(mReducedStoichiometry);
       mReducedStoichiometry.resize(mLinkZero.getNumIndependent(), mReducedStoichiometry.numCols(), true);
     }
   else
     {
       mLinkZero = mpModel->getL0();
       mReducedStoichiometry = mpModel->getRedStoi();
     }
 
   return true;
 }

CMCAMethod * CMCAMethod::createMethod ( CCopasiMethod::SubType subType = CCopasiMethod::mcaMethodReder )

static

Create a MCA method. Note: the returned object has to be released after use with delete

Definition at line 34 of file CMCAMethod.cpp.

References CMCAMethod().

Referenced by CMCATask::createMethod().

 {
   return new CMCAMethod();
 }

bool CMCAMethod::elevateChildren ( )

virtual

This methods must be called to elevate subgroups to derived objects. The default implementation does nothing.

Returns: bool success

Reimplemented from CCopasiParameterGroup.

Definition at line 154 of file CMCAMethod.cpp.

References initializeParameter().

 {
   initializeParameter();
   return true;
 }

const CModel * CMCAMethod::getModel ( ) const

Get the Model

Definition at line 797 of file CMCAMethod.cpp.

References mpModel.

Referenced by CMCAResultSubwidget::loadAll().

 {
   return this->mpModel;
 }

const CMatrix<C_FLOAT64>& CMCAMethod::getScaledConcentrationCC ( ) const

inline

Definition at line 80 of file CMCAMethod.h.

References mScaledConcCC.

81 {return mScaledConcCC;}

CMCAMethod::mScaledConcCC

CMatrix< C_FLOAT64 > mScaledConcCC

Definition: CMCAMethod.h:193

const CArrayAnnotation* CMCAMethod::getScaledConcentrationCCAnn ( ) const

inline

Definition at line 92 of file CMCAMethod.h.

References mScaledConcCCAnn.

Referenced by CMCAResultSubwidget::loadConcentrationCCs(), and CMCATask::printResult().

93 {return mScaledConcCCAnn;}

CMCAMethod::mScaledConcCCAnn

CArrayAnnotation * mScaledConcCCAnn

Definition: CMCAMethod.h:194

const CMatrix<C_FLOAT64>& CMCAMethod::getScaledElasticities ( ) const

inline

Definition at line 77 of file CMCAMethod.h.

References mScaledElasticities.

78 {return mScaledElasticities;}

CMCAMethod::mScaledElasticities

CMatrix< C_FLOAT64 > mScaledElasticities

Definition: CMCAMethod.h:190

const CArrayAnnotation* CMCAMethod::getScaledElasticitiesAnn ( ) const

inline

Definition at line 88 of file CMCAMethod.h.

References mScaledElasticitiesAnn.

Referenced by CMCAResultSubwidget::loadElasticities(), and CMCATask::printResult().

89 {return mScaledElasticitiesAnn;}

CMCAMethod::mScaledElasticitiesAnn

CArrayAnnotation * mScaledElasticitiesAnn

Definition: CMCAMethod.h:191

const CMatrix<C_FLOAT64>& CMCAMethod::getScaledFluxCC ( ) const

inline

Definition at line 83 of file CMCAMethod.h.

References mScaledFluxCC.

84 {return mScaledFluxCC;}

CMCAMethod::mScaledFluxCC

CMatrix< C_FLOAT64 > mScaledFluxCC

Definition: CMCAMethod.h:196

const CArrayAnnotation* CMCAMethod::getScaledFluxCCAnn ( ) const

inline

Definition at line 96 of file CMCAMethod.h.

References mScaledFluxCCAnn.

Referenced by CMCAResultSubwidget::loadFluxCCs(), and CMCATask::printResult().

97 {return mScaledFluxCCAnn;}

CMCAMethod::mScaledFluxCCAnn

CArrayAnnotation * mScaledFluxCCAnn

Definition: CMCAMethod.h:197

const CSteadyStateMethod::ReturnCode& CMCAMethod::getSteadyStateStatus ( ) const

inline

Definition at line 147 of file CMCAMethod.h.

References mSSStatus.

Referenced by CMCAResultSubwidget::loadAll(), and CMCATask::printResult().

148 {return mSSStatus;}

CMCAMethod::mSSStatus

CSteadyStateMethod::ReturnCode mSSStatus

Definition: CMCAMethod.h:217

const CMatrix<C_FLOAT64>& CMCAMethod::getUnscaledConcentrationCC ( ) const

inline

Definition at line 67 of file CMCAMethod.h.

References mUnscaledConcCC.

68 {return mUnscaledConcCC;}

CMCAMethod::mUnscaledConcCC

CMatrix< C_FLOAT64 > mUnscaledConcCC

Definition: CMCAMethod.h:184

const CArrayAnnotation* CMCAMethod::getUnscaledConcentrationCCAnn ( ) const

inline

Definition at line 90 of file CMCAMethod.h.

References mUnscaledConcCCAnn.

Referenced by CMCAResultSubwidget::loadConcentrationCCs(), and CMCATask::printResult().

91 {return mUnscaledConcCCAnn;}

CMCAMethod::mUnscaledConcCCAnn

CArrayAnnotation * mUnscaledConcCCAnn

Definition: CMCAMethod.h:185

const CMatrix<C_FLOAT64>& CMCAMethod::getUnscaledElasticities ( ) const

inline

Definition at line 62 of file CMCAMethod.h.

References mUnscaledElasticities.

63 {return mUnscaledElasticities;}

CMCAMethod::mUnscaledElasticities

CMatrix< C_FLOAT64 > mUnscaledElasticities

Definition: CMCAMethod.h:181

const CArrayAnnotation* CMCAMethod::getUnscaledElasticitiesAnn ( ) const

inline

Definition at line 86 of file CMCAMethod.h.

References mUnscaledElasticitiesAnn.

Referenced by CMCAResultSubwidget::loadElasticities(), and CMCATask::printResult().

87 {return mUnscaledElasticitiesAnn;}

CMCAMethod::mUnscaledElasticitiesAnn

CArrayAnnotation * mUnscaledElasticitiesAnn

Definition: CMCAMethod.h:182

const CMatrix<C_FLOAT64>& CMCAMethod::getUnscaledFluxCC ( ) const

inline

Definition at line 72 of file CMCAMethod.h.

References mUnscaledFluxCC.

73 {return mUnscaledFluxCC;}

CMCAMethod::mUnscaledFluxCC

CMatrix< C_FLOAT64 > mUnscaledFluxCC

Definition: CMCAMethod.h:187

const CArrayAnnotation* CMCAMethod::getUnscaledFluxCCAnn ( ) const

inline

Definition at line 94 of file CMCAMethod.h.

References mUnscaledFluxCCAnn.

Referenced by CMCAResultSubwidget::loadFluxCCs(), and CMCATask::printResult().

95 {return mUnscaledFluxCCAnn;}

CMCAMethod::mUnscaledFluxCCAnn

CArrayAnnotation * mUnscaledFluxCCAnn

Definition: CMCAMethod.h:188

void CMCAMethod::initializeParameter ( )

private

Initialize the method parameter

Definition at line 138 of file CMCAMethod.cpp.

References CCopasiParameterGroup::assertParameter(), CCopasiParameter::BOOL, CCopasiParameterGroup::getParameter(), CCopasiParameter::getValue(), mpUseReeder, mpUseSmallbone, CCopasiParameter::Value::pBOOL, CCopasiParameter::Value::pUDOUBLE, CCopasiParameterGroup::removeParameter(), CCopasiParameterGroup::setValue(), and CCopasiParameter::UDOUBLE.

Referenced by CMCAMethod(), and elevateChildren().

 {
   CCopasiParameter *pParm;
 
   assertParameter("Modulation Factor", CCopasiParameter::UDOUBLE, 1.0e-009);
 
   if ((pParm = getParameter("MCA.ModulationFactor")) != NULL)
     {
       setValue("Modulation Factor", *pParm->getValue().pUDOUBLE);
       removeParameter("MCA.ModulationFactor");
     }
 
   mpUseReeder = assertParameter("Use Reeder", CCopasiParameter::BOOL, true)->getValue().pBOOL;
   mpUseSmallbone = assertParameter("Use Smallbone", CCopasiParameter::BOOL, true)->getValue().pBOOL;
 }

void CMCAMethod::initObjects ( )

private

Definition at line 77 of file CMCAMethod.cpp.

References mScaledConcCC, mScaledConcCCAnn, mScaledElasticities, mScaledElasticitiesAnn, mScaledFluxCC, mScaledFluxCCAnn, mUnscaledConcCC, mUnscaledConcCCAnn, mUnscaledElasticities, mUnscaledElasticitiesAnn, mUnscaledFluxCC, mUnscaledFluxCCAnn, CArrayAnnotation::setDescription(), CArrayAnnotation::setDimensionDescription(), CArrayAnnotation::setMode(), and CArrayAnnotation::VECTOR.

Referenced by CMCAMethod().

 {
   CArrayAnnotation *
   tmp = new CArrayAnnotation("Unscaled elasticities", this,
                              new CCopasiMatrixInterface<CMatrix<C_FLOAT64> >(&mUnscaledElasticities), true);
   tmp->setMode(CArrayAnnotation::VECTOR);
   tmp->setDescription("Unscaled elasticity matrix");
   tmp->setDimensionDescription(0, "Reactions (reduced system)");
   tmp->setDimensionDescription(1, "Species (reduced system)");
   mUnscaledElasticitiesAnn = tmp;
 
   tmp = new CArrayAnnotation("Unscaled concentration control coefficients", this,
                              new CCopasiMatrixInterface<CMatrix<C_FLOAT64> >(&mUnscaledConcCC), true);
   tmp->setMode(CArrayAnnotation::VECTOR);
   tmp->setDescription("Unscaled concentration control coefficients");
   tmp->setDimensionDescription(0, "Species (reduced system)");
   tmp->setDimensionDescription(1, "Reactions (reduced system)");
   mUnscaledConcCCAnn = tmp;
 
   tmp = new CArrayAnnotation("Unscaled flux control coefficients", this,
                              new CCopasiMatrixInterface<CMatrix<C_FLOAT64> >(&mUnscaledFluxCC), true);
   tmp->setMode(CArrayAnnotation::VECTOR);
   tmp->setDescription("Unscaled flux control coefficients");
   tmp->setDimensionDescription(0, "Reactions (reduced system)");
   tmp->setDimensionDescription(1, "Reactions (reduced system)");
   mUnscaledFluxCCAnn = tmp;
 
   tmp = new CArrayAnnotation("Scaled elasticities", this,
                              new CCopasiMatrixInterface<CMatrix<C_FLOAT64> >(&mScaledElasticities), true);
   tmp->setMode(CArrayAnnotation::VECTOR);
   tmp->setDescription("Scaled elasticity matrix");
   tmp->setDimensionDescription(0, "Reactions (reduced system)");
   tmp->setDimensionDescription(1, "Species (reduced system)");
   mScaledElasticitiesAnn = tmp;
 
   tmp = new CArrayAnnotation("Scaled concentration control coefficients", this,
                              new CCopasiMatrixInterface<CMatrix<C_FLOAT64> >(&mScaledConcCC), true);
   tmp->setMode(CArrayAnnotation::VECTOR);
   tmp->setDescription("Scaled concentration control coefficients");
   tmp->setDimensionDescription(0, "Species (reduced system)");
   tmp->setDimensionDescription(1, "Reactions (reduced system)");
   mScaledConcCCAnn = tmp;
 
   tmp = new CArrayAnnotation("Scaled flux control coefficients", this,
                              new CCopasiMatrixInterface<CMatrix<C_FLOAT64> >(&mScaledFluxCC), true);
   tmp->setMode(CArrayAnnotation::VECTOR);
   tmp->setDescription("Scaled flux control coefficients");
   tmp->setDimensionDescription(0, "Reactions (reduced system)");
   tmp->setDimensionDescription(1, "Reactions (reduced system)");
   mScaledFluxCCAnn = tmp;
 }

bool CMCAMethod::isValidProblem ( const CCopasiProblem * pProblem )

virtual

Check if the method is suitable for this problem

Returns: bool suitability of the method

Reimplemented from CCopasiMethod.

Definition at line 803 of file CMCAMethod.cpp.

References CModelEntity::ASSIGNMENT, CCopasiVector< T >::begin(), CStateTemplate::beginIndependent(), CCopasiVector< T >::end(), CStateTemplate::endIndependent(), CCopasiMessage::ERROR, CModel::getCompartments(), CCopasiProblem::getModel(), CModel::getNumIndependentReactionMetabs(), CModel::getStateTemplate(), CCopasiMethod::isValidProblem(), mpUseReeder, and mpUseSmallbone.

 {
   if (!CCopasiMethod::isValidProblem(pProblem)) return false;
 
   const CMCAProblem * pP = dynamic_cast<const CMCAProblem *>(pProblem);
 
   if (!pP)
     {
       //not a TrajectoryProblem
       CCopasiMessage(CCopasiMessage::ERROR, "Problem is not an MCA problem.");
       return false;
     }
 
   CModel * pModel = pP->getModel();
 
   if (pModel == NULL)
     return false;
 
   // Check if the model contains an ODE.
   size_t NumODE =
     pModel->getStateTemplate().endIndependent() - pModel->getStateTemplate().beginIndependent();
 
   if (pModel->getNumIndependentReactionMetabs() < NumODE)
     {
       CCopasiMessage(CCopasiMessage::ERROR, "MCA is not applicable for a system with explicit ODEs.");
       return false;
     }
 
   // Check if the model has a compartment with an assignment
   CCopasiVector< CCompartment >::const_iterator it = pModel->getCompartments().begin();
   CCopasiVector< CCompartment >::const_iterator end = pModel->getCompartments().end();
 
   for (; it != end; ++it)
     if ((*it)->getStatus() == CModelEntity::ASSIGNMENT)
       {
         CCopasiMessage(CCopasiMessage::ERROR, "MCA is not applicable for a system with changing volumes.");
         return false;
       }
 
   if (!*mpUseReeder && !*mpUseSmallbone)
     {
       CCopasiMessage(CCopasiMessage::ERROR, "At least one of the algorithm Reeder or Smallbone must be selected.");
       return false;
     }
 
   return true;
 }

C_INT32 CMCAMethod::load ( CReadConfig & configBuffer )

Read SSMCAUnscaled from configuration file

Definition at line 745 of file CMCAMethod.cpp.

References C_INT32, CReadConfig::getVariable(), CReadConfig::LOOP, and mSSReder.

 {
   C_INT32 Fail = 0;
 
   if ((Fail = configBuffer.getVariable("SSMCAUnscaled", "C_INT16",
                                        (void *) & mSSReder,
                                        CReadConfig::LOOP)))
     return Fail;
 
   return Fail;
 }

bool CMCAMethod::process ( )

virtual

Definition at line 757 of file CMCAMethod.cpp.

References CalculateMCA(), and mSteadyStateResolution.

Referenced by CMCATask::process().

 {
   // check if current state is a steady state
   // if not, calculate TimeMCA only
   if (1 /*mIsSteadyState*/)
     {
       CalculateMCA(mSteadyStateResolution);
     }
   else
     {
       //   CalculateTimeMCA(mSteadyStateResolution);
     }
 
   return true;
 }

void CMCAMethod::resizeAllMatrices ( )

virtual

Resizes all result matrices and updates the corresponding array annotations. This needs to be called before output initialization (in case the output references parts of the matrix) from the task initialization, but it also needs to be called before selecting elements of the matrices in the object selection dialog. The model needs to be set before calling this.

Definition at line 160 of file CMCAMethod.cpp.

References CModel::getMetabolitesX(), CModel::getNumDependentReactionMetabs(), CModel::getNumIndependentReactionMetabs(), CModel::getReactions(), CModel::getTotSteps(), mElasticityDependencies, mpModel, mScaledConcCC, mScaledConcCCAnn, mScaledElasticities, mScaledElasticitiesAnn, mScaledFluxCC, mScaledFluxCCAnn, mUnscaledConcCC, mUnscaledConcCCAnn, mUnscaledElasticities, mUnscaledElasticitiesAnn, mUnscaledFluxCC, mUnscaledFluxCCAnn, CMatrix< CType >::numCols(), CMatrix< CType >::numRows(), CMatrix< CType >::resize(), CArrayAnnotation::resize(), and CArrayAnnotation::setCopasiVector().

Referenced by CMCATask::updateMatrices().

 {
   if (!mpModel) return;
 
   mUnscaledElasticities.resize(mpModel->getTotSteps(),
                                mpModel->getNumIndependentReactionMetabs() + mpModel->getNumDependentReactionMetabs());
   mUnscaledElasticitiesAnn->resize();
   mUnscaledElasticitiesAnn->setCopasiVector(0, &mpModel->getReactions());
   mUnscaledElasticitiesAnn->setCopasiVector(1, &mpModel->getMetabolitesX());
 
   mUnscaledConcCC.resize(mpModel->getNumIndependentReactionMetabs() + mpModel->getNumDependentReactionMetabs(),
                          mpModel->getTotSteps());
   mUnscaledConcCCAnn->resize();
   mUnscaledConcCCAnn->setCopasiVector(0, &mpModel->getMetabolitesX());
   mUnscaledConcCCAnn->setCopasiVector(1, &mpModel->getReactions());
 
   mUnscaledFluxCC.resize(mpModel->getTotSteps(), mpModel->getTotSteps());
   mUnscaledFluxCCAnn->resize();
   mUnscaledFluxCCAnn->setCopasiVector(0, &mpModel->getReactions());
   mUnscaledFluxCCAnn->setCopasiVector(1, &mpModel->getReactions());
 
   mScaledElasticities.resize(mUnscaledElasticities.numRows(), mUnscaledElasticities.numCols());
   mScaledElasticitiesAnn->resize();
   mScaledElasticitiesAnn->setCopasiVector(0, &mpModel->getReactions());
   mScaledElasticitiesAnn->setCopasiVector(1, &mpModel->getMetabolitesX());
 
   // Reactions are columns, species are rows
   mScaledConcCC.resize(mUnscaledConcCC.numRows(), mUnscaledConcCC.numCols());
   mScaledConcCCAnn->resize();
   mScaledConcCCAnn->setCopasiVector(0, &mpModel->getMetabolitesX());
   mScaledConcCCAnn->setCopasiVector(1, &mpModel->getReactions());
 
   // Reactions are columns and rows
   mScaledFluxCC.resize(mUnscaledFluxCC.numRows(), mUnscaledFluxCC.numCols());
   mScaledFluxCCAnn->resize();
   mScaledFluxCCAnn->setCopasiVector(0, &mpModel->getReactions());
   mScaledFluxCCAnn->setCopasiVector(1, &mpModel->getReactions());
 
   mElasticityDependencies.resize(mUnscaledElasticities.numRows(), mUnscaledElasticities.numCols());
 }

bool CMCAMethod::scaleMCA	(	const bool &	status,
		C_FLOAT64	res
	)

Scales the coefficients (i.e. Kacser format, rather than Reder)

Definition at line 440 of file CMCAMethod.cpp.

References CMatrix< CType >::array(), CCopasiVector< T >::begin(), C_FLOAT64, CCopasiVector< T >::end(), CSteadyStateMethod::found, CModel::getMetabolitesX(), CModel::getNumDependentReactionMetabs(), CModel::getNumIndependentReactionMetabs(), CModel::getReactions(), CModelEntity::getValue(), mElasticityDependencies, mpModel, mScaledConcCC, mScaledElasticities, mScaledFluxCC, mSSStatus, mUnscaledConcCC, mUnscaledElasticities, mUnscaledFluxCC, and CMatrix< CType >::numCols().

Referenced by CalculateMCA().

 {
   assert(mpModel);
   // The number of metabs determined by reaction.
   size_t numSpeciesReaction =
     mpModel->getNumIndependentReactionMetabs() + mpModel->getNumDependentReactionMetabs();
   CCopasiVector< CMetab > & metabs = mpModel->getMetabolitesX();
   CCopasiVector< CMetab >::const_iterator itSpecies = metabs.begin();
   CCopasiVector< CMetab >::const_iterator endSpecies = itSpecies + numSpeciesReaction;
 
   CCopasiVector< CReaction > & reacs = mpModel->getReactions();
   CCopasiVector< CReaction >::const_iterator itReaction;
   CCopasiVector< CReaction >::const_iterator endReaction = reacs.end();
 
   size_t col;
 
   // Scale Elasticities
   C_FLOAT64 * pUnscaled;
   C_FLOAT64 * pScaled;
   bool * pElasticityDependency;
 
   // Reactions are rows, species are columns
   for (col = 0; itSpecies != endSpecies; ++itSpecies, col++)
     {
       C_FLOAT64 VolumeInv = 1.0 / (*itSpecies)->getCompartment()->getValue();
       C_FLOAT64 Number = (*itSpecies)->getValue();
 
       for (itReaction = reacs.begin(),
            pUnscaled = mUnscaledElasticities.array() + col,
            pScaled = mScaledElasticities.array() + col,
            pElasticityDependency = mElasticityDependencies.array() + col;
            itReaction != endReaction;
            ++itReaction,
            pUnscaled += numSpeciesReaction,
            pScaled += numSpeciesReaction,
            pElasticityDependency += numSpeciesReaction)
         {
           if (!*pElasticityDependency)
             {
               *pScaled = 0.0;
             }
           else if (fabs((*itReaction)->getFlux() * VolumeInv) >= res)
             {
               *pScaled = *pUnscaled * Number / (*itReaction)->getParticleFlux();
             }
           else
             {
               *pScaled = (((*itReaction)->getFlux() < 0.0) ? - std::numeric_limits<C_FLOAT64>::infinity() : std::numeric_limits<C_FLOAT64>::infinity());
             }
         }
     }
 
   //if we are not in a steady state we cannot calculate CCs
   if (mSSStatus != CSteadyStateMethod::found ||
       !status)
     {
       mScaledConcCC = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
       mScaledFluxCC = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
 
       return false;
     }
 
   // Scale ConcCC
   // Reactions are columns, species are rows
   itSpecies = metabs.begin();
   pUnscaled = mUnscaledConcCC.array();
   pScaled = mScaledConcCC.array();
 
   for (; itSpecies != endSpecies; ++itSpecies)
     {
       // In rare occasions the concentration might not be updated
       (*(*itSpecies)->getConcentrationReference()->getRefresh())();
       C_FLOAT64 alt = fabs((*itSpecies)->getConcentration());
 
       for (itReaction = reacs.begin(); itReaction != endReaction; ++itReaction, ++pUnscaled, ++pScaled)
         {
           if (alt >= res)
             *pScaled =
               *pUnscaled * (*itReaction)->getParticleFlux() / (*itSpecies)->getValue();
           else
             *pScaled = *pUnscaled * std::numeric_limits<C_FLOAT64>::infinity();
         }
     }
 
   CCopasiVector< CReaction >::const_iterator itReactionCol;
   itReaction = reacs.begin();
   pUnscaled = mUnscaledFluxCC.array();
   pScaled = mScaledFluxCC.array();
 
   for (; itReaction != endReaction; ++itReaction)
     {
       const CCompartment * pCompartment = (*itReaction)->getLargestCompartment();
       C_FLOAT64 Resolution = res * pCompartment->getValue();
 
       C_FLOAT64 Scale = (*itReaction)->getFlux();
       C_FLOAT64 tmp = 0.0;
       C_FLOAT64 eq = 0.0;
 
       // We use the summation theorem to verify the scaling
       C_FLOAT64 *pSum = pUnscaled;
       C_FLOAT64 *pSumEnd = pSum + mScaledFluxCC.numCols();
 
       for (itReactionCol = reacs.begin(); pSum != pSumEnd; ++itReactionCol, ++pSum)
         {
           tmp += *pSum * (*itReactionCol)->getFlux();
           eq += fabs(*pSum);
         }
 
       eq /= mScaledFluxCC.numCols();
 
       if (fabs(tmp) < Resolution && eq >= Resolution)
         {
           Scale = std::numeric_limits< C_FLOAT64 >::infinity();
         }
 
       tmp = 0.0;
 
       for (itReactionCol = reacs.begin(); itReactionCol != endReaction; ++itReactionCol, ++pUnscaled, ++pScaled)
         {
           // In the diagonal the scaling factors cancel.
           if (eq < res)
             {
               *pScaled = (itReactionCol != itReaction) ? 0.0  : 1.0;
             }
           else if (itReactionCol == itReaction)
             {
               *pScaled = *pUnscaled;
             }
           else if (fabs(Scale) >= res)
             {
               *pScaled = *pUnscaled * (*itReactionCol)->getFlux() / Scale;
             }
           else
             {
               const CCompartment * pCompartmentCol = (*itReactionCol)->getLargestCompartment();
               C_FLOAT64 ScaleCol = (pCompartmentCol == NULL) ?
                                    fabs((*itReactionCol)->getFlux()) :
                                    fabs((*itReactionCol)->getFlux() / pCompartmentCol->getValue());
 
               if (fabs(ScaleCol) <= res)
                 {
                   *pScaled = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
                 }
               else
                 {
                   *pScaled = (((*itReaction)->getFlux() < 0.0) ?
                               - std::numeric_limits<C_FLOAT64>::infinity() :
                               std::numeric_limits<C_FLOAT64>::infinity());
                 }
             }
 
           tmp += *pScaled;
         }
     }
 
   return true;
 }

void CMCAMethod::setFactor ( C_FLOAT64 factor )

Definition at line 787 of file CMCAMethod.cpp.

References mFactor.

 {
   this->mFactor = factor;
 }

void CMCAMethod::setModel ( CModel * model )

Set the Model

Definition at line 663 of file CMCAMethod.cpp.

References mpModel.

Referenced by CMCATask::printResult(), and CMCATask::updateMatrices().

 {
   mpModel = model;
 }

void CMCAMethod::setSteadyStateResolution ( C_FLOAT64 factor )

Definition at line 792 of file CMCAMethod.cpp.

References mSteadyStateResolution.

 {
   this->mSteadyStateResolution = resolution;
 }

void CMCAMethod::setSteadyStateTask ( CSteadyStateTask * pSteadyStateTask )

Definition at line 773 of file CMCAMethod.cpp.

References CSteadyStateTask::getResult(), mpSteadyStateTask, mSSStatus, and CSteadyStateMethod::notFound.

Referenced by CMCATask::process().

 {
   mpSteadyStateTask = pSteadyStateTask;
 
   if (mpSteadyStateTask != NULL)
     {
       mSSStatus = mpSteadyStateTask->getResult();
     }
   else
     {
       mSSStatus = CSteadyStateMethod::notFound;
     }
 }