COPASI API  4.16.103
CInterpolation.cpp
Go to the documentation of this file.
1 // Copyright (C) 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 //#include "copasi.h"
7 //#include "copasi/utilities/CVersion.h"
8 //#include "copasi/utilities/CMatrix.h"
9 //#include "copasi/utilities/CCopasiVector.h"
10 
11 #include "CInterpolation.h"
12 #include <vector>
13 #include <set>
14 #include <new>
15 #include <iostream>
16 
17 /*******************************/
18 /* Code for class CStateRecord */
19 /*******************************/
20 
22  mpValues(NULL),
23  mSize(0)
24 {}
25 
26 CStateRecord::CStateRecord(size_t speciesNum):
27  mpValues(new C_FLOAT64[speciesNum + 2]),
28  mSize(speciesNum + 2)
29 {}
30 
32  mpValues(new C_FLOAT64[src.mSize]),
33  mSize(src.mSize)
34 {
35  memcpy(mpValues, src.mpValues, sizeof(C_FLOAT64)*mSize);
36 }
37 
39 {
41 }
42 
44 {
45  if (this != &rhs)
46  {
47  if (mSize != rhs.mSize)
48  {
50  mpValues = new C_FLOAT64[rhs.mSize];
51  mSize = rhs.mSize;
52  }
53 
54  memcpy(mpValues, rhs.mpValues, sizeof(C_FLOAT64)*mSize);
55  }
56 
57  return *this;
58 }
59 
61 {
62  return mpValues;
63 }
64 
66 {
67  return mpValues[0];
68 }
69 
71 {
72  return mpValues[mSize - 1];
73 }
74 
75 const size_t & CStateRecord::getArrayLen() const
76 {
77  return mSize;
78 }
79 
80 //Friend Function
81 void printRecord(const CStateRecord& record)
82 {
83  for (size_t i = 0; i < record.mSize; i++)
84  {
85  std::cout << record.mpValues[i] << " ";
86  }
87 
88  std::cout << std::endl;
89  return;
90 }
91 
92 /************************************/
93 /* Code for Class of CInterpolation */
94 /************************************/
96  mpX(NULL),
97  mpY(NULL),
98  mStateNum(0),
99  mBeginRecordIndex(0),
100  mpState(NULL),
101  mHaveRecordNum(0),
102  mpInterpolationState(NULL),
103  mShiftBeginIndex(false)
104 {}
105 
106 CInterpolation::CInterpolation(size_t stateNum, size_t speciesNum):
107  mpX(new C_FLOAT64[stateNum]),
108  mpY(new C_FLOAT64[stateNum]),
109  mStateNum(stateNum),
110  mBeginRecordIndex(0),
111  mHaveRecordNum(0),
112  mShiftBeginIndex(false)
113 {
114  mpInterpolationState = new CStateRecord(speciesNum);
115 
116  mpState = static_cast< CStateRecord* >
117  (operator new[](sizeof(CStateRecord) * stateNum));
118 
119  for (size_t i = 0; i < stateNum; i++)
120  new(&mpState[i]) CStateRecord(speciesNum);
121 
122  return;
123 }
124 
126 {
127  pdeletev(mpX);
128  pdeletev(mpY);
129 
131 
132  for (int i = mStateNum - 1; i >= 0; i--)
133  {
134  mpState[i].~CStateRecord();
135  }
136 
137  if (mpState != NULL)
138  operator delete[](mpState);
139 }
140 
142  const C_FLOAT64 *y)
143 {
144  mHaveRecordNum++;
145 
147  {
149 
150  if (!mShiftBeginIndex)
151  mShiftBeginIndex = true;
152  }
153 
154  if (mShiftBeginIndex)
156 
157  size_t newStateIndex = mBeginRecordIndex
158  + mHaveRecordNum - 1;
159 
160  newStateIndex %= mStateNum;
161  //std::cout << "newStateIndex: " << newStateIndex << std::endl;
162 
163  C_FLOAT64 *pArray = mpState[newStateIndex].getArray();
164  const size_t size = mpState[newStateIndex].getArrayLen();
165 
166  pArray[0] = time;
167 
168  for (size_t i = 1; i < size; i++, y++)
169  {
170  ++pArray;
171  *pArray = *y;
172  }
173 
174  return;
175 }
176 
178 {
179  mBeginRecordIndex = 0;
180  mHaveRecordNum = 0;
181  mShiftBeginIndex = false;
182 }
183 
185 {
186  //Before use this function, data have been stored in
187  //arraies mpX, mpY. Lengths of array is given by
188  //mHaveRecordNum
189 
190  if ((x > mpX[mHaveRecordNum - 1])
191  || (x < mpX[0])) //If x is within the range
192  //Show Error !!!
193  std::cout << "x " << x << " is out of range" << std::endl;
194 
195  if (mHaveRecordNum == 1)
196  return mpY[0]; //only one point recorded
197 
198  size_t i, j;
199 
200  for (i = 0; i < mHaveRecordNum; i++) //x is in *mpX
201  {
202  if (x == mpX[i])
203  return mpY[i];
204  }
205 
206  //Calculate interpolation value
207  C_FLOAT64 numer, tmpNumer = 1;
208  C_FLOAT64 denom;
209 
210  for (i = 0; i < mHaveRecordNum; i++)
211  tmpNumer *= x - mpX[i];
212 
213  C_FLOAT64 y = 0;
214 
215  for (i = 0; i < mHaveRecordNum; i++)
216  {
217  numer = tmpNumer / (x - mpX[i]);
218  denom = 1;
219 
220  for (j = 0; j < mHaveRecordNum; j++)
221  {
222  if (i != j)
223  denom *= mpX[i] - mpX[j];
224  }
225 
226  y += mpY[i] * numer / denom;
227  }
228 
229  return y;
230 }
231 
233 {
234  size_t i = mBeginRecordIndex;
235 
236  //(1)Load time into mpY
237  //(2)Load sum of slow reaction propensities
238  // into mpX
239 
240  for (size_t j = 0; j < mHaveRecordNum; j++)
241  {
242  mpY[j] = mpState[i].getTime();
243  mpX[j] = mpState[i].getSlowPropensitySum();
244 
245  if (++i > (mStateNum - 1))
246  i -= mStateNum;
247  }
248 
249  //(3)do interpolation at propensity=0
251  pArray[0] = interpolation(0.0);
252 }
253 
255 {
256  C_FLOAT64 *pInterArray = mpInterpolationState->getArray();
257  C_FLOAT64 **ppArray = new C_FLOAT64*[mHaveRecordNum];
259 
260  //(1)Load time into array mpX
261  //(2)Load pointers of objects CStateRecord into ppArray
262  size_t i = mBeginRecordIndex;
263 
264  for (size_t j = 0; j < mHaveRecordNum; j++)
265  {
266  ppArray[j] = mpState[i].getArray();
267  mpX[j] = mpState[i].getTime();
268 
269  if (++i > mStateNum - 1)
270  i -= mStateNum;
271  }
272 
273  //(3)Load propensities into mpY
274  //(4)Do interpolation
275  size_t speciesNum = (mpInterpolationState->getArrayLen()) - 2;
276 
277  for (size_t species = 1; species <= speciesNum; species++)
278  {
279  for (size_t j = 0; j < mHaveRecordNum; j++)
280  mpY[j] = ppArray[j][species];
281 
282  pInterArray[species] = interpolation(time);
283  }
284 }
285 
287 {
289  calculateState();
290  return mpInterpolationState;
291 }
292 
293 //Friend Function
295 {
296  std::cout << "State Num: " << interp.mStateNum << std::endl;
297  std::cout << "Record Num: " << interp.mHaveRecordNum << std::endl;
298  std::cout << "Begin Record Index: " << interp.mBeginRecordIndex << std::endl;
299  std::cout << "Shift or not: ";
300 
301  if (interp.mShiftBeginIndex)
302  std::cout << "Yes" << std::endl;
303  else
304  std::cout << "No" << std::endl;
305 
306  int site;
307 
308  for (size_t i = 0; i < interp.mHaveRecordNum; i++)
309  {
310  std::cout << "Record #:" << i << std::endl;
311  site = (interp.mBeginRecordIndex + i) % interp.mStateNum;
312  printRecord(interp.mpState[site]);
313  }
314 
315  std::cout << "Finish, aha" << std::endl;
316  return;
317 }
CStateRecord * mpInterpolationState
CStateRecord & operator=(const CStateRecord &rhs)
#define pdelete(p)
Definition: copasi.h:215
void printInterpolation(const CInterpolation &interp)
void recordState(const C_FLOAT64 time, const C_FLOAT64 *y)
CStateRecord * mpState
size_t mBeginRecordIndex
const C_FLOAT64 & getSlowPropensitySum() const
void printRecord(const CStateRecord &record)
#define pdeletev(p)
Definition: copasi.h:216
void calculateFiringTime()
C_FLOAT64 interpolation(C_FLOAT64 x)
C_FLOAT64 * getArray() const
const C_FLOAT64 & getTime() const
CStateRecord * getInterpolationState()
C_FLOAT64 * mpY
#define C_FLOAT64
Definition: copasi.h:92
C_FLOAT64 * mpValues
C_FLOAT64 * mpX
const size_t & getArrayLen() const