COPASI API  4.16.103
CStepMatrix.cpp
Go to the documentation of this file.
1 // Begin CVS Header
2 // $Source: /Volumes/Home/Users/shoops/cvs/copasi_dev/copasi/elementaryFluxModes/CStepMatrix.cpp,v $
3 // $Revision: 1.13 $
4 // $Name: $
5 // $Author: heilmand $
6 // $Date: 2010/08/02 15:12:41 $
7 // End CVS Header
8 
9 // Copyright (C) 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 "copasi.h"
20 
21 #include "CStepMatrix.h"
22 #include "CStepMatrixColumn.h"
23 
24 #include "utilities/CMatrix.h"
25 
27  CVector< CStepMatrixColumn * >(0),
28  mRows(0),
29  mPivot(0),
30  mFirstUnconvertedRow(0)
31 {}
32 
34  CVector< CStepMatrixColumn * >(0),
35  mRows(rows),
36  mPivot(rows),
37  mFirstUnconvertedRow(0)
38 {
39  size_t * pPivot = mPivot.array();
40 
41  for (size_t i = 0; i < mRows; ++i, ++pPivot)
42  {
43  *pPivot = i;
44  }
45 }
46 
48  CVector< CStepMatrixColumn * >(0),
49  mRows(nullspaceMatrix.numRows()),
50  mPivot(nullspaceMatrix.numRows()),
51  mFirstUnconvertedRow(0)
52 {
53  size_t Cols = nullspaceMatrix.numCols();
54 
55  resize(Cols);
56  iterator it = array();
57  mInsert = mBeyond = it + Cols;
58 
59  CVector< CStepMatrixColumn * > Columns(Cols);
60  CStepMatrixColumn ** pColumn = Columns.array();
61  CStepMatrixColumn ** pColumnEnd = pColumn + Cols;
62 
63  for (; pColumn != pColumnEnd; ++pColumn, ++it)
64  {
65  *pColumn = new CStepMatrixColumn(mRows);
66 
67  (*pColumn)->setIterator(it);
68  *it = *pColumn;
69  }
70 
71  size_t i;
72  size_t j;
73  const C_INT64 * pValue = nullspaceMatrix.array();
74  size_t * pPivot = mPivot.array();
75 
76  std::vector< size_t > NegativeRows;
77  bool hasNegative;
78  bool hasPositive;
79 
80  for (i = 0; i < mRows; ++i, ++pPivot)
81  {
82  *pPivot = i;
83 
84  hasNegative = false;
85  hasPositive = false;
86 
87  for (j = 0; j < Cols; ++j, ++pValue)
88  {
89  if (*pValue > 0)
90  {
91  hasPositive = true;
92  }
93  else if (*pValue < 0)
94  {
95  hasNegative = true;
96  }
97  }
98 
99  if ((!hasNegative && hasPositive))
100  {
101  convertRow(i, nullspaceMatrix);
102  }
103  }
104 
105  // We need to add the information of the unconverted rows of nullspace matrix
106  // to the columns
107 
108  if (nullspaceMatrix.size() != 0 &&
110  {
111  pValue = & nullspaceMatrix(mFirstUnconvertedRow, 0);
112  }
113  else
114  {
115  pValue = NULL;
116  }
117 
118  for (i = mFirstUnconvertedRow; i < mRows; ++i)
119  {
120  pColumn = Columns.array();
121 
122  for (j = 0; j < Cols; ++j, ++pValue, ++pColumn)
123  {
124  (*pColumn)->push_front(*pValue);
125  }
126  }
127 }
128 
130 {
131  iterator it = array();
132 
133  for (; it != mInsert; ++it)
134  {
135  if (*it != NULL)
136  {
137  delete *it;
138  }
139  }
140 }
141 
143 {
145 
146  iterator it = array();
147 
148  for (; it != mInsert; ++it)
149  {
150  assert(*it != NULL);
151 
152  assert((*it)->getMultiplier() >= 0);
153 
154  if ((*it)->getMultiplier() > 0)
155  {
156  (*it)->unsetBit(Index);
157  }
158 
159  (*it)->truncate();
160  }
161 
163 }
164 
165 /*
166 void CStepMatrix::convertMatrix() {
167 
168 }
169 */
170 
172 {
173  return mFirstUnconvertedRow;
174 }
175 
177 {
178  return mRows - mFirstUnconvertedRow;
179 }
180 
182  CStepMatrixColumn const * pPositive,
183  CStepMatrixColumn const * pNegative)
184 {
185  CStepMatrixColumn * pColumn = new CStepMatrixColumn(set, pPositive, pNegative);
186 
187  add(pColumn);
188 
189  return pColumn;
190 }
191 
193 {
194  delete pColumn;
195 }
196 
197 bool CStepMatrix::splitColumns(std::vector< CStepMatrixColumn * > & PositiveColumns,
198  std::vector< CStepMatrixColumn * > & NegativeColumns,
199  std::vector< CStepMatrixColumn * > & NullColumns)
200 {
201  assert(PositiveColumns.empty());
202  assert(NegativeColumns.empty());
203  assert(NullColumns.empty());
204 
205  iterator it = array();
206 
207  for (; it != mInsert; ++it)
208  {
209  assert(*it != NULL);
210 
211  const C_INT64 & Value = (*it)->getMultiplier();
212 
213  if (Value > 0)
214  {
215  PositiveColumns.push_back(*it);
216  }
217  else if (Value < 0)
218  {
219  NegativeColumns.push_back(*it);
220  }
221  else
222  {
223  NullColumns.push_back(*it);
224  }
225  }
226 
227  if (NegativeColumns.empty())
228  {
229  // We do not have any linear combinations, thus we can convert immediately
230  convertRow();
231 
232  return false;
233  }
234 
235  return true;
236 }
237 
238 void CStepMatrix::removeInvalidColumns(std::vector< CStepMatrixColumn * > & invalidColumns)
239 {
240  std::vector< CStepMatrixColumn * >::iterator it = invalidColumns.begin();
241  std::vector< CStepMatrixColumn * >::iterator end = invalidColumns.end();
242 
243  for (; it != end; ++it)
244  {
245  removeColumn(*it);
246  }
247 }
248 
250  CVector<size_t> & indexes) const
251 {
252  pColumn->getAllUnsetBitIndexes(indexes);
253 
254  // Apply the QR pivot
255  size_t * pIndex = indexes.array();
256  size_t * pIndexEnd = pIndex + indexes.size();
257 
258  for (; pIndex != pIndexEnd; ++pIndex)
259  *pIndex = mPivot[*pIndex];
260 
261  //DebugFile << "@CSM: " << indexes << std::endl;
262 }
263 
265  CVector< size_t > & indexes) const
266 {
267  const CZeroSet & ZeroSet = pColumn->getZeroSet();
268 
269  indexes.resize(ZeroSet.getNumberOfUnsetBits());
270  size_t * pIndex = indexes.array();
271  size_t * pIndexEnd = pIndex + indexes.size();
272 
273  CZeroSet::CIndex Bit = 0;
274  size_t Index = 0;
275 
276  for (; pIndex != pIndexEnd; ++Bit, ++Index)
277  {
278  if (!ZeroSet.isSet(Bit))
279  {
280  // Apply pivot.
281  *pIndex = mPivot[Index];
282  pIndex++;
283  }
284  }
285 
286  return;
287 }
288 
290 {
291  return array();
292 }
293 
295 {
296  return mInsert;
297 }
298 
299 void CStepMatrix::convertRow(const size_t & index,
300  CMatrix< C_INT64 > & nullspaceMatrix)
301 {
303 
304  iterator it = array();
305 
306  C_INT64 * pValue = & nullspaceMatrix(index, 0);
307 
308  if (mFirstUnconvertedRow != index)
309  {
310  C_INT64 * pFirstUnconvertedValue = & nullspaceMatrix(mFirstUnconvertedRow, 0);
311 
312  for (; it != mInsert; ++it, ++pValue, ++pFirstUnconvertedValue)
313  {
314  // At this point the matrix is compact since no columns have been destroyed,
315  // i.e., we do not need to check whether *it != NULL.
316  assert(*pValue >= 0);
317 
318  if (*pValue > 0)
319  {
320  (*it)->unsetBit(Index);
321  }
322 
323  *pValue = *pFirstUnconvertedValue;
324  }
325 
326  // We need to remember the reordering.
327  size_t tmp = mPivot[index];
330  }
331  else
332  {
333  for (; it != mInsert; ++it, ++pValue)
334  {
335  if (*pValue != 0)
336  {
337  (*it)->unsetBit(Index);
338  }
339  }
340  }
341 
343 }
344 
346 {
347  iterator from = array();
348  iterator to = array();
349 
350  for (; from != mInsert; ++from)
351  {
352  if (*from != NULL)
353  {
354  (*from)->setIterator(to);
355  *to = *from;
356 
357  ++to;
358  }
359  }
360 
361  mInsert = to;
362 }
363 
364 std::ostream & operator << (std::ostream & os, const CStepMatrix & m)
365 {
366  os << m.mPivot << std::endl;
367 
368  CZeroSet::CIndex Bit;
369 
372 
373  for (it = m.begin(); it != end; ++it)
374  {
375  os << **it << std::endl;
376  }
377 
378  return os;
379 }
void removeInvalidColumns(std::vector< CStepMatrixColumn * > &invalidColumns)
void push_front(const C_INT64 &value)
void getUnsetBitIndexes(const CStepMatrixColumn *pColumn, CVector< size_t > &indexes) const
size_t getNumberOfUnsetBits() const
Definition: CZeroSet.h:90
size_t getFirstUnconvertedRow() const
const CZeroSet & getZeroSet() const
void getAllUnsetBitIndexes(const CStepMatrixColumn *pColumn, CVector< size_t > &indexes) const
CStepMatrixColumn *const * const_iterator
Definition: CStepMatrix.h:29
void add(CStepMatrixColumn *pColumn)
Definition: CStepMatrix.h:50
std::ostream & operator<<(std::ostream &os, const CStepMatrix &m)
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
#define C_INT64
Definition: copasi.h:88
size_t mFirstUnconvertedRow
Definition: CStepMatrix.h:134
const_iterator begin() const
bool splitColumns(std::vector< CStepMatrixColumn * > &PositiveColumns, std::vector< CStepMatrixColumn * > &NegativeColumns, std::vector< CStepMatrixColumn * > &NullColumns)
size_t mRows
Definition: CStepMatrix.h:130
size_t getNumUnconvertedRows() const
const_iterator end() const
void removeColumn(CStepMatrixColumn *pColumn)
iterator mInsert
Definition: CStepMatrix.h:136
void convertRow()
bool isSet(const CIndex &index) const
Definition: CZeroSet.h:80
size_t size() const
Definition: CVector.h:100
void getAllUnsetBitIndexes(CVector< size_t > &indexes) const
CType * array()
Definition: CVector.h:139
virtual size_t size() const
Definition: CMatrix.h:132
void compact()
iterator mBeyond
Definition: CStepMatrix.h:138
CVector< size_t > mPivot
Definition: CStepMatrix.h:132
virtual size_t numCols() const
Definition: CMatrix.h:144
CStepMatrixColumn * addColumn(const CZeroSet &set, const CStepMatrixColumn *pPositive, const CStepMatrixColumn *pNegative)
virtual CType * array()
Definition: CMatrix.h:337