31 mProgressCounterMax(0),
34 mProgressCounter2Max(0),
35 mhProgressCounter2(0),
38 mExpandedStoiTranspose(0, 0),
42 mContinueCombination(true)
52 mProgressCounterMax(0),
55 mProgressCounter2Max(0),
56 mhProgressCounter2(0),
59 mExpandedStoiTranspose(0, 0),
63 mContinueCombination(true)
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),
84 mContinueCombination(src.mContinueCombination)
112 if (pTask == NULL)
return false;
116 if (
mpModel == NULL)
return false;
124 std::stack< CStepMatrixColumn * > KernelStack;
127 for (
unsigned int col = 0; col < KernelMatrix.
numCols(); col++)
131 for (
unsigned int row = 0; row < KernelMatrix.
numRows(); row++)
133 KernelCol->
push_front(KernelMatrix(row, col));
136 KernelStack.push(KernelCol);
167 bool Continue =
true;
182 std::vector< CStepMatrixColumn * > PositiveColumns;
183 std::vector< CStepMatrixColumn * > NegativeColumns;
184 std::vector< CStepMatrixColumn * > NullColumns;
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);
247 const std::vector< CStepMatrixColumn * > NullColumns)
279 if (pPositive != NULL && pNegative != NULL)
293 std::vector< CStepMatrixColumn * >::iterator it =
mNewColumns.begin();
294 std::vector< CStepMatrixColumn * >::iterator end =
mNewColumns.end();
296 for (; it != end; ++it)
299 (*it)->getZeroSet() >= Intersection)
326 std::vector< CStepMatrixColumn * > InvalidColumns;
328 std::vector< CStepMatrixColumn * >::const_iterator NullIt = nullColumns.begin();
329 std::vector< CStepMatrixColumn * >::const_iterator NullEnd = nullColumns.end();
331 for (; NullIt != NullEnd; ++NullIt)
333 if (!((*NullIt)->getZeroSet()).isExtremeRay(
mNewColumns))
335 InvalidColumns.push_back(*NullIt);
352 size_t ReactionCounter = 0;
354 for (; itReaction != endReaction; ++itReaction, ++ReactionCounter)
356 if ((*itReaction)->isReversible())
368 size_t NumReactions = Stoi.
numCols();
372 size_t NumSpecies = Stoi.
numRows();
374 size_t Dim =
std::min(NumExpandedReactions, NumSpecies);
387 C_INT64 *pExpandedStoiTranspose;
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;
394 for (; pStoi != pStoiEnd; ++pExpandedStoiTransposeColumn)
396 pStoiRowEnd = pStoi + NumReactions;
397 pExpandedStoiTranspose = pExpandedStoiTransposeColumn;
400 for (; pStoi < pStoiRowEnd; ++pStoi, pExpandedStoiTranspose += NumSpecies, ++itReactionExpansion)
403 if (itReactionExpansion->second ==
false)
405 *pExpandedStoiTranspose = (
C_INT64) - floor(*pStoi + 0.5);
408 ++itReactionExpansion;
409 pExpandedStoiTranspose += NumSpecies;
412 *pExpandedStoiTranspose = (
C_INT64) floor(*pStoi + 0.5);
431 for (; it != end; ++it)
435 size_t NumReactions = Indexes.
size();
438 if (NumReactions == 2 &&
448 size_t * pIndex = Indexes.
array();
449 size_t * pIndexEnd = pIndex + NumReactions;
452 for (; pIndex != pIndexEnd; ++pIndex, pARow += NumSpecies)
463 size_t NumCols = Kernel.
numCols();
470 C_INT64 * pColumnEnd = pColumn + NumCols;
472 for (; pColumn != pColumnEnd; ++pColumn)
474 std::map< size_t, C_FLOAT64 > Reactions;
475 bool Reversible =
true;
477 pIndex = Indexes.
array();
478 C_INT64 * pFluxMultiplier = pColumn;
480 for (; pIndex != pIndexEnd; ++pIndex, pFluxMultiplier += NumCols)
482 if (*pFluxMultiplier < 0)
486 else if (*pFluxMultiplier < 0)
493 Reactions[ReactionForward.first] =
494 (
C_FLOAT64)((ReactionForward.second ==
true) ? *pFluxMultiplier : -*pFluxMultiplier);
502 if (pIndex != pIndexEnd)
515 bool Problems =
false;
519 size_t Size = values.
size();
520 size_t Columns = values.
numCols();
523 C_FLOAT64 * pColumnEnd = pColumn + Columns;
527 for (; pColumn < pColumnEnd; ++pColumn)
529 size_t Multiplier = 1;
530 size_t m00, m01, m10, m11;
531 size_t maxden = 10000000;
537 for (pValue = pColumn; pValue < pValueEnd; pValue += Columns)
541 if (x < 100 * std::numeric_limits< C_FLOAT64 >::epsilon())
571 while (m10 * (ai = (
size_t) x) + m11 <= maxden)
593 if (fabs(fabs(*pValue) - ((
C_FLOAT64) m00) / ((
C_FLOAT64) m10)) > 100.0 * std::numeric_limits< C_FLOAT64 >::epsilon())
595 ai = (maxden - m11) / m10;
596 m00 = m00 * ai + m01;
597 m10 = m10 * ai + m11;
608 Multiplier *= m10 / GCD1;
611 for (pValue = pColumn; pValue < pValueEnd; pValue += Columns)
613 *pValue *= Multiplier;
625 size_t * pIndex = indexes.
array();
626 size_t * pIndexEnd = pIndex + indexes.
size();
628 for (; pIndex != pIndexEnd; ++pIndex)
640 size_t * pIndex = indexes.
array();
641 size_t * pIndexEnd = pIndex + indexes.
size();
643 for (; pIndex != pIndexEnd; ++pIndex)
652 std::vector< CFluxMode >::iterator itMode =
mpFluxModes->begin();
653 std::vector< CFluxMode >::iterator endMode =
mpFluxModes->end();
655 for (; itMode != endMode; ++itMode)
657 if (itMode->isReversed(mode))
673 size_t NumRows = matrix.
numRows();
674 size_t NumCols = matrix.
numCols();
681 size_t * pPivot = rowPivot.
array();
683 for (
size_t i = 0; i < NumRows; i++, pPivot++)
694 C_INT64 * pColumnEnd = pColumn + NumCols;
697 C_INT64 * pActiveRowStart = pColumn;
698 C_INT64 * pActiveRowEnd = pColumnEnd;
708 size_t CurrentRowIndex = 0;
709 size_t CurrentColumnIndex = 0;
710 size_t NonZeroIndex = 0;
713 IgnoredColumns =
false;
714 bool * pIgnoredColumn = IgnoredColumns.
array();
715 size_t IgnoredColumnCount = 0;
718 for (; pColumn < pColumnEnd; ++pColumn, ++CurrentColumnIndex, ++pIgnoredColumn)
720 assert(CurrentColumnIndex == CurrentRowIndex + IgnoredColumnCount);
723 pRow = pActiveRowStart + CurrentColumnIndex;
724 NonZeroIndex = CurrentRowIndex;
726 for (; pRow < pRowEnd; pRow += NumCols, ++NonZeroIndex)
732 if (NonZeroIndex >= NumRows)
734 *pIgnoredColumn =
true;
735 IgnoredColumnCount++;
740 if (NonZeroIndex != CurrentRowIndex)
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));
748 size_t tmp = RowPivot[CurrentRowIndex];
749 RowPivot[CurrentRowIndex] = RowPivot[NonZeroIndex];
750 RowPivot[NonZeroIndex] = tmp;
752 C_INT64 TMP = Identity[CurrentRowIndex];
753 Identity[CurrentRowIndex] = Identity[NonZeroIndex];
754 Identity[NonZeroIndex] = TMP;
757 if (*(pActiveRowStart + CurrentColumnIndex) < 0)
759 for (pRow = pActiveRowStart; pRow < pActiveRowEnd; ++pRow)
764 Identity[CurrentRowIndex] *= -1;
768 pRow = pActiveRowStart + NumCols;
769 pIdentity = Identity.
array() + CurrentRowIndex + 1;
771 C_INT64 ActiveRowValue = *(pActiveRowStart + CurrentColumnIndex);
772 *(pActiveRowStart + CurrentColumnIndex) = Identity[CurrentRowIndex];
774 for (; pRow < pRowEnd; pRow += NumCols, ++pIdentity)
776 C_INT64 RowValue = *(pRow + CurrentColumnIndex);
781 *(pRow + CurrentColumnIndex) = 0;
789 C_INT64 alpha = ActiveRowValue / GCD1;
790 C_INT64 beta = RowValue / GCD1;
793 pActiveRow = pActiveRowStart;
797 GCD1 =
abs64(*pIdentity);
799 for (; pActiveRow < pActiveRowEnd; ++pActiveRow, ++pCurrent)
804 *pCurrent = alpha * *pCurrent - beta * *pActiveRow;
808 (GCD2 =
abs64(*pCurrent)) > 0)
818 pActiveRow = pActiveRowStart;
821 for (; pActiveRow < pActiveRowEnd; ++pActiveRow, ++pCurrent)
828 pActiveRowStart += NumCols;
829 pActiveRowEnd += NumCols;
836 assert(CurrentColumnIndex == CurrentRowIndex + IgnoredColumnCount);
837 assert(CurrentColumnIndex == NumCols);
844 size_t KernelRows = NumRows;
845 size_t KernelCols = NumRows + IgnoredColumnCount - NumCols;
847 assert(KernelCols > 0);
849 kernel.
resize(KernelRows, KernelCols);
855 C_INT64 * pKernelIntEnd = pKernelInt + Kernel.
size();
857 pActiveRowStart = matrix[CurrentRowIndex];
858 pActiveRowEnd = matrix[NumRows];
861 pIdentity = Identity.
array() + CurrentRowIndex;
864 for (; pActiveRowStart < pActiveRowEnd; pActiveRowStart += NumCols, ++pKernelColumn, ++pIdentity)
866 pKernelInt = pKernelColumn;
867 pIgnoredColumn = IgnoredColumns.
array();
869 pRow = pActiveRowStart;
870 pRowEnd = pRow + NumCols;
876 for (; pRow < pRowEnd; ++pRow, ++pIgnoredColumn)
883 *pKernelInt = - *pRow;
884 pKernelInt += KernelCols;
889 for (; pRow < pRowEnd; ++pRow, ++pIgnoredColumn)
897 pKernelInt += KernelCols;
902 pIdentity = Identity.
array() + CurrentRowIndex;
903 pKernelInt = Kernel[CurrentRowIndex];
905 for (; pKernelInt < pKernelIntEnd; pKernelInt += KernelCols + 1, ++pIdentity)
907 *pKernelInt = *pIdentity;
911 pPivot = RowPivot.
array();
912 pRow = Kernel.
array();
913 pRowEnd = pRow + Kernel.
size();
915 for (; pRow < pRowEnd; ++pPivot, pRow += KernelCols)
917 memcpy(kernel[*pPivot], pRow, KernelCols *
sizeof(
C_INT64));
929 size_t NumReactions = Indexes.
size();
935 size_t * pIndex = Indexes.
array();
936 size_t * pIndexEnd = pIndex + NumReactions;
939 for (; pIndex != pIndexEnd; ++pIndex, pARow += NumSpecies)
956 std::stack<CStepMatrixColumn *> kernelStack)
959 unsigned int col, row;
961 while (!kernelStack.empty())
963 KernelCol = kernelStack.top();
964 #ifdef COPASI_DEBUG_TRACE
965 DebugFile <<
"Current Column: " << *KernelCol << std::endl;
966 #endif // COPASI_DEBUG_TRACE
968 #ifdef COPASI_DEBUG_TRACE
969 DebugFile <<
"New Rank Kernel " << RankKernel << std::endl;
970 #endif // COPASI_DEBUG_TRACE
974 fluxModeMatrix->
add(KernelCol);
976 else if (RankKernel.
numCols() > 1)
978 for (col = 0; col < RankKernel.
numCols(); col++)
982 for (row = 0; row < RankKernel.
numRows(); row++)
987 kernelStack.push(KernelCol);
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)
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
virtual size_t numRows() const
void getAllUnsetBitIndexes(const CStepMatrixColumn *pColumn, CVector< size_t > &indexes) const
CStepMatrixColumn *const * const_iterator
std::vector< CFluxMode > * mpFluxModes
void add(CStepMatrixColumn *pColumn)
void resize(size_t size, const bool ©=false)
unsigned C_INT32 mProgressCounterMax
void findRemoveInvalidColumns(const std::vector< CStepMatrixColumn * > &nullColumns)
const CMatrix< C_FLOAT64 > & getRedStoi() const
virtual bool progressItem(const size_t &handle)
virtual bool initialize()
const_iterator begin() const
bool splitColumns(std::vector< CStepMatrixColumn * > &PositiveColumns, std::vector< CStepMatrixColumn * > &NegativeColumns, std::vector< CStepMatrixColumn * > &NullColumns)
unsigned C_INT32 mProgressCounter2Max
size_t getNumUnconvertedRows() const
const size_t & getNumberOfSetBits() const
std::vector< CType * >::const_iterator const_iterator
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
unsigned C_INT32 mProgressCounter
virtual bool finishItem(const size_t &handle)
CStepMatrix * mpStepMatrix
virtual void resize(size_t rows, size_t cols, const bool ©=false)
CBitPatternMethod(const CCopasiContainer *pParent=NULL)
bool isExtremeRay(const std::vector< CStepMatrixColumn * > &columns) const
size_t mhProgressCounter2
CMatrix< C_INT64 > mExpandedStoiTranspose
CMatrix< C_INT64 > performRankTest(CStepMatrixColumn *pIntersectColumn)
virtual size_t size() const
CCopasiVectorNS< CReaction > & getReactions()
std::vector< std::pair< size_t, bool > > mReactionForward
static CZeroSet intersection(const CZeroSet &set1, const CZeroSet &set2)
virtual bool initialize()
CVector< size_t > mReactionPivot
virtual size_t numCols() const
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
bool mContinueCombination
CModel * getModel() const
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