COPASI API  4.16.103
CSSAMethod.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2013 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 /**
12  * CSSAMethod class.
13  * Methods to conduct Stoichiometric Stability Analysis.-
14  */
15 
16 #include "copasi.h"
17 
18 #include "lapack/blaswrap.h"
19 #include "lapack/lapackwrap.h"
20 
26 
27 #include "model/CModel.h"
28 #include "model/CState.h"
30 
31 /**
32  * Copy constructor.
33  * @param "const CSSAMethod &" src
34  */
36  const CCopasiContainer * pParent):
37  CEFMAlgorithm(src, pParent)
38 {}
39 
41  CEFMAlgorithm(CCopasiMethod::stoichiometricStabilityAnalysis, pParent)
42 {}
43 
44 // taken from CEFMAlgorithm, modified to make reactions reversible, implicitely,
45 // i.e. without changing the model itself.
47 {
48  CEFMTask * pTask = dynamic_cast< CEFMTask *>(getObjectParent());
49 
50  if (pTask == NULL) return false;
51 
52  mpModel = pTask->getProblem()->getModel();
53 
54  if (mpModel == NULL) return false;
55 
56  mFluxModes.clear();
58 
59  /* ModelStoi is the transpose of the models stoichiometry matrix */
61 
62  unsigned C_INT32 row, numRows = ModelStoi.numRows();
63  unsigned C_INT32 col, numCols = ModelStoi.numCols();
64 
65  /* Get the reactions from the model */
66  /* Note: We have as many reactions as we have rows in ModelStoi */
68 
69  /* Reversible reaction counter */
70  mReversible = 0;
71 
72  for (row = 0; row < numRows; row++)
73  {
74  if (Reaction[row]->isReversible())
75  mReversible++;
76  }
77 
78  unsigned C_INT32 numOrigRows = numRows;
79  numRows += mReversible;
80  mNumReactions = numRows;
81 
82  mReversible = 0;
83  /* Size the stoichiometry matrix passed to the algorithm */
84  mStoi.resize(numRows);
85  std::vector< std::vector< C_FLOAT64 > >::iterator it = mStoi.begin();
86  std::vector< std::vector< C_FLOAT64 > >::iterator end = mStoi.end();
87 
88  for (; it != end; ++it)
89  it->resize(numCols);
90 
91  /* Vector to keep track of the rearrangements necessary to put the */
92  /* reversible reactions to the top of stoichiometry matrix */
93  mReorderedReactions.resize(numRows);
94  mIsBackwardReaction.resize(numRows);
95 
96  unsigned C_INT32 Insert;
97  unsigned C_INT32 InsertReversible = 0;
98  unsigned C_INT32 InsertIrreversible = numRows - 1;
99 
100  /* Build the transpose of the stoichiometry matrix, */
101 
102  /* sort reversible reactions to the top, count them */
103 
104  /* and keep track of the rearrangement */
105 
106  //for (row = 0; row < numOrigRows; row++)
107  for (row = 0; row < numOrigRows; row++)
108  {
109  if (Reaction[row]->isReversible())
110  {
111  Insert = InsertReversible++;
112  mIsBackwardReaction[Insert] = true;
113 
114  mReversible++;
115  }
116  else
117  {
118  Insert = InsertIrreversible--;
119  mIsBackwardReaction[Insert] = false;
120  }
121 
122  mReorderedReactions[Insert] = Reaction[row];
123 
124  for (col = 0; col < numCols; col++)
125  {
126  mStoi[Insert][col] = ModelStoi(row, col);
127 
128  if (Reaction[row]->isReversible())
129  mStoi[InsertReversible - 1][col] = - ModelStoi(row, col);
130  }
131  }
132 
133  mStep = 0;
134  mMaxStep = numCols;
135 
136  if (mpCallBack)
137  mhSteps =
138  mpCallBack->addItem("Current Step",
140  & mStep,
141  & mMaxStep);
142 
143  return true;
144 }
145 
147 {
148  return false;
149 }
150 
151 bool
153 {
154  if (!initialize())
155  {
156  if (mpCallBack)
158 
159  return false;
160  }
161 
162  if (!mpModel)
163  {
164  if (mpCallBack)
166 
167  return false;
168  }
169 
170  //initialize matrices for calculation;
171  // if (mpModel->compileIfNecessary(NULL))
172  // {
176  //}
177  // else
178  // return false;
179 
180  if (!testForMixingStability())
181  {
182  if (mpCallBack)
184 
185  return false;
186  }
187 
188  if (mpCallBack)
190 
191  return true;
192 }
193 
194 TriLogic
196 {return mIsMixingStable[indexEC];}
197 
198 bool
200 {return mIsBackwardReaction[index];}
201 
202 bool
204 {
206 
207  C_INT32 numECs = mExtremeCurrents.size();
208 
209  C_INT32 ijmax = mTransformedSubJacobians[0].numRows();
210 
211  mIsMixingStable.clear();
212 
213  C_FLOAT64 outarray[ijmax * ijmax];
214  memset(outarray, 0, ijmax * ijmax * sizeof(C_FLOAT64));
215 
216  for (int ecIndex = 0; ecIndex < numECs; ++ecIndex)
217  {
218  C_FLOAT64 *inarray = mTransformedSubJacobians[ecIndex].array();
219 
220  // add matrix and its extremeCurrent - upper triangle only
221  for (int i = 0; i < ijmax; ++i)
222  for (int j = i; j < ijmax; ++j)
223  outarray[i + j * ijmax] = inarray[i + j * ijmax] + inarray[j + i * ijmax];
224 
225  // get the eigenvalues AND eigenvectors:
226  // input to dsyev_()
227  char jobz = 'V';
228  char uplo = 'U';
229  C_INT32 n = ijmax;
230  C_FLOAT64 *A = outarray;
231  C_INT32 lda = ijmax;
232  C_INT32 lwork = 64 * ijmax;
233  // dsyev_() output
234  C_FLOAT64 eigenvalues[ijmax];
235  C_FLOAT64 work[lwork];
236  C_INT32 info;
237 
238  dsyev_(&jobz, &uplo, &n, A, &lda, eigenvalues, work, &lwork, &info);
239 
240  if (info)
241  return false; // maybe raise an exception here, to discriminate from the sane case
242 
243  CVector<C_FLOAT64> EC = mExtremeCurrents[ecIndex];
244 
245  // set the state to 'stable' (1 - stable; 0 - unknown; -1 - not mixing stable)
246  TriLogic state = TriTrue;
247 
248  int partMetabs[ijmax];
249  memset(&partMetabs, 0, ijmax * sizeof(int));
250 
251  for (int j = 0; j < EC.size(); ++j)
252  {
253  if (EC[j])
254  {
255  // if one reaction's kinetic function is not Mass Action Law,
256  // this method is not valid and we can say nothing definite
257  // about the stability of this EC.
258  std::string type = mReorderedReactions[j]->getFunction()->getObjectName();
259 
260  if (type.find("Mass action") == std::string::npos &&
261  type.find("Constant flux") == std::string::npos)
262  {
263  state = TriUnspecified;
264  break;
265  }
266 
267  for (int i = 0; i < ijmax; ++i)
268  {
269  if (mStoi[j][i])
270  {
271  partMetabs[i] = 1;
272  break;
273  }
274  }
275  }
276  }
277 
278  if (state == TriTrue)
279  {
280  // eigenvalues are sorted in rising order, and if the last ev is
281  // not positive, we have still to look after zero evs.
282  if (eigenvalues[ijmax - 1] > 0)
283  state = TriFalse;
284  else
285  {
286  int zeroes = 0;
287 
288  for (int i = 0; i < ijmax && state == 1; ++i)
289  {
290  // only count zero eigenvals if the eigenvector contains a reaction
291  // that is participating in this extreme current.
292  if (eigenvalues[i] == 0)
293  {
294  for (int j = 0; j < ijmax; ++ j)
295  if (partMetabs[j] && outarray[i * ijmax + j])
296  if (++zeroes)
297  break;
298 
299  if (zeroes == 2)
300  state = TriFalse;
301  }
302  }
303  }
304  }
305 
306  mIsMixingStable.push_back(state);
307  }
308 
309  return true;
310 }
311 
312 bool
314 {
315  mTransformedSubJacobians.clear();
316 
317  // parameters for the double matrix multiplication with dgemm_
318  char cN = 'N'; // don't transpose matrices
319  C_INT m = mStoichiometry.numRows(); // the leading dimension of the Stoi.M. and of the resulting M.
320  C_INT n = mStoichiometry.numCols(); // the dimensions of matrix diagE, also used as 'k'.
321 
322  C_FLOAT64 alpha = 1.;
323  C_FLOAT64 beta = 0.;
324 
325  CMatrix< C_FLOAT64 > product;
326  product.resize(m, m);
327 
328  std::vector< CVector<C_FLOAT64> >::iterator iter;
329 
330  for (iter = mExtremeCurrents.begin(); iter != mExtremeCurrents.end(); ++iter)
331  {
332  CMatrix<C_FLOAT64> diagE = diag(*iter);
333 
334  // compute N*Ei*kappa^T:
335  // As the representation of the matrices in their arrays is in
336  // transposed (mathematical) order, we compute kappa*Ei*N^T,
337  // which corresponds to
338  // (mTransposedKineticMatrix*diagE)*mStoichiometry
339  C_FLOAT64 dummy[m * n];
340  dgemm_(&cN, &cN,
341  &m, &n, &n, &alpha,
343  diagE.array(), &n,
344  &beta, dummy, &m);
345  dgemm_(&cN, &cN,
346  &m, &m, &n, &alpha,
347  dummy, &m,
348  mStoichiometry.array(), &n,
349  &beta, product.array(), &m);
350 
351  mTransformedSubJacobians.push_back(product);
352  }
353 
354  return true;
355 }
356 
357 bool
359 {
360  if (!mpModel) return false;
361 
362  //rebuild the stoichiometric matrix with reversible reactions divided into two.
363 
364  // build the _transposed_ of stoichiometry, with expanded reversible reactions
365  int numRows = mStoi[0].size();
366  int numCols = mStoi.size();
367 
368  mStoichiometry.resize(numRows, numCols);
369  memset(mStoichiometry.array(), 0, numRows * numCols * sizeof(C_FLOAT64));
370 
371  for (int i = 0; i < numRows; ++i)
372  for (int j = 0; j < numCols; ++j)
373  mStoichiometry(i, j) = mStoi[j][i];
374 
375  return true;
376 }
377 
378 bool
380 {
381  if (!mpModel) return false;
382 
384 
385  C_FLOAT64 num_reactions = mNumReactions;
386 
387  mTransposedKineticMatrix.resize((C_INT64)num_reactions, (C_INT64)num_metabolites);
388  memset(mTransposedKineticMatrix.array(), 0, (int)(num_metabolites * num_reactions * sizeof(C_FLOAT64)));
389 
391 
392  for (int ireac = 0; ireac < mIsBackwardReaction.size(); ++ireac)
393  {
395 
396  if (mIsBackwardReaction[ireac])
397  substrates = mReorderedReactions[ireac]->getChemEq().getProducts();
398  else
399  substrates = mReorderedReactions[ireac]->getChemEq().getSubstrates();
400 
401  for (int jsubst = 0; jsubst < substrates.size(); ++jsubst)
402  {
403  CCopasiVector< CMetab >::iterator found = find(metabOrder.begin(),
404  metabOrder.end(),
405  substrates[jsubst]->getMetabolite());
406 
407  if ((*found)->isDependent()) continue;
408 
409  int kmetab = std::distance(metabOrder.begin(), found);
410  mTransposedKineticMatrix[ireac][kmetab] = substrates[jsubst]->getMultiplicity();
411  }
412  }
413 
414  return true;
415 }
416 
417 bool
419 {
420  mExtremeCurrents.clear();
421 
423 
424  std::vector< CFluxMode > fluxmodes = getFluxModes();
425 
426  CVector< C_FLOAT64 > extremecurrent;
427  extremecurrent.resize(mNumReactions);
428 
429  std::vector< CFluxMode >::iterator iter;
430 
431  for (iter = fluxmodes.begin(); iter != fluxmodes.end(); ++iter)
432  {
433  memset(extremecurrent.array(), 0, extremecurrent.size()*sizeof(C_FLOAT64));
434  unsigned C_INT32 i;
435 
436  for (i = 0; i < iter->size(); ++i)
437  {
438  extremecurrent[iter->getReactionIndex(i)] = iter->getMultiplier(i);
439  }
440  }
441 
442  return true;
443 }
bool process(CProcessReport *handler)
Definition: CSSAMethod.cpp:146
#define C_INT
Definition: copasi.h:115
std::vector< bool > mIsBackwardReaction
Definition: CSSAMethod.h:62
bool decomposeJacobian()
Definition: CSSAMethod.cpp:313
CCopasiProblem * getProblem()
std::vector< TriLogic > mIsMixingStable
Definition: CSSAMethod.h:79
unsigned C_INT32 mMaxStep
virtual size_t size() const
std::vector< std::vector< C_FLOAT64 > > mStoi
bool isReactionReversed(int index) const
Definition: CSSAMethod.cpp:199
iterator begin()
virtual size_t numRows() const
Definition: CMatrix.h:138
TriLogic
Definition: copasi.h:125
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
#define C_INT64
Definition: copasi.h:88
#define C_INT32
Definition: copasi.h:90
const CMatrix< C_FLOAT64 > & getRedStoi() const
Definition: CModel.cpp:1154
std::vector< CVector< C_FLOAT64 > > mExtremeCurrents
Definition: CSSAMethod.h:74
bool testForMixingStability()
Definition: CSSAMethod.cpp:203
CMatrix< C_FLOAT64 > mTransposedKineticMatrix
Definition: CSSAMethod.h:70
CTSSATask * pTask
iterator end()
bool buildKineticMatrix()
Definition: CSSAMethod.cpp:379
void calculateFluxModes()
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
virtual void cleanup()
TriLogic isMixingStable(int indexEC)
Definition: CSSAMethod.cpp:195
bool calculate()
Definition: CSSAMethod.cpp:152
bool initialize()
Definition: CSSAMethod.cpp:46
virtual bool finishItem(const size_t &handle)
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 buildExtremeCurrents()
Definition: CSSAMethod.cpp:418
CCopasiVectorNS< CReaction > mReactions
Definition: CSSAMethod.h:53
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
std::vector< CMatrix< C_FLOAT64 > > mTransformedSubJacobians
Definition: CSSAMethod.h:72
bool buildStoichiometry()
Definition: CSSAMethod.cpp:358
CMatrix< CType > diag(CVector< CType > vector)
Definition: CSSAMethod.h:103
size_t size() const
Definition: CVector.h:100
int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a, integer *lda, doublereal *w, doublereal *work, integer *lwork, integer *info)
CModel * mpModel
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
CMatrix< C_FLOAT64 > mStoichiometry
Definition: CSSAMethod.h:56
unsigned C_INT32 mStep
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
C_INT32 mNumReactions
Definition: CSSAMethod.h:65
CSSAMethod(const CSSAMethod &src, const CCopasiContainer *pParent=NULL)
Definition: CSSAMethod.cpp:35
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
virtual size_t numCols() const
Definition: CMatrix.h:144
size_t mReversible
CProcessReport * mpCallBack
virtual CType * array()
Definition: CMatrix.h:337
CModel * getModel() const
CCopasiContainer * getObjectParent() const