COPASI API  4.16.103
CIndexedPriorityQueue.cpp
Go to the documentation of this file.
1 // Begin CVS Header
2 // $Source: /Volumes/Home/Users/shoops/cvs/copasi_dev/copasi/utilities/CIndexedPriorityQueue.cpp,v $
3 // $Revision: 1.19 $
4 // $Name: $
5 // $Author: shoops $
6 // $Date: 2011/03/07 19:34:55 $
7 // End CVS Header
8 
9 // Copyright (C) 2011 - 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 // Copyright (C) 2001 - 2007 by Pedro Mendes, Virginia Tech Intellectual
20 // Properties, Inc. and EML Research, gGmbH.
21 // All rights reserved.
22 
23 #include "copasi.h"
24 #include "CCopasiMessage.h"
25 #include "CIndexedPriorityQueue.h"
26 
28 {}
29 
31 {}
32 
34 {
35  return mHeap[0].mKey;
36 }
37 
39 {
40  return mHeap[0].mIndex;
41 }
42 
43 // juergen: added 26 July, 2002
45 {
46  size_t t;
47 
48  // check if index is valid
49  if (index >= mIndexPointer.size()) return C_INVALID_INDEX;
50 
51  if ((mIndexPointer[index] != C_INVALID_INDEX) && (mIndexPointer[index] != (mHeap.size() - 1))) // if the node with the given index exists in the tree
52  {// remove the node with the given index from the tree
53  swapNodes(t = mIndexPointer[index], mHeap.size() - 1);
54  mHeap.pop_back();
55  mIndexPointer[index] = -1;
56  heapify(t);
57  }
58  else if (mIndexPointer[index] == (mHeap.size() - 1)) // last node in the heap
59  {
60  mHeap.pop_back();
61  mIndexPointer[index] = -1;
62  }
63 
64  return 0;
65 }
66 
67 // juergen: added 26 July, 2002
68 size_t CIndexedPriorityQueue::insertStochReaction(const size_t index, const C_FLOAT64 key)
69 {
70  size_t pos;
71 
72  // check if index is valid
73  if (index >= mIndexPointer.size()) return - 1;
74 
75  // first the node is inserted at the end of the heap
76  mIndexPointer[index] = mHeap.size();
77  PQNode heap_node(index, key);
78  mHeap.push_back(heap_node);
79  // bubble the node up the tree to the right position !
80  pos = mIndexPointer[index];
81 
82  while ((pos > 0) && (mHeap[parent(pos)].mKey > key))
83  {
84  swapNodes(pos, parent(pos));
85  pos = parent(pos);
86  }
87 
88  return 0;
89 }
90 
91 // juergen: added 26 July, 2002
92 void CIndexedPriorityQueue::initializeIndexPointer(const size_t numberOfReactions)
93 {
94  size_t i;
95 
96  for (i = 0; i < numberOfReactions; i++)
97  {
98  mIndexPointer.push_back(C_INVALID_INDEX);
99  }
100 }
101 
102 size_t CIndexedPriorityQueue::pushPair(const size_t index, const C_FLOAT64 key)
103 {
104  // Add an element to the priority queue. This merely pushes an item onto
105  // the back of the vector corresponding to the heap, and pushes the index
106  // onto the back of the vector corresponding to the index structure.
107  // Implicit is the assumption that this is done in contiguous ascending order
108  // for the index; i.e. in index order 0,1,2,3...
109  // N.B. The structures are not yet ordered as an indexed priority queue; this must
110  // be done using the buildHeap() method
111 
112  // First check that the index corresponds to the heap size before insertion
113  if (static_cast<unsigned int>(index) != mHeap.size())
114  {
115  CCopasiMessage(CCopasiMessage::ERROR, "Error inserting pair into priority queue");
116  return - 1;
117  }
118 
119  PQNode heap_node(index, key);
120  mHeap.push_back(heap_node);
121  // at first, position == index
122  size_t position = index; // for clarity
123  mIndexPointer.push_back(position);
124  return 0;
125 }
126 
128 {
129  for (size_t i = mHeap.size() / 2 - 1; i != C_INVALID_INDEX; i--)
130  {
131  heapify(i);
132  }
133 }
134 
136 {mHeap.clear(); mIndexPointer.clear();}
137 
138 void CIndexedPriorityQueue::updateNode(const size_t index, const C_FLOAT64 new_key)
139 {
140  size_t pos = mIndexPointer[index];
141  // cout << "Setting heap at " << pos << " to be " << new_key << endl;
142  mHeap[pos].mKey = new_key;
143  updateAux(pos);
144 }
145 
146 void CIndexedPriorityQueue::swapNodes(const size_t pos1, const size_t pos2)
147 {
148  // cout << "Swapping node " << pos1 << "(" << mHeap[pos1].mKey << ") with node ";
149  // cout << pos2 << "(" << mHeap[pos2].mKey << ")\n";
150  C_FLOAT64 tempkey = mHeap[pos1].mKey;
151  size_t index1 = mHeap[pos1].mIndex;
152  size_t index2 = mHeap[pos2].mIndex;
153  // Put the contents of the node at pos2 into the node at pos1
154  mHeap[pos1].mIndex = index2;
155  mHeap[pos1].mKey = mHeap[pos2].mKey;
156  // Put the contents formerly in the node at pos1 into the node at pos2
157  mHeap[pos2].mIndex = index1;
158  mHeap[pos2].mKey = tempkey;
159  // Swap the pointers in the index vector
160  mIndexPointer[index1] = pos2;
161  mIndexPointer[index2] = pos1;
162 }
163 
164 void CIndexedPriorityQueue::heapify(const size_t current)
165 {
166  size_t left = leftChild(current);
167  size_t right = rightChild(current);
168  size_t highest_priority = current;
169 
170  // cout << "Heapifying " << current << "(Currently " << mHeap[current].mKey << ")" << endl;
171  if ((static_cast<unsigned int>(left) < mHeap.size()) && (mHeap[left].mKey < mHeap[current].mKey))
172  {
173  highest_priority = left;
174  }
175 
176  if ((static_cast<unsigned int>(right) < mHeap.size()) && (mHeap[right].mKey < mHeap[highest_priority].mKey))
177  {
178  highest_priority = right;
179  }
180 
181  if (highest_priority != current)
182  {
183  swapNodes(current, highest_priority);
184  heapify(highest_priority);
185  }
186 }
187 
188 void CIndexedPriorityQueue::updateAux(const size_t pos)
189 {
190  size_t parent_pos = parent(pos);
191  C_FLOAT64 keyval = mHeap[pos].mKey;
192 
193  if (parent_pos != C_INVALID_INDEX &&
194  keyval < mHeap[parent_pos].mKey)
195  {
196  swapNodes(pos, parent_pos);
197  updateAux(parent_pos);
198  }
199  else
200  {
201  size_t left = leftChild(pos);
202  size_t right = rightChild(pos);
203  C_FLOAT64 min = 0.0;
204  size_t min_pos = 0;
205 
206  if (left < mHeap.size())
207  {
208  min = mHeap[left].mKey;
209  min_pos = left;
210  }
211 
212  C_FLOAT64 val; // = mHeap[right].mKey; //!!!
213 
214  if ((right < mHeap.size()) && ((val = mHeap[right].mKey) < min))
215  {
216  min = val;
217  min_pos = right;
218  }
219 
220  if ((min_pos > 0) && (keyval > min))
221  {
222  // swapNodes(mHeap[pos].mIndex, min_pos);
223  swapNodes(pos, min_pos);
224  updateAux(min_pos);
225  }
226  }
227 }
228 
229 #ifdef TEST_PRIORITY_QUEUE
230 
231 int main(int argc, char **argv)
232 {
233  // Generates a vector of pairs, with a size given by the first argument. Each element is added to a priority queue.
234  if (argc != 2)
235  {
236  cout << "Usage: " << argv[0] << " <number of pairs to generate>" << endl;
237  return - 1;
238  }
239 
240  int count = atoi(argv[1]);
241  cout << "Creating priority queue of size " << count << endl;
242  std::vector<C_FLOAT64> invec;
244  CRandom *rand = new CRandom(1);
245  C_FLOAT64 rndval;
246  cout << "Input vector:\n";
247 
248  for (int i = 0; i < count; i++)
249  {
250  rndval = rand->getUniformRandom();
251  invec.push_back(rndval);
252  cout << "element " << i << ":" << rndval << endl;
253  pq.pushPair(i, invec[i]);
254  }
255 
256  cout << "Building heap\n";
257  pq.buildHeap();
258  cout << "Done building heap\n";
259  // Display the priority queue
260  cout << "\nPriority Queue:\n";
261 
262  for (int i = 0; i < count; i++)
263  {
264  cout << "Queue: ";
265 
266  for (int j = 0; j < count; j++) cout << " " << j << "-" << pq[j];
267 
268  cout << endl;
269  cout << "Position: " << i;
270  cout << " Index = " << pq.topIndex();
271  cout << " Key = " << pq.topKey() << endl;
272  pq.updateNode(pq.topIndex(), 10000000);
273  }
274 
275  return 0;
276 }
277 
278 #endif // TEST_PRIORITY_QUEUE
279 
280 std::ostream & operator<<(std::ostream &os, const PQNode & d)
281 {
282  os << "(" << d.mIndex << ", " << d.mKey << ")";
283  return os;
284 }
285 
286 std::ostream & operator<<(std::ostream &os, const CIndexedPriorityQueue & d)
287 {
288  size_t i;
289 
290  os << "PQ: " << std::endl;
291 
292  std::vector <PQNode>::const_iterator it;
293  os << " mHeap: " << std::endl;
294 
295  for (it = d.mHeap.begin(); it != d.mHeap.end(); it++)
296  os << *it << std::endl;
297 
298  os << " mIndexPointer: " << std::endl;
299 
300  for (i = 0; i < d.mIndexPointer.size(); i++)
301  os << d.mIndexPointer[i] << " ";
302 
303  os << std::endl;
304 
305  os << std::endl;
306 
307  return os;
308 }
std::vector< size_t > mIndexPointer
int main(int argc, char **argv)
#define C_INVALID_INDEX
Definition: copasi.h:222
std::vector< PQNode > mHeap
size_t leftChild(const size_t pos) const
std::ostream & operator<<(std::ostream &os, const PQNode &d)
size_t insertStochReaction(const size_t index, const C_FLOAT64 key)
void heapify(const size_t pos)
void updateAux(const size_t position)
size_t removeStochReaction(const size_t index)
void updateNode(const size_t index, const C_FLOAT64 key)
#define C_FLOAT64
Definition: copasi.h:92
C_FLOAT64 mKey
void swapNodes(const size_t index1, const size_t index2)
size_t rightChild(const size_t pos) const
size_t parent(const size_t pos) const
size_t pushPair(const size_t index, const C_FLOAT64 key)
#define min(a, b)
Definition: f2c.h:175
void initializeIndexPointer(const size_t numberOfReactions)