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