COPASI API  4.16.103
CStepMatrixColumn.cpp
Go to the documentation of this file.
1 // Begin CVS Header
2 // $Source: /Volumes/Home/Users/shoops/cvs/copasi_dev/copasi/elementaryFluxModes/CStepMatrixColumn.cpp,v $
3 // $Revision: 1.12 $
4 // $Name: $
5 // $Author: ssahle $
6 // $Date: 2012/04/22 14:51:18 $
7 // End CVS Header
8 
9 // Copyright (C) 2012 - 2010 by Pedro Mendes, Virginia Tech Intellectual
10 // Properties, Inc., University of Heidelberg, and The University
11 // of Manchester.
12 // All rights reserved.
13 
14 // Copyright (C) 2008 by Pedro Mendes, Virginia Tech Intellectual
15 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
16 // and The University of Manchester.
17 // All rights reserved.
18 
19 #include <stdlib.h>
20 #include <cmath>
21 #include <limits>
22 
23 #include "copasi.h"
24 
25 #include "CStepMatrixColumn.h"
26 #include "CBitPatternTreeMethod.h"
27 
29  mZeroSet(size),
30  mReaction(),
31  mIterator(NULL)
32 {}
33 
35  CStepMatrixColumn const * pPositive,
36  CStepMatrixColumn const * pNegative):
37  mZeroSet(set),
38  mReaction(),
39  mIterator(NULL)
40 {
41  C_INT64 PosMult = -pNegative->getMultiplier();
42  C_INT64 NegMult = pPositive->getMultiplier();
43 
44  C_INT64 GCD1 = abs64(PosMult);
45  C_INT64 GCD2 = abs64(NegMult);
46 
47  // Divide PosMult and NegMult by GCD(PosMult, NegMult);
48  CBitPatternTreeMethod::GCD(GCD1, GCD2);
49 
50  if (GCD1 != 1)
51  {
52  PosMult /= GCD1;
53  NegMult /= GCD1;
54  }
55 
56  // -1 is used to identify the start of the GCD search.
57  GCD1 = -1;
58 
59  mReaction.resize(pPositive->mReaction.size());
60  std::vector< C_INT64 >::iterator it = mReaction.begin();
61  std::vector< C_INT64 >::iterator end = mReaction.end();
62 
63  std::vector< C_INT64 >::const_iterator itPos = pPositive->mReaction.begin();
64  std::vector< C_INT64 >::const_iterator itNeg = pNegative->mReaction.begin();
65 
66  for (; it != end; ++it, ++itPos, ++itNeg)
67  {
68  // TODO We need to check that we do not have numerical overflow
69  *it = PosMult * *itPos + NegMult * *itNeg;
70 
71  if (*it == 0 || GCD1 == 1)
72  {
73  continue;
74  }
75 
76  if (GCD1 == -1)
77  {
78  GCD1 = abs64(*it);
79  continue;
80  }
81 
82  GCD2 = abs64(*it);
83 
84  CBitPatternTreeMethod::GCD(GCD1, GCD2);
85  }
86 
87  if (GCD1 > 1)
88  {
89  for (it = mReaction.begin(); it != end; ++it)
90  {
91  *it /= GCD1;
92  }
93  }
94 }
95 
97 {
98  assert(*mIterator = this);
99 
100  *mIterator = NULL;
101 }
102 
104 {
105  return mZeroSet;
106 }
107 
108 std::vector< C_INT64 > & CStepMatrixColumn::getReaction()
109 {
110  return mReaction;
111 }
112 
114 {
115  size_t Size = mZeroSet.getNumberOfBits();
116  indexes.resize(Size);
117  size_t * pIndex = indexes.array();
118 
119  CZeroSet::CIndex Index;
120  size_t i = 0;
121  size_t imax = Size - mReaction.size();
122 
123  for (; i < imax; ++i, ++Index)
124  {
125  if (!mZeroSet.isSet(Index))
126  {
127  *pIndex = i;
128  pIndex++;
129  }
130  }
131 
132  for (i = mReaction.size(); i > 0;)
133  {
134  --i;
135 
136  if (mReaction[i] != 0)
137  {
138  *pIndex = (mReaction.size() - i - 1) + imax;
139  pIndex++;
140  }
141  }
142 
143  Size = pIndex - indexes.array();
144  indexes.resize(Size, true);
145 
146  //DebugFile << *this << std::endl;
147  //DebugFile << "@CSMC: " << indexes << std::endl;
148 }
149 
151 {
152  mReaction.insert(mReaction.begin(), value);
153 }
154 
156 {
157  mReaction.pop_back();
158 }
159 
160 std::ostream & operator << (std::ostream & os, const CStepMatrixColumn & c)
161 {
162  os << ' ';
163 
164  size_t Size = c.mZeroSet.getNumberOfBits();
165  CZeroSet::CIndex Index;
166  size_t i = 0;
167  size_t imax = Size - c.mReaction.size();
168 
169  for (; i < imax; ++i, ++Index)
170  {
171  if (c.mZeroSet.isSet(Index))
172  {
173  os << "*\t";
174  }
175  else
176  {
177  os << ".\t";
178  }
179  }
180 
181  for (i = c.mReaction.size(); i > 0;)
182  {
183  --i;
184 
185  os << c.mReaction[i] << "\t";;
186  }
187 
188  return os;
189 }
size_t getNumberOfBits() const
Definition: CZeroSet.h:95
void push_front(const C_INT64 &value)
static void GCD(C_INT64 &m, C_INT64 &n)
const CZeroSet & getZeroSet() const
std::vector< C_INT64 > mReaction
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
#define C_INT64
Definition: copasi.h:88
#define abs64
Definition: copasi.h:94
bool isSet(const CIndex &index) const
Definition: CZeroSet.h:80
void getAllUnsetBitIndexes(CVector< size_t > &indexes) const
CStepMatrixColumn ** mIterator
CType * array()
Definition: CVector.h:139
CStepMatrixColumn(const size_t &size=0)
const C_INT64 & getMultiplier() const
std::vector< C_INT64 > & getReaction()
std::ostream & operator<<(std::ostream &os, const CStepMatrixColumn &c)