COPASI API  4.16.103
CEFMAlgorithm.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) 2002 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 /**
16  * CEFMAlgorithm class.
17  * Used to calculate elementary flux modes
18  *
19  * Created for COPASI by Stefan Hoops 2002-05-08
20  * (C) Stefan Hoops 2002
21  */
22 
23 #include "copasi.h"
24 
25 #include "CEFMAlgorithm.h"
26 #include "CEFMProblem.h"
27 #include "CEFMTask.h"
28 #include "CFluxMode.h"
29 #include "CTableauMatrix.h"
30 
31 #include "model/CModel.h"
34 
36  CCopasiNode< size_t >(),
37  mTableauLines()
38 {}
39 
41  CCopasiNode< size_t >(src),
42  mTableauLines(src.mTableauLines)
43 {}
44 
46  const CTableauMatrix & matrix):
47  CCopasiNode< size_t >(index),
48  mTableauLines()
49 {
50  update(matrix);
51 }
52 
54 {}
55 
57 {
58  mTableauLines.clear();
59 
60  std::list< const CTableauLine * >::const_iterator it = matrix.begin();
61  std::list< const CTableauLine * >::const_iterator end = matrix.end();
62  size_t TableauLineIndex = 0;
63 
64  for (; it != end; ++it, ++TableauLineIndex)
65  {
66  if ((*it)->getMultiplier(TableauLineIndex) != 0.0)
67  {
68  mTableauLines.push_back(TableauLineIndex);
69  }
70  }
71 }
72 
74  CEFMMethod(CCopasiTask::fluxMode, CCopasiMethod::EFMAlgorithm, pParent),
75  mpModel(NULL),
76  mStoi(),
77  mReversible(0),
78  mpCurrentTableau(NULL),
79  mpNextTableau(NULL),
80  mIndexSet()
81 {initObjects();}
82 
84  CEFMMethod(CCopasiTask::fluxMode, subType, pParent),
85  mpModel(NULL),
86  mStoi(),
87  mReversible(0),
88  mpCurrentTableau(NULL),
89  mpNextTableau(NULL),
90  mIndexSet()
91 {initObjects();}
92 
94  const CCopasiContainer * pParent):
95  CEFMMethod(src, pParent),
96  mpModel(NULL),
97  mStoi(),
98  mReversible(0),
99  mpCurrentTableau(NULL),
100  mpNextTableau(NULL),
101  mIndexSet()
102 {initObjects();}
103 
105 {
109 }
110 
112 {
114 }
115 
117 {
118  if (!CEFMMethod::initialize())
119  {
120  return false;
121  }
122 
123  CEFMTask * pTask = dynamic_cast< CEFMTask *>(getObjectParent());
124 
125  if (pTask == NULL) return false;
126 
127  mpModel = pTask->getProblem()->getModel();
128 
129  if (mpModel == NULL) return false;
130 
131  mpFluxModes->clear();
132 
133  /* ModelStoi is the transpose of the models stoichiometry matrix */
135 
136  size_t row, numRows = ModelStoi.numRows();
137  size_t col, numCols = ModelStoi.numCols();
138 
139  /* Size the stoichiometry matrix passed to the algorithm */
140  mStoi.resize(numRows);
141  std::vector< std::vector< C_FLOAT64 > >::iterator it = mStoi.begin();
142  std::vector< std::vector< C_FLOAT64 > >::iterator end = mStoi.end();
143 
144  for (; it != end; ++it)
145  it->resize(numCols);
146 
147  /* Get the reactions from the model */
148  /* Note: We have as many reactions as we have rows in ModelStoi */
149  const CCopasiVectorNS < CReaction > & Reaction = mpModel->getReactions();
150 
151  /* Vector to keep track of the rearrangements necessary to put the */
152  /* reversible reactions to the top of stoichiometry matrix */
153  mpReorderedReactions->resize(numRows);
154 
155  /* Reversible reaction counter */
156  mReversible = 0;
157  size_t Insert;
158  size_t InsertReversible = 0;
159  size_t InsertIrreversible = numRows - 1;
160 
161  /* Build the transpose of the stoichiometry matrix, */
162 
163  /* sort reversible reactions to the top, count them */
164 
165  /* and keep track of the rearrangement */
166  for (row = 0; row < numRows; row++)
167  {
168  if (Reaction[row]->isReversible())
169  {
170  Insert = InsertReversible++;
171  mReversible++;
172  }
173  else
174  Insert = InsertIrreversible--;
175 
176  (*mpReorderedReactions)[Insert] = Reaction[row];
177 
178  for (col = 0; col < numCols; col++)
179  mStoi[Insert][col] = ModelStoi(row, col);
180  }
181 
182  mStep = 0;
183  mStepProcess = 0;
184  mMaxStep = (unsigned C_INT32) numCols;
185 
186  if (mpCallBack)
187  mhSteps =
188  mpCallBack->addItem("Current Step",
189  mStepProcess,
190  & mMaxStep);
191 
192  return true;
193 }
194 
196 {
197  if (!initialize())
198  {
199  if (mpCallBack)
201 
202  return false;
203  }
204 
206 
207  return true;
208 }
209 
211 {
212  bool Continue = true;
213 
214  if (mStoi.size())
215  {
216  /* initialize the current tableau matrix */
219 
220  /* Do the iteration */
221  mIndexSet.resize(mMaxStep);
222 
223  for (mStep = 0; mStep < mMaxStep; mStep++)
224  mIndexSet[mStep] = mStep;
225 
226  while (findMinimalCombinationIndex() && Continue)
227  {
229  mStepProcess++;
230 
231  if (mpCallBack)
232  Continue &= mpCallBack->progressItem(mhSteps);
233 
234  static_cast<CCopasiTask *>(getObjectParent())->output(COutputInterface::DURING);
235  }
236 
237  /* Build the elementary flux modes to be returned */
238  if (Continue)
239  buildFluxModes();
240 
241  /* Delete the current / final tableau matrix */
243  }
244 
245  if (mpCallBack)
246  Continue &= mpCallBack->finishItem(mhSteps);
247 }
248 
250 {
251  std::list< const CTableauLine * >::iterator a;
252  std::list< const CTableauLine * >::iterator b;
253  C_FLOAT64 ma, mb;
254 #ifdef COPASI_DEBUG_TRACE
255  DebugFile << *mpCurrentTableau << std::endl;
256 #endif //COPASI_DEBUG_TRACE
257 
259 
260  /* Move all lines with zeros in the step column to the new tableau matrix */
261  /* and remove them from the current tableau matrix */
262  a = mpCurrentTableau->begin();
263 
264  bool Continue = true;
265 
266  unsigned C_INT32 Counter, MaxCounter;
267  size_t hCounter;
268 
269  Counter = 0;
270  MaxCounter = (unsigned C_INT32) mpCurrentTableau->size();
271 
272  if (mpCallBack)
273  hCounter =
274  mpCallBack->addItem("Current Line",
275  Counter,
276  & MaxCounter);
277 
278  while (a != mpCurrentTableau->end() && Continue)
279  if ((*a)->getMultiplier(mStep) == 0.0)
280  {
281  /* We have to make sure that "a" points to the next element in the */
282  /* list after the removal of itself */
283 
284  if (a == mpCurrentTableau->begin())
285  {
286  mpNextTableau->addLine(*a, false);
288  a = mpCurrentTableau->begin();
289  }
290  else
291  {
292  /* We have to remember the previous element so that a++ points to */
293  /* past the removed one */
294  b = a--;
295  mpNextTableau->addLine(*b, false);
297  a++;
298  }
299 
300  Counter++;
301 
302  if (mpCallBack)
303  Continue &= mpCallBack->progressItem(hCounter);
304  }
305  else
306  a++;
307 
308  C_FLOAT64 Sign;
309 
310  /* Now we create all linear combinations of the remaining lines in the */
311  /* current tableau */
312  a = mpCurrentTableau->begin();
313 
314  //std::cout << "Tableau Marker A" << std::endl << std::endl;
315 
316  while (a != mpCurrentTableau->end() && Continue)
317  {
318  b = a;
319  b++;
320 
321  /* We make sure that "mb" is positive */
322  mb = (*a)->getMultiplier(mStep);
323 
324  if (mb < 0.0)
325  {
326  mb *= -1.0;
327  Sign = 1.0;
328  }
329  else
330  Sign = -1.0;
331 
332  //std::cout << "Tableau Marker B" << std::endl << std::endl;
333 
334  while (b != mpCurrentTableau->end() && Continue)
335  {
336  ma = Sign * (*b)->getMultiplier(mStep);
337 
338  /* The multiplier "ma" for irreversible reactions must be positive */
339  if (ma > 0.0 || (*a)->isReversible())
340  mpNextTableau->addLine(new CTableauLine(ma, **a, mb, **b));
341 
342  /*CTableauLine * debugLine = new CTableauLine(ma, **a, mb, **b);
343  if(debugLine->isReversible())
344  std::cout << "Reversible Rxn" << std::endl;
345  else std::cout << "Irreversible Rxn" << std::endl;
346  std::cout << "Flux Score: " << debugLine->getScore() << "Flux Mode Vector: "
347  << debugLine->getFluxMode() << std::endl;*/
348 
349  b++;
350 
351  if (mpCallBack)
352  Continue &= mpCallBack->proceed();
353  }
354 
355  // We no longer need a since all linear combinations have been build;
357  a = mpCurrentTableau->begin();
358 
359  Counter++;
360 
361  if (mpCallBack)
362  Continue &= mpCallBack->progressItem(hCounter);
363  }
364 
365  if (mpCallBack)
366  Continue &= mpCallBack->finishItem(hCounter);
367 
368  /* Assign the next tableau to the current tableau and cleanup */
370 
372 
373  mpNextTableau = NULL;
374 }
375 
377 {
378  mpFluxModes->clear();
379 
380  std::list< const CTableauLine * >::iterator a = mpCurrentTableau->begin();
381  std::list< const CTableauLine * >::iterator end = mpCurrentTableau->end();
382 
383  while (a != end)
384  {
385  mpFluxModes->push_back(CFluxMode(*a));
386  a++;
387  }
388 }
389 
391 {
392  double minCombine = std::numeric_limits<double>::infinity();
393  double combine = 0;
394  size_t minIndex = 0;
395  size_t counter;
396 
397  if (mIndexSet.size() == 0)
398  return false;
399  else if (mIndexSet.size() == 1)
400  {
401  mStep = (unsigned C_INT32) mIndexSet[0];
402  mIndexSet.pop_back();
403  return true;
404  }
405 
406  for (counter = 0; counter < mIndexSet.size(); counter++)
407  {
408  combine = calculateCombinations(mIndexSet[counter]);
409 
410  if (combine < minCombine)
411  {
412  minCombine = combine;
413  minIndex = counter;
414  }
415 
416  if (combine == 0)
417  break;
418  }
419 
420  mStep = (unsigned C_INT32) mIndexSet[minIndex];
421  mIndexSet.erase(mIndexSet.begin() + minIndex);
422 
423  return true;
424 }
425 
427 {
428  double posIrr = 0;
429  double negIrr = 0;
430  double rev = 0;
431 
432  //Reversible reactions
433  std::list< const CTableauLine * >::const_iterator it = mpCurrentTableau->begin();
434  std::list< const CTableauLine * >::const_iterator end = mpCurrentTableau->end();
435 
436  for (; it != end; ++it)
437  {
438  if ((*it)->isReversible() && (*it)->getMultiplier(index) != 0.0)
439  {
440  rev++;
441  }
442  else if ((*it)->getMultiplier(index) < 0.0)
443  {
444  negIrr++;
445  }
446  else if ((*it)->getMultiplier(index) > 0.0)
447  {
448  posIrr++;
449  }
450  }
451 
452  return (posIrr + rev) * (negIrr + rev);
453 }
double calculateCombinations(size_t index)
#define pdelete(p)
Definition: copasi.h:215
CCopasiProblem * getProblem()
unsigned C_INT32 mMaxStep
std::vector< std::vector< C_FLOAT64 > > mStoi
CEFMAlgorithm(const CCopasiContainer *pParent=NULL)
std::list< const CTableauLine * >::iterator begin()
size_t size() const
std::vector< CFluxMode > * mpFluxModes
Definition: CEFMMethod.h:91
std::list< const CTableauLine * >::iterator end()
CTableauMatrix * mpNextTableau
#define C_INT32
Definition: copasi.h:90
virtual bool calculate()
virtual bool progressItem(const size_t &handle)
virtual bool initialize()
Definition: CEFMMethod.cpp:109
CTSSATask * pTask
virtual bool proceed()
#define DESTRUCTOR_TRACE
Definition: copasi.h:206
bool findMinimalCombinationIndex()
std::vector< size_t > mIndexSet
void calculateFluxModes()
virtual bool initialize()
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
std::vector< const CReaction * > * mpReorderedReactions
Definition: CEFMMethod.h:96
void removeLine(const std::list< const CTableauLine * >::iterator line)
size_t numRows() const
Definition: CMatrix.h:767
virtual bool finishItem(const size_t &handle)
void update(const CTableauMatrix &matrix)
void addLine(const CTableauLine *src, const bool &check=true)
CTableauMatrix * mpCurrentTableau
void calculateNextTableau()
CModel * mpModel
#define C_FLOAT64
Definition: copasi.h:92
unsigned C_INT32 mStep
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
unsigned C_INT32 mStepProcess
const CMatrix< C_FLOAT64 > & getStoi() const
Definition: CModel.cpp:1160
size_t mReversible
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
CProcessReport * mpCallBack
CModel * getModel() const
CCopasiContainer * getObjectParent() const