COPASI API  4.16.103
CBitPatternMethod.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 #include <stdlib.h>
7 #include <cmath>
8 
9 #include "copasi.h"
10 
11 #include "CBitPatternMethod.h"
12 #include "CEFMProblem.h"
13 #include "CEFMTask.h"
14 #include "CStepMatrix.h"
15 #include "CBitPatternTree.h"
16 
17 #include "model/CModel.h"
18 #include "model/CChemEqInterface.h"
21 
22 #include "lapack/blaswrap.h"
23 #include "lapack/lapackwrap.h"
24 
25 #define DEBUG_MATRIX
26 
28  CEFMMethod(CCopasiTask::fluxMode, CCopasiMethod::EFMBitPatternAlgorithm, pParent),
29  mpModel(NULL),
30  mProgressCounter(0),
31  mProgressCounterMax(0),
32  mhProgressCounter(0),
33  mProgressCounter2(0),
34  mProgressCounter2Max(0),
35  mhProgressCounter2(0),
36  mReactionForward(),
37  mReactionPivot(0),
38  mExpandedStoiTranspose(0, 0),
39  mpStepMatrix(NULL),
40  mMinimumSetSize(0),
41  mStep(0),
42  mContinueCombination(true)
43 {
44  initObjects();
45 }
46 
48  const CCopasiContainer * pParent):
49  CEFMMethod(CCopasiTask::fluxMode, subType, pParent),
50  mpModel(NULL),
51  mProgressCounter(0),
52  mProgressCounterMax(0),
53  mhProgressCounter(0),
54  mProgressCounter2(0),
55  mProgressCounter2Max(0),
56  mhProgressCounter2(0),
57  mReactionForward(),
58  mReactionPivot(0),
59  mExpandedStoiTranspose(0, 0),
60  mpStepMatrix(NULL),
61  mMinimumSetSize(0),
62  mStep(0),
63  mContinueCombination(true)
64 {
65  initObjects();
66 }
67 
69  const CCopasiContainer * pParent):
70  CEFMMethod(src, pParent),
71  mpModel(src.mpModel),
72  mProgressCounter(src.mProgressCounter),
73  mProgressCounterMax(src.mProgressCounterMax),
74  mhProgressCounter(src.mhProgressCounter),
75  mProgressCounter2(src.mProgressCounter2),
76  mProgressCounter2Max(src.mProgressCounter2Max),
77  mhProgressCounter2(src.mhProgressCounter2),
78  mReactionForward(src.mReactionForward),
79  mReactionPivot(src.mReactionPivot),
80  mExpandedStoiTranspose(src.mExpandedStoiTranspose),
81  mpStepMatrix(src.mpStepMatrix),
82  mMinimumSetSize(src.mMinimumSetSize),
83  mStep(src.mStep),
84  mContinueCombination(src.mContinueCombination)
85 {
86  initObjects();
87 }
88 
90 {
91 }
92 
94 {
96 }
97 
99 {
100  if (!CEFMMethod::initialize())
101  {
102  return false;
103  }
104 
106  mReactionForward.clear();
107 
108  mContinueCombination = true;
109 
110  CEFMTask * pTask = dynamic_cast< CEFMTask *>(getObjectParent());
111 
112  if (pTask == NULL) return false;
113 
114  mpModel = pTask->getProblem()->getModel();
115 
116  if (mpModel == NULL) return false;
117 
118  // We first build the kernel matrix
119  CMatrix< C_INT64 > KernelMatrix;
120  buildKernelMatrix(KernelMatrix);
121 
122  mMinimumSetSize = KernelMatrix.numCols() - 2;
123 
124  std::stack< CStepMatrixColumn * > KernelStack;
125  CStepMatrixColumn * KernelCol;
126 
127  for (unsigned int col = 0; col < KernelMatrix.numCols(); col++)
128  {
129  KernelCol = new CStepMatrixColumn(KernelMatrix.numRows());
130 
131  for (unsigned int row = 0; row < KernelMatrix.numRows(); row++)
132  {
133  KernelCol->push_front(KernelMatrix(row, col));
134  }
135 
136  KernelStack.push(KernelCol);
137  }
138 
139  // Now we create the initial step matrix
140  mpStepMatrix = new CStepMatrix(KernelMatrix);
141  /*mpStepMatrix = new CStepMatrix(KernelMatrix.numRows());
142  buildFluxModeMatrix(mpStepMatrix, KernelStack);
143 
144  DebugFile << *mpStepMatrix << std::endl;*/
145 
146  /*
147  mpStepMatrix = FluxModeMatrix;
148  mpStepMatrix.convertMatrix();
149  */
150 
151  //mpStepMatrix = new CStepMatrix(KernelMatrix);
152 
153  mProgressCounter = 0;
155 
156  if (mpCallBack)
158  mpCallBack->addItem("Current Step",
161 
162  return true;
163 }
164 
166 {
167  bool Continue = true;
168 
169  if (!initialize())
170  {
171  if (mpCallBack)
173 
174  return false;
175  }
176 
177  while (mpStepMatrix->getNumUnconvertedRows() > 0 &&
178  Continue)
179  {
181 
182  std::vector< CStepMatrixColumn * > PositiveColumns;
183  std::vector< CStepMatrixColumn * > NegativeColumns;
184  std::vector< CStepMatrixColumn * > NullColumns;
185 
186  if (mpStepMatrix->splitColumns(PositiveColumns,
187  NegativeColumns,
188  NullColumns))
189  {
190  // Process each step.
191  // Iterate over all combinations and add/remove columns to the step matrix
192  mProgressCounter2 = 0;
193  mProgressCounter2Max = (unsigned C_INT32)(PositiveColumns.size() * NegativeColumns.size());
194 
195  if (mpCallBack)
197  mpCallBack->addItem("Combinations",
200 
201  for (unsigned int i = 0; i < NegativeColumns.size(); i++)
202  for (unsigned int j = 0; j < PositiveColumns.size(); j++)
203  combine(PositiveColumns[j], NegativeColumns[i], NullColumns);
204 
205  if (mpCallBack)
207 
208  Continue &= mContinueCombination;
209 
210  if (Continue)
211  {
212  // We can now destroy all negative columns, which removes them from
213  // the step matrix.
214  mpStepMatrix->removeInvalidColumns(NegativeColumns);
215 
216  // Remove columns of the step matrix which are no longer extreme rays.
217  findRemoveInvalidColumns(NullColumns);
218 
219  // We compact the step matrix which has empty columns due to the removal of columns above.
221 
222  // Now we can convert the processed row.
224  }
225  }
226 
228 
229  if (mpCallBack)
231  }
232 
233  if (Continue)
234  {
235  buildFluxModes();
236  }
237 
238  if (mpCallBack)
240 
241  return true;
242 }
243 
244 //TODO: Change combine method to accept columns instead of trees.
246  const CStepMatrixColumn * pNegative,
247  const std::vector< CStepMatrixColumn * > NullColumns)
248 {
250  {
252  }
253 
255  {
256  return;
257  }
258 
259  CZeroSet Intersection = CZeroSet::intersection(pPositive->getZeroSet(),
260  pNegative->getZeroSet());
261 
262  // Adjacency test
263  if (Intersection.getNumberOfSetBits() < mMinimumSetSize)
264  {
265  return;
266  }
267 
268  //INSERT RANK TEST HERE
269 
270  CStepMatrixColumn * RankTestColumn = new CStepMatrixColumn(Intersection, pPositive, pNegative);
271  CMatrix<C_INT64> rtKernel = performRankTest(RankTestColumn);
272 
273  if (rtKernel.numCols() > 1)
274  {
275  return;
276  }
277 
278  // Assert columns are not null
279  if (pPositive != NULL && pNegative != NULL)
280  {
281  // We need to check whether the existing matrix contains already a leaf which is a superset
282  if (Intersection.isExtremeRay(NullColumns))
283  {
284  // We are sure that the previous Null Tree did not contain any
285  // super sets, however the new columns may.
286 
287  // We check whether the new columns do not already contain a superset
288  if (Intersection.isExtremeRay(mNewColumns))
289  {
290  CStepMatrixColumn * pColumn = mpStepMatrix->addColumn(Intersection, pPositive, pNegative);
291 
292  // Remove all new column which are no longer extreme rays
293  std::vector< CStepMatrixColumn * >::iterator it = mNewColumns.begin();
294  std::vector< CStepMatrixColumn * >::iterator end = mNewColumns.end();
295 
296  for (; it != end; ++it)
297  {
298  if ((*it) != NULL &&
299  (*it)->getZeroSet() >= Intersection)
300  {
301  // We need to remove the column from the step matrix, which can be done by deleting it.
303  *it = NULL;
304  }
305  }
306 
307  mNewColumns.push_back(pColumn);
308  }
309  }
310 
312 
313  if (mpCallBack)
315  }
316 }
317 
318 void CBitPatternMethod::findRemoveInvalidColumns(const std::vector< CStepMatrixColumn * > & nullColumns)
319 {
320  if (mNewColumns.empty())
321  {
322  return;
323  }
324 
325  // Determine the columns which became invalid.
326  std::vector< CStepMatrixColumn * > InvalidColumns;
327 
328  std::vector< CStepMatrixColumn * >::const_iterator NullIt = nullColumns.begin();
329  std::vector< CStepMatrixColumn * >::const_iterator NullEnd = nullColumns.end();
330 
331  for (; NullIt != NullEnd; ++NullIt)
332  {
333  if (!((*NullIt)->getZeroSet()).isExtremeRay(mNewColumns))
334  {
335  InvalidColumns.push_back(*NullIt);
336  }
337  }
338 
339  mpStepMatrix->removeInvalidColumns(InvalidColumns);
340  mNewColumns.clear();
341 }
342 
344 {
345  // Calculate the kernel matrix
346  // of the reduced stoichiometry matrix to get the kernel matrix for the:
347  // Nullspace Approach to Determine the Elementary Modes of Chemical Reaction Systems (Wagner 2004)
348 
351 
352  size_t ReactionCounter = 0;
353 
354  for (; itReaction != endReaction; ++itReaction, ++ReactionCounter)
355  {
356  if ((*itReaction)->isReversible())
357  {
358  // mpReorderedReactons->push_back(*itReaction);
359  mReactionForward.push_back(std::make_pair(ReactionCounter, false));
360  }
361 
362  mpReorderedReactions->push_back(*itReaction);
363  mReactionForward.push_back(std::make_pair(ReactionCounter, true));
364  }
365 
366  const CMatrix< C_FLOAT64 > & Stoi = mpModel->getRedStoi();
367 
368  size_t NumReactions = Stoi.numCols();
369 
370  size_t NumExpandedReactions = mReactionForward.size();
371 
372  size_t NumSpecies = Stoi.numRows();
373 
374  size_t Dim = std::min(NumExpandedReactions, NumSpecies);
375 
376  if (Dim == 0)
377  {
378  return;
379  }
380 
381  mExpandedStoiTranspose.resize(NumExpandedReactions, NumSpecies);
382 
383  const C_FLOAT64 *pStoi = Stoi.array();
384  const C_FLOAT64 *pStoiEnd = pStoi + Stoi.size();
385  const C_FLOAT64 *pStoiRowEnd;
386 
387  C_INT64 *pExpandedStoiTranspose;
388  C_INT64 *pExpandedStoiTransposeColumn = mExpandedStoiTranspose.array();
389 
390  std::vector< const CReaction * >::const_iterator itReactionPivot;
391  std::vector< const CReaction * >::const_iterator endReactionPivot;
392  std::vector< std::pair< size_t, bool > >::const_iterator itReactionExpansion;
393 
394  for (; pStoi != pStoiEnd; ++pExpandedStoiTransposeColumn)
395  {
396  pStoiRowEnd = pStoi + NumReactions;
397  pExpandedStoiTranspose = pExpandedStoiTransposeColumn;
398  itReactionExpansion = mReactionForward.begin();
399 
400  for (; pStoi < pStoiRowEnd; ++pStoi, pExpandedStoiTranspose += NumSpecies, ++itReactionExpansion)
401  {
402  // TODO We should check the we have integer stoichiometry.
403  if (itReactionExpansion->second == false)
404  {
405  *pExpandedStoiTranspose = (C_INT64) - floor(*pStoi + 0.5);
406 
407  // Advance the iterators
408  ++itReactionExpansion;
409  pExpandedStoiTranspose += NumSpecies;
410  }
411 
412  *pExpandedStoiTranspose = (C_INT64) floor(*pStoi + 0.5);
413  }
414  }
415 
416  // Calculate the kernel of the matrix
417  CMatrix< C_INT64 > ExpandedStoiTranspose(mExpandedStoiTranspose);
418  CalculateKernel(ExpandedStoiTranspose, kernelInt, mReactionPivot);
419 
420  return;
421 }
422 
424 {
427 
428  CVector< size_t > Indexes;
429  size_t NumSpecies = mExpandedStoiTranspose.numCols();
430 
431  for (; it != end; ++it)
432  {
433  getUnsetBitIndexes(*it, Indexes);
434 
435  size_t NumReactions = Indexes.size();
436 
437  // Remove trivial modes, i.e., reversible reactions
438  if (NumReactions == 2 &&
439  (*mpReorderedReactions)[mReactionForward[Indexes[0]].first] ==
440  (*mpReorderedReactions)[mReactionForward[Indexes[1]].first])
441  {
442  continue;
443  }
444 
445  // Build the stoichiometry matrix reduced to the reactions participating in the current mode.
446  CMatrix< C_INT64 > A(NumReactions, NumSpecies);
447 
448  size_t * pIndex = Indexes.array();
449  size_t * pIndexEnd = pIndex + NumReactions;
450  C_INT64 * pARow = A.array();
451 
452  for (; pIndex != pIndexEnd; ++pIndex, pARow += NumSpecies)
453  {
454  memcpy(pARow, &mExpandedStoiTranspose(*pIndex, 0), NumSpecies * sizeof(C_INT64));
455  }
456 
457  // Calculate the kernel of the matrix
458  CMatrix< C_INT64 > ExpandedStoiTranspose(A);
459  CMatrix< C_INT64 > Kernel;
460  CVector< size_t > Pivot;
461  CalculateKernel(ExpandedStoiTranspose, Kernel, Pivot);
462 
463  size_t NumCols = Kernel.numCols();
464 
465  // Now we create the flux mode as we have the multiplier and reaction indexes.
466  // We need to invert the sign of the multiplier for reactions which are not forward.
467  // A flux mode is reversible if all reactions are reversible;
468 
469  C_INT64 * pColumn = Kernel.array();
470  C_INT64 * pColumnEnd = pColumn + NumCols;
471 
472  for (; pColumn != pColumnEnd; ++pColumn)
473  {
474  std::map< size_t, C_FLOAT64 > Reactions;
475  bool Reversible = true;
476 
477  pIndex = Indexes.array();
478  C_INT64 * pFluxMultiplier = pColumn;
479 
480  for (; pIndex != pIndexEnd; ++pIndex, pFluxMultiplier += NumCols)
481  {
482  if (*pFluxMultiplier < 0)
483  {
484  break;
485  }
486  else if (*pFluxMultiplier < 0)
487  {
488  continue;
489  }
490 
491  std::pair< size_t, bool > & ReactionForward = mReactionForward[*pIndex];
492 
493  Reactions[ReactionForward.first] =
494  (C_FLOAT64)((ReactionForward.second == true) ? *pFluxMultiplier : -*pFluxMultiplier);
495 
496  if (!(*mpReorderedReactions)[ReactionForward.first]->isReversible())
497  {
498  Reversible = false;
499  }
500  }
501 
502  if (pIndex != pIndexEnd)
503  {
504  continue;
505  }
506 
507  addMode(CFluxMode(Reactions, Reversible));
508  }
509  }
510 }
511 
512 #ifdef XXXX
513 void CBitPatternMethod::convertToIntegers(CMatrix< C_FLOAT64 > & values)
514 {
515  bool Problems = false;
516 
517  static const C_FLOAT64 limit = 10.0 / std::numeric_limits< C_INT32 >::max();
518 
519  size_t Size = values.size();
520  size_t Columns = values.numCols();
521 
522  C_FLOAT64 * pColumn = values.array();
523  C_FLOAT64 * pColumnEnd = pColumn + Columns;
524  C_FLOAT64 * pValue = values.array();
525  C_FLOAT64 * pValueEnd = pColumn + Size;
526 
527  for (; pColumn < pColumnEnd; ++pColumn)
528  {
529  size_t Multiplier = 1;
530  size_t m00, m01, m10, m11;
531  size_t maxden = 10000000;
532  C_INT32 GCD1, GCD2;
533  size_t ai;
534 
535  C_FLOAT64 x;
536 
537  for (pValue = pColumn; pValue < pValueEnd; pValue += Columns)
538  {
539  x = fabs(*pValue);
540 
541  if (x < 100 * std::numeric_limits< C_FLOAT64 >::epsilon())
542  {
543  continue;
544  }
545 
546  /*
547  * Find rational approximation to given real number
548  * David Eppstein / UC Irvine / 8 Aug 1993
549  *
550  * With corrections from:
551  * Arno Formella, May 2008
552  * Stefan Hoops, Sept 2009
553  *
554  * Based on the theory of continued fractions
555  * if x = a1 + 1/(a2 + 1/(a3 + 1/(a4 + ...)))
556  * then best approximation is found by truncating this series
557  * (with some adjustments in the last term).
558  *
559  * Note the fraction can be recovered as the first column of the matrix
560  * (a1 1 ) (a2 1 ) (a3 1 ) ...
561  * (1 0 ) (1 0 ) (1 0)
562  * Instead of keeping the sequence of continued fraction terms,
563  * we just keep the last partial product of these matrices.
564  */
565 
566  /* initialize matrix */
567  m00 = m11 = 1;
568  m01 = m10 = 0;
569 
570  /* loop finding terms until denom gets too big */
571  while (m10 * (ai = (size_t) x) + m11 <= maxden)
572  {
573  size_t t;
574  t = m00 * ai + m01;
575  m01 = m00;
576  m00 = t;
577 
578  t = m10 * ai + m11;
579  m11 = m10;
580  m10 = t;
581 
582  if (fabs(x - (C_FLOAT64) ai) < limit)
583  break; // SH: We reached the numerical precision of the machine;
584 
585  x = 1 / (x - (C_FLOAT64) ai);
586  }
587 
588  if ((C_FLOAT64) m10 * (C_FLOAT64) ai + (C_FLOAT64) m11 > (C_FLOAT64) maxden)
589  {
590  Problems = true;
591  }
592 
593  if (fabs(fabs(*pValue) - ((C_FLOAT64) m00) / ((C_FLOAT64) m10)) > 100.0 * std::numeric_limits< C_FLOAT64 >::epsilon())
594  {
595  ai = (maxden - m11) / m10;
596  m00 = m00 * ai + m01;
597  m10 = m10 * ai + m11;
598  }
599 
600  // Find the greatest common divisor (GCD) of the multiplier and the current denominator.
601  // Euclidean algorithm
602  GCD1 = m10;
603  GCD2 = Multiplier;
604 
605  GCD(GCD1, GCD2);
606 
607  // Calculate the least common multiplier: LCM = v1 * v2 / GCD(v1, v2)
608  Multiplier *= m10 / GCD1;
609  }
610 
611  for (pValue = pColumn; pValue < pValueEnd; pValue += Columns)
612  {
613  *pValue *= Multiplier;
614  }
615  }
616 }
617 #endif //
618 
620  CVector<size_t> & indexes) const
621 {
622  mpStepMatrix->getAllUnsetBitIndexes(pColumn, indexes);
623 
624  // Apply the QR pivot
625  size_t * pIndex = indexes.array();
626  size_t * pIndexEnd = pIndex + indexes.size();
627 
628  for (; pIndex != pIndexEnd; ++pIndex)
629  *pIndex = mReactionPivot[*pIndex];
630 
631  //DebugFile << "@CBPM: " << indexes << std::endl;
632 }
633 
635  CVector< size_t > & indexes) const
636 {
637  mpStepMatrix->getUnsetBitIndexes(pColumn, indexes);
638 
639  // Apply the QR pivot
640  size_t * pIndex = indexes.array();
641  size_t * pIndexEnd = pIndex + indexes.size();
642 
643  for (; pIndex != pIndexEnd; ++pIndex)
644  {
645  *pIndex = mReactionPivot[*pIndex];
646  }
647 }
648 
649 // private
651 {
652  std::vector< CFluxMode >::iterator itMode = mpFluxModes->begin();
653  std::vector< CFluxMode >::iterator endMode = mpFluxModes->end();
654 
655  for (; itMode != endMode; ++itMode)
656  {
657  if (itMode->isReversed(mode))
658  {
659  return;
660  }
661  }
662 
663  mpFluxModes->push_back(mode);
664  return;
665 }
666 
667 // static
669  CMatrix< C_INT64 > & kernel,
670  CVector< size_t > & rowPivot)
671 {
672  // Gaussian elimination
673  size_t NumRows = matrix.numRows();
674  size_t NumCols = matrix.numCols();
675 
676  assert(NumRows > 0);
677  assert(NumCols > 0);
678 
679  // Initialize row pivots
680  rowPivot.resize(NumRows);
681  size_t * pPivot = rowPivot.array();
682 
683  for (size_t i = 0; i < NumRows; i++, pPivot++)
684  {
685  *pPivot = i;
686  }
687 
688  CVector< size_t > RowPivot(rowPivot);
689 
690  CVector< C_INT64 > Identity(NumRows);
691  Identity = 1;
692 
693  C_INT64 * pColumn = matrix.array();
694  C_INT64 * pColumnEnd = pColumn + NumCols;
695 
696  C_INT64 * pActiveRow;
697  C_INT64 * pActiveRowStart = pColumn;
698  C_INT64 * pActiveRowEnd = pColumnEnd;
699 
700  C_INT64 * pRow;
701  C_INT64 * pRowEnd = pColumn + matrix.size();
702 
703  C_INT64 * pCurrent;
704  C_INT64 * pIdentity;
705 
706  CVector< C_INT64 > SwapTmp(NumCols);
707 
708  size_t CurrentRowIndex = 0;
709  size_t CurrentColumnIndex = 0;
710  size_t NonZeroIndex = 0;
711 
712  CVector< bool > IgnoredColumns(NumCols);
713  IgnoredColumns = false;
714  bool * pIgnoredColumn = IgnoredColumns.array();
715  size_t IgnoredColumnCount = 0;
716 
717  // For each column
718  for (; pColumn < pColumnEnd; ++pColumn, ++CurrentColumnIndex, ++pIgnoredColumn)
719  {
720  assert(CurrentColumnIndex == CurrentRowIndex + IgnoredColumnCount);
721 
722  // Find non zero entry in current column.
723  pRow = pActiveRowStart + CurrentColumnIndex;
724  NonZeroIndex = CurrentRowIndex;
725 
726  for (; pRow < pRowEnd; pRow += NumCols, ++NonZeroIndex)
727  {
728  if (*pRow != 0)
729  break;
730  }
731 
732  if (NonZeroIndex >= NumRows)
733  {
734  *pIgnoredColumn = true;
735  IgnoredColumnCount++;
736 
737  continue;
738  }
739 
740  if (NonZeroIndex != CurrentRowIndex)
741  {
742  // Swap rows
743  memcpy(SwapTmp.array(), matrix[CurrentRowIndex], NumCols * sizeof(C_INT64));
744  memcpy(matrix[CurrentRowIndex], matrix[NonZeroIndex], NumCols * sizeof(C_INT64));
745  memcpy(matrix[NonZeroIndex], SwapTmp.array(), NumCols * sizeof(C_INT64));
746 
747  // Record pivot
748  size_t tmp = RowPivot[CurrentRowIndex];
749  RowPivot[CurrentRowIndex] = RowPivot[NonZeroIndex];
750  RowPivot[NonZeroIndex] = tmp;
751 
752  C_INT64 TMP = Identity[CurrentRowIndex];
753  Identity[CurrentRowIndex] = Identity[NonZeroIndex];
754  Identity[NonZeroIndex] = TMP;
755  }
756 
757  if (*(pActiveRowStart + CurrentColumnIndex) < 0)
758  {
759  for (pRow = pActiveRowStart; pRow < pActiveRowEnd; ++pRow)
760  {
761  *pRow *= -1;
762  }
763 
764  Identity[CurrentRowIndex] *= -1;
765  }
766 
767  // For each row
768  pRow = pActiveRowStart + NumCols;
769  pIdentity = Identity.array() + CurrentRowIndex + 1;
770 
771  C_INT64 ActiveRowValue = *(pActiveRowStart + CurrentColumnIndex);
772  *(pActiveRowStart + CurrentColumnIndex) = Identity[CurrentRowIndex];
773 
774  for (; pRow < pRowEnd; pRow += NumCols, ++pIdentity)
775  {
776  C_INT64 RowValue = *(pRow + CurrentColumnIndex);
777 
778  if (RowValue == 0)
779  continue;
780 
781  *(pRow + CurrentColumnIndex) = 0;
782 
783  // compute GCD(*pActiveRowStart, *pRow)
784  C_INT64 GCD1 = abs64(ActiveRowValue);
785  C_INT64 GCD2 = abs64(RowValue);
786 
787  GCD(GCD1, GCD2);
788 
789  C_INT64 alpha = ActiveRowValue / GCD1;
790  C_INT64 beta = RowValue / GCD1;
791 
792  // update rest of row
793  pActiveRow = pActiveRowStart;
794  pCurrent = pRow;
795  *pIdentity *= alpha;
796 
797  GCD1 = abs64(*pIdentity);
798 
799  for (; pActiveRow < pActiveRowEnd; ++pActiveRow, ++pCurrent)
800  {
801  // Assert that we do not have a numerical overflow.
802  assert(fabs(((C_FLOAT64) alpha) * ((C_FLOAT64) * pCurrent) - ((C_FLOAT64) beta) * ((C_FLOAT64) * pActiveRow)) < std::numeric_limits< C_INT64 >::max());
803 
804  *pCurrent = alpha * *pCurrent - beta * *pActiveRow;
805 
806  // We check that the row values do not have any common divisor.
807  if (GCD1 > 1 &&
808  (GCD2 = abs64(*pCurrent)) > 0)
809  {
810  GCD(GCD1, GCD2);
811  }
812  }
813 
814  if (GCD1 > 1)
815  {
816  *pIdentity /= GCD1;
817 
818  pActiveRow = pActiveRowStart;
819  pCurrent = pRow;
820 
821  for (; pActiveRow < pActiveRowEnd; ++pActiveRow, ++pCurrent)
822  {
823  *pCurrent /= GCD1;
824  }
825  }
826  }
827 
828  pActiveRowStart += NumCols;
829  pActiveRowEnd += NumCols;
830  CurrentRowIndex++;
831 
832  // DebugFile << matrix << std::endl;
833  // DebugFile << Identity << std::endl;
834  }
835 
836  assert(CurrentColumnIndex == CurrentRowIndex + IgnoredColumnCount);
837  assert(CurrentColumnIndex == NumCols);
838 
839  // std::cout << matrix << std::endl;
840  // std::cout << IgnoredColumns << std::endl;
841  // std::cout << Identity << std::endl;
842  // std::cout << RowPivot << std::endl << std::endl;
843 
844  size_t KernelRows = NumRows;
845  size_t KernelCols = NumRows + IgnoredColumnCount - NumCols;
846 
847  assert(KernelCols > 0);
848 
849  kernel.resize(KernelRows, KernelCols);
850 
851  CMatrix< C_INT64 > Kernel(KernelRows, KernelCols);
852  Kernel = 0;
853 
854  C_INT64 * pKernelInt = Kernel.array();
855  C_INT64 * pKernelIntEnd = pKernelInt + Kernel.size();
856 
857  pActiveRowStart = matrix[CurrentRowIndex];
858  pActiveRowEnd = matrix[NumRows];
859 
860  // Null space matrix identity part
861  pIdentity = Identity.array() + CurrentRowIndex;
862  C_INT64 * pKernelColumn = Kernel.array();
863 
864  for (; pActiveRowStart < pActiveRowEnd; pActiveRowStart += NumCols, ++pKernelColumn, ++pIdentity)
865  {
866  pKernelInt = pKernelColumn;
867  pIgnoredColumn = IgnoredColumns.array();
868 
869  pRow = pActiveRowStart;
870  pRowEnd = pRow + NumCols;
871 
872  if (*pIdentity < 0)
873  {
874  *pIdentity *= -1;
875 
876  for (; pRow < pRowEnd; ++pRow, ++pIgnoredColumn)
877  {
878  if (*pIgnoredColumn)
879  {
880  continue;
881  }
882 
883  *pKernelInt = - *pRow;
884  pKernelInt += KernelCols;
885  }
886  }
887  else
888  {
889  for (; pRow < pRowEnd; ++pRow, ++pIgnoredColumn)
890  {
891  if (*pIgnoredColumn)
892  {
893  continue;
894  }
895 
896  *pKernelInt = *pRow;
897  pKernelInt += KernelCols;
898  }
899  }
900  }
901 
902  pIdentity = Identity.array() + CurrentRowIndex;
903  pKernelInt = Kernel[CurrentRowIndex];
904 
905  for (; pKernelInt < pKernelIntEnd; pKernelInt += KernelCols + 1, ++pIdentity)
906  {
907  *pKernelInt = *pIdentity;
908  }
909 
910  // Undo the reordering introduced by Gaussian elimination to the kernel matrix.
911  pPivot = RowPivot.array();
912  pRow = Kernel.array();
913  pRowEnd = pRow + Kernel.size();
914 
915  for (; pRow < pRowEnd; ++pPivot, pRow += KernelCols)
916  {
917  memcpy(kernel[*pPivot], pRow, KernelCols * sizeof(C_INT64));
918  }
919 
920  return true;
921 }
922 
924 {
925  CVector<size_t> Indexes;
926 
927  getAllUnsetBitIndexes(pIntersectColumn, Indexes);
928 
929  size_t NumReactions = Indexes.size();
930  size_t NumSpecies = mExpandedStoiTranspose.numCols();
931 
932  // Build the stoichiometry matrix reduced to the reactions participating in the current mode.
933  CMatrix< C_INT64 > A(NumReactions, NumSpecies);
934 
935  size_t * pIndex = Indexes.array();
936  size_t * pIndexEnd = pIndex + NumReactions;
937  C_INT64 * pARow = A.array();
938 
939  for (; pIndex != pIndexEnd; ++pIndex, pARow += NumSpecies)
940  {
941  memcpy(pARow, &mExpandedStoiTranspose(*pIndex, 0), NumSpecies * sizeof(C_INT64));
942  }
943 
944  // Calculate the kernel of the matrix
945  CMatrix< C_INT64 > ExpandedStoiTranspose(A);
946  CMatrix< C_INT64 > Kernel;
947  CVector< size_t > Pivot;
948 
949  //DebugFile << ExpandedStoiTranspose << std::endl;
950  CalculateKernel(ExpandedStoiTranspose, Kernel, Pivot);
951 
952  return Kernel;
953 }
954 
956  std::stack<CStepMatrixColumn *> kernelStack)
957 {
958  CStepMatrixColumn * KernelCol;
959  unsigned int col, row;
960 
961  while (!kernelStack.empty())
962  {
963  KernelCol = kernelStack.top();
964 #ifdef COPASI_DEBUG_TRACE
965  DebugFile << "Current Column: " << *KernelCol << std::endl;
966 #endif // COPASI_DEBUG_TRACE
967  CMatrix<C_INT64> RankKernel = performRankTest(KernelCol);
968 #ifdef COPASI_DEBUG_TRACE
969  DebugFile << "New Rank Kernel " << RankKernel << std::endl;
970 #endif // COPASI_DEBUG_TRACE
971 
972  if (RankKernel.numCols() == 1)
973  {
974  fluxModeMatrix->add(KernelCol);
975  }
976  else if (RankKernel.numCols() > 1)
977  {
978  for (col = 0; col < RankKernel.numCols(); col++)
979  {
980  KernelCol = new CStepMatrixColumn(RankKernel.numRows());
981 
982  for (row = 0; row < RankKernel.numRows(); row++)
983  {
984  KernelCol->push_front(RankKernel(row, col));
985  }
986 
987  kernelStack.push(KernelCol);
988  }
989  }
990 
991  kernelStack.pop();
992  }
993 }
994 
995 /*
996 Future Implementation:
997 After building left kernel of the stoichiometric matrix 'SP', parse through each vector column 'vi' in SP.
998  SP contains "pre-flux modes" which must be checked based on rank test.
999 Build the stoichiometric transpose submatrix NTvi based on involved reactions in vi.
1000 Calculate the left kernel of NTvi, 'K'.
1001 If #cols in K = 1:
1002  Move vi from SP into new step matrix 'SF'.
1003 Else if #cols in K > 1:
1004  Add each unique nontrivial kernel vector into SP.
1005 
1006 The end result is that SF will contain only flux modes, but not all possible flux modes.
1007 Through this method we never have to remove any elements from the matrix.
1008 
1009 After this we then build the initial step matrix from SF.
1010  */
void removeInvalidColumns(std::vector< CStepMatrixColumn * > &invalidColumns)
void buildKernelMatrix(CMatrix< C_INT64 > &kernel)
void combine(const CStepMatrixColumn *pPositive, const CStepMatrixColumn *pNegative, const std::vector< CStepMatrixColumn * > NullColumns)
static void GCD(C_INT64 &m, C_INT64 &n)
#define pdelete(p)
Definition: copasi.h:215
CCopasiProblem * getProblem()
void push_front(const C_INT64 &value)
void getUnsetBitIndexes(const CStepMatrixColumn *pColumn, CVector< size_t > &indexes) const
void buildFluxModeMatrix(CStepMatrix *fluxModeMatrix, std::stack< CStepMatrixColumn * > kernelStack)
size_t getFirstUnconvertedRow() const
const CZeroSet & getZeroSet() const
iterator begin()
virtual size_t numRows() const
Definition: CMatrix.h:138
void getAllUnsetBitIndexes(const CStepMatrixColumn *pColumn, CVector< size_t > &indexes) const
CStepMatrixColumn *const * const_iterator
Definition: CStepMatrix.h:29
std::vector< CFluxMode > * mpFluxModes
Definition: CEFMMethod.h:91
void add(CStepMatrixColumn *pColumn)
Definition: CStepMatrix.h:50
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
#define C_INT64
Definition: copasi.h:88
unsigned C_INT32 mProgressCounterMax
#define C_INT32
Definition: copasi.h:90
void findRemoveInvalidColumns(const std::vector< CStepMatrixColumn * > &nullColumns)
virtual bool calculate()
const CMatrix< C_FLOAT64 > & getRedStoi() const
Definition: CModel.cpp:1154
virtual bool progressItem(const size_t &handle)
virtual bool initialize()
Definition: CEFMMethod.cpp:109
const_iterator begin() const
bool splitColumns(std::vector< CStepMatrixColumn * > &PositiveColumns, std::vector< CStepMatrixColumn * > &NegativeColumns, std::vector< CStepMatrixColumn * > &NullColumns)
CTSSATask * pTask
unsigned C_INT32 mProgressCounter2Max
virtual bool proceed()
iterator end()
size_t getNumUnconvertedRows() const
const size_t & getNumberOfSetBits() const
Definition: CZeroSet.h:85
std::vector< CType * >::const_iterator const_iterator
Definition: CCopasiVector.h:57
void getAllUnsetBitIndexes(const CStepMatrixColumn *pColumn, CVector< size_t > &indexes) const
const_iterator end() const
void removeColumn(CStepMatrixColumn *pColumn)
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
void addMode(const CFluxMode &mode)
std::vector< const CReaction * > * mpReorderedReactions
Definition: CEFMMethod.h:96
void convertRow()
#define abs64
Definition: copasi.h:94
unsigned C_INT32 mProgressCounter
virtual bool finishItem(const size_t &handle)
CStepMatrix * mpStepMatrix
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
size_t size() const
Definition: CVector.h:100
CBitPatternMethod(const CCopasiContainer *pParent=NULL)
bool isExtremeRay(const std::vector< CStepMatrixColumn * > &columns) const
Definition: CZeroSet.cpp:97
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
CMatrix< C_INT64 > mExpandedStoiTranspose
CMatrix< C_INT64 > performRankTest(CStepMatrixColumn *pIntersectColumn)
virtual size_t size() const
Definition: CMatrix.h:132
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
std::vector< std::pair< size_t, bool > > mReactionForward
void compact()
static CZeroSet intersection(const CZeroSet &set1, const CZeroSet &set2)
Definition: CZeroSet.h:137
virtual bool initialize()
CVector< size_t > mReactionPivot
virtual size_t numCols() const
Definition: CMatrix.h:144
CStepMatrixColumn * addColumn(const CZeroSet &set, const CStepMatrixColumn *pPositive, const CStepMatrixColumn *pNegative)
void getUnsetBitIndexes(const CStepMatrixColumn *pColumn, CVector< size_t > &indexes) const
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
CProcessReport * mpCallBack
virtual CType * array()
Definition: CMatrix.h:337
CModel * getModel() const
#define min(a, b)
Definition: f2c.h:175
static bool CalculateKernel(CMatrix< C_INT64 > &matrix, CMatrix< C_INT64 > &kernel, CVector< size_t > &rowPivot)
CCopasiContainer * getObjectParent() const
std::vector< CStepMatrixColumn * > mNewColumns
unsigned C_INT32 mProgressCounter2
#define max(a, b)
Definition: f2c.h:176