COPASI API  4.16.103
CLayoutEngine.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2014 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 <algorithm>
7 #include <cmath>
8 #include <iostream>
9 #include "copasi.h"
10 #include "CLayoutEngine.h"
12 
14  : mpLayout(l),
15  mSecondOrder(false)
16 {
17  if (!mpLayout) return;
18 
20 
21  if (mSecondOrder) mVariables.resize(mVariables.size() * 2, 0);
22 
23  mRhs.resize(mVariables.size(), 0);
24 
25  //for RK4
26  mRhsA.resize(mVariables.size(), 0);
27  mRhsB.resize(mVariables.size(), 0);
28  mRhsC.resize(mVariables.size(), 0);
29  mVar2.resize(mVariables.size(), 0);
30 
31  //init LSODA
32  mLsodaStatus = 1;
33  //mState = 1;
34  //mJType = 2;
35  mErrorMsg.str("");
36 
37  mData.pMethod = this;
38  mData.dim = mVariables.size();
39  //mAtol = mpModel->initializeAtolVector(*mpAbsoluteTolerance, *mpReducedModel);
40  //mY = mpState->beginIndependent();
41 
42  //mRtol = *mpRelativeTolerance;
43 
44  mTime = 0;
45 
46  mDWork.resize(22 + mData.dim * std::max<C_INT>(16, mData.dim + 9) + 3 * 0);
47  mDWork[4] = mDWork[5] = mDWork[6] = mDWork[7] = mDWork[8] = mDWork[9] = 0.0;
48  mDWork[7] = 0; //minstep
49  mIWork.resize(20 + mData.dim);
50  mIWork[4] = mIWork[6] = mIWork[9] = 0;
51 
52  mIWork[5] = 100; //*mpMaxInternalSteps;
53  mIWork[7] = 12;
54  mIWork[8] = 5;
56 }
57 
58 void CLayoutEngine::calcRHS(std::vector<double> & state, double* rhs)
59 {
60  std::vector<double> forces; forces.resize(mpLayout->getNumVariables());
61  calcForces(state, forces);
62 
63  unsigned int i, imax = mpLayout->getNumVariables();
64 
65  for (i = 0; i < imax; ++i)
66  {
67  if (mSecondOrder)
68  {
69  rhs[i + imax] = (forces[i] * 0.04 - state[i + imax] * 0.05) / mpLayout->getMassVector()[i];
70  rhs[i] = state[i + imax];
71  }
72  else
73  {
74  //first order
75  rhs[i] = forces[i] / mpLayout->getMassVector()[i];
76  }
77  }
78 }
79 
80 void CLayoutEngine::calcForces(std::vector<double> & state, std::vector<double> & forces)
81 {
82  if (!mpLayout) return;
83 
84  //do numerical derivations...
85  mpLayout->setState(state);
86 
87  //assymmetric algorithm
88  /*double pot0 = mpLayout->getPotential();
89 
90  double store, newpot;
91  unsigned int i, imax = mpLayout->getNumVariables();
92  for (i=0; i<imax; ++i)
93  {
94  store = mVariables[i];
95  mVariables[i]+=1.0;
96  mpLayout->setState(mVariables);
97  newpot = mpLayout->getPotential();
98  mVariables[i] = store;
99  mForce[i] = pot0-newpot;
100  }*/
101 
102  //symmetric algorithm
103  double pot0;
104  double store, newpot;
105  unsigned int i, imax = mpLayout->getNumVariables();
106 
107  for (i = 0; i < imax; ++i)
108  {
109  store = state[i];
110  //std::cout << "var " << store;
111  state[i] -= 0.5;
112  mpLayout->setState(state);
113  pot0 = mpLayout->getPotential();
114  //std::cout << "pot0 " << pot0 << std::endl;
115  state[i] += 1.0;
116  mpLayout->setState(state);
117  newpot = mpLayout->getPotential();
118  state[i] = store;
119  forces[i] = pot0 - newpot;
120  //std::cout << "force " << i << " = " << forces[i] << std::endl;
121  }
122 }
123 
125 {
126  if (!mpLayout) return -1.0;
127 
128  unsigned int i, imax = mVariables.size();
129 
130  //Euler
131 
133  double pot = mpLayout->getPotential();
134 
135  calcRHS(mVariables, &mRhs[0]);
136 
137  double dt = 3;
138 
139  for (i = 0; i < imax; ++i)
140  {
141  mVariables[i] += mRhs[i] * dt;
142  }
143 
144  double newpot;
145 
146  for (;;)
147  {
149  newpot = mpLayout->getPotential();
150 
151  if (newpot < pot) break;
152 
153  if (dt < 1e-3) break;
154 
155  //step back
156  dt *= 0.5;
157 
158  for (i = 0; i < imax; ++i)
159  {
160  mVariables[i] -= mRhs[i] * dt;
161  }
162  }
163 
164  //std::cout << dt << " " << newpot << std::endl;
165  return newpot;
166 }
167 
169 {
170  if (!mpLayout) return;
171 
172  const double dt = 0.2;
173 
174  unsigned int i, imax = mVariables.size();
175 
176  //Euler
177  /*calcRHS(mVariables, &mRhs[0]);
178 
179  for (i = 0; i < imax; ++i)
180  {
181  mVariables[i] += mRhs[i] * 0.05;
182  }*/
183 
184  //RK4
185  calcRHS(mVariables, &mRhs[0]);
186 
187  for (i = 0; i < imax; ++i)
188  {
189  mVar2[i] = mVariables[i] + mRhs[i] * 0.5 * dt;
190  }
191 
192  calcRHS(mVar2, &mRhsA[0]);
193 
194  for (i = 0; i < imax; ++i)
195  {
196  mVar2[i] = mVariables[i] + mRhsA[i] * 0.5 * dt;
197  }
198 
199  calcRHS(mVar2, &mRhsB[0]);
200 
201  for (i = 0; i < imax; ++i)
202  {
203  mVar2[i] = mVariables[i] + mRhsB[i] * dt;
204  }
205 
206  calcRHS(mVar2, &mRhsC[0]);
207 
208  for (i = 0; i < imax; ++i)
209  {
210  mVariables[i] += dt / 6 * (mRhs[i] + 2 * mRhsA[i] + 2 * mRhsB[i] + mRhsC[i]);
211  }
212 
213  //LSODA
214  /*
215 
216  C_FLOAT64 EndTime = mTime + dt;
217  C_INT ITOL = 1; // mRtol scalar, mAtol scalar
218  C_INT one = 1;
219  long int two = 2;
220  long int three = 3;
221  double reltol = 0.1;
222  double abstol = 0.5;
223 
224  C_INT DSize = mDWork.size();
225  C_INT ISize = mIWork.size();
226 
227  //std::cout << "jhhjk" << &mVariables[0] << std::endl;
228 
229  mLSODA(&EvalF, // 1. evaluate F
230  &mData.dim, // 2. number of variables
231  &mVariables[0], // 3. the array of current concentrations
232  &mTime, // 4. the current time
233  &EndTime, // 5. the final time
234  &ITOL, // 6. error control
235  &reltol, // 7. relative tolerance array
236  &abstol, // 8. absolute tolerance array
237  &two, // 9. output by overshoot & interpolation
238  &mLsodaStatus, // 10. the state control variable
239  &one, // 11. further options (one)
240  &mDWork[0], // 12. the double work array
241  &DSize, // 13. the double work array size
242  &mIWork[0], // 14. the int work array
243  &ISize, // 15. the int work array size
244  NULL, // 16. evaluate J (not given)
245  &two); // 17. the type of jacobian calculate (2)
246 
247  if (mLsodaStatus != 2)
248  std::cout << mLsodaStatus << std::endl;
249  */
250 
252 }
253 
254 void CLayoutEngine::EvalF(const C_INT * n, const C_FLOAT64 * t, const C_FLOAT64 * y, C_FLOAT64 * ydot)
255 {static_cast<Data *>((void *) n)->pMethod->evalF(t, y, ydot);}
256 
257 void CLayoutEngine::evalF(const C_FLOAT64 * /* t */, const C_FLOAT64 * /* y */, C_FLOAT64 * ydot)
258 {
259  //std::cout << *t << " " << y << " " << &mVariables[0] << std::endl;
260 
261 // assert(y == &mVariables[0]);
262 
263  calcRHS(mVariables, ydot);
264 
265  return;
266 }
#define C_INT
Definition: copasi.h:115
void evalF(const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
std::vector< double > mRhs
Definition: CLayoutEngine.h:47
CLayoutEngine * pMethod
Definition: CLayoutEngine.h:22
virtual const std::vector< double > & getMassVector() const
void setOstream(std::ostream &os)
std::vector< double > mRhsA
Definition: CLayoutEngine.h:50
std::vector< double > mVariables
Definition: CLayoutEngine.h:46
virtual const std::vector< double > & getInitialValues() const =0
std::ostringstream mErrorMsg
Definition: CLayoutEngine.h:74
static void EvalF(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
std::vector< C_INT > mIWork
Definition: CLayoutEngine.h:69
C_FLOAT64 mTime
Definition: CLayoutEngine.h:90
CAbstractLayoutInterface * mpLayout
Definition: CLayoutEngine.h:42
std::vector< C_FLOAT64 > mDWork
Definition: CLayoutEngine.h:64
void calcRHS(std::vector< double > &state, double *rhs)
virtual double getPotential()=0
#define C_FLOAT64
Definition: copasi.h:92
std::vector< double > mRhsB
Definition: CLayoutEngine.h:51
unsigned int getNumVariables() const
virtual bool setState(const std::vector< double > &vars)=0
CLayoutEngine(CAbstractLayoutInterface *l, bool so)
void calcForces(std::vector< double > &state, std::vector< double > &forces)
C_INT mLsodaStatus
Definition: CLayoutEngine.h:79
std::vector< double > mVar2
Definition: CLayoutEngine.h:53
std::vector< double > mRhsC
Definition: CLayoutEngine.h:52
void stepIntegration()