COPASI API  4.16.103
COptMethodCoranaWalk.cpp
Go to the documentation of this file.
1 // Copyright (C) 2012 - 2013 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
5
6 // written by Pedro Mendes, January 2012
7 // adapted from the Corana Simulated Annealing algorithm
8 // this does only the top level iteration of that algorithm
9
10 #include <cmath>
11
12 #include "copasi.h"
13
14 #include "COptMethodCoranaWalk.h"
15 #include "COptProblem.h"
16 #include "COptItem.h"
18
22
23 #define STORED 2
24 #define NS 5
25 #define K 1.0
26
28  COptMethod(CCopasiTask::optimization, CCopasiMethod::CoranaWalk, pParent)
29 {
30  addParameter("Temperature", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0);
31  addParameter("Iterations", CCopasiParameter::UINT, (unsigned C_INT32) 100);
32  addParameter("Random Number Generator", CCopasiParameter::UINT, (unsigned C_INT32) CRandom::mt19937);
33  addParameter("Seed", CCopasiParameter::UINT, (unsigned C_INT32) 0);
34
35  initObjects();
36 }
37
39  const CCopasiContainer * pParent):
40  COptMethod(src, pParent)
41 {initObjects();}
42
44 {cleanup();}
45
47 {
49 }
50
52 {
53  if (!initialize())
54  {
55  if (mpCallBack)
57
58  return false;
59  }
60
61  size_t i, j;
62
63  size_t h, a;
64  C_FLOAT64 xc, p, c, New, minstep;
65  bool processing;
66
67  // set the minimum step size as being the average of step sizes
68  // or 100*DBL_EPSILON if average is zero
69  for (i = 0, minstep = 0.0; i < mVariableSize; i++)
70  {
71  minstep += mCurrent[i];
72  }
73
74  if (minstep > std::numeric_limits< C_FLOAT64 >::epsilon())
75  minstep /= i;
76  else
77  minstep = 100 * std::numeric_limits< C_FLOAT64 >::epsilon();
78
79  // initial point is first guess but we have to make sure that we
80  // are within the parameter domain
81  for (i = 0; i < mVariableSize; i++)
82  {
83  const COptItem & OptItem = *(*mpOptItem)[i];
84
85  switch (OptItem.checkConstraint(OptItem.getStartValue()))
86  {
87  case - 1:
88  mCurrent[i] = *OptItem.getLowerBoundValue();
89  break;
90
91  case 1:
92  mCurrent[i] = *OptItem.getUpperBoundValue();
93  break;
94
95  case 0:
96  mCurrent[i] = OptItem.getStartValue();
97  break;
98  }
99
100  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
101
102  // The step must not contain any zeroes
103  mStep[i] = std::max(fabs(mCurrent[i]), minstep);
104  }
105
106  // find the objective function value at the start
108  // count it
109  mCurrentIteration = 1;
110
111  if (!isnan(mEvaluationValue))
112  {
113  // and store that value
116
117  // We found a new best value lets report it.
119
120  if (mpCallBack)
122  }
123
124  mAccepted = 0;
125
126  processing = true; // we want to do some work
127
128  do
129  {
130  for (j = 0; j < NS && mContinue && processing; j++) // adjustment in all directions
131  {
132  for (h = 0; h < mVariableSize && mContinue && processing; h++) // adjustment in one direction
133  {
134  // Calculate the step
135  xc = (2.0 * mpRandom->getRandomCC() - 1) * mStep[h];
136  New = mCurrent[h] + xc;
137
138  // Set the new parameter value
139  (*(*mpSetCalculateVariable)[h])(New);
140
141  // Check all parametric constraints
143  {
144  // Undo since not accepted
145  (*(*mpSetCalculateVariable)[h])(mCurrent[h]);
146  continue;
147  }
148
149  // evaluate the function
150  evaluate();
151  // count it
153
154  if (mpCallBack)
156
157  // Check all functional constraints
159  {
160  // Undo since not accepted
161  (*(*mpSetCalculateVariable)[h])(mCurrent[h]);
162
163  continue;
164  }
165
166  // here we have a completely feasible point
167  // keep if energy is reduced
169  {
170  // only one value has changed...
171  mCurrent[h] = New;
173  i++; // a new point
174  mAccepted[h]++; // a new point in this coordinate
175
177  {
178  // and store that value
181
182  // We found a new best value lets report it.
184  }
185  }
186  else
187  {
188  // keep with probability p, if energy is increased
189  p = exp((mCurrentValue - mEvaluationValue) / (K * mTemperature));
190
191  if (p > mpRandom->getRandomCO())
192  {
193  // only one value has changed...
194  mCurrent[h] = New;
196  i++; // a new point
197  mAccepted[h]++; // a new point in this coordinate
198  }
199  else
200  // Undo since not accepted
201  (*(*mpSetCalculateVariable)[h])(mCurrent[h]);
202  }
203
204  // check if it is time to stop
205  if (mCurrentIteration == mIterations) processing = false;
206  }
207  }
208
209  // update the step sizes
210  for (a = 0; a < mVariableSize; a++)
211  {
212  c = (C_FLOAT64) mAccepted[a] / (C_FLOAT64) NS;
213
214  if (c > 0.6)
215  mStep[a] *= 1 + 5 * (c - 0.6);
216  else if (c < 0.4)
217  mStep[a] /= 1 + 5 * (0.4 - c);
218
219  mAccepted[a] = 0;
220  }
221
222  // the equilibrium energy
223
224  if (processing)
225  {
227
228  for (a = 0; a < mVariableSize; a++)
229  (*(*mpSetCalculateVariable)[a])(mCurrent[a]);
230
232  }
233
234  if (mpCallBack)
236  }
237  while (processing && mContinue);
238
239  if (mpCallBack)
241
242  return true;
243 }
244
246 {
247  return true;
248 }
249
251 {
252  // We do not need to check whether the parametric constraints are fulfilled
253  // since the parameters are created within the bounds.
254
257
258  // When we leave the either functional domain
259  // we set the objective value +Inf
261  mEvaluationValue = std::numeric_limits<C_FLOAT64>::infinity();
262
263  return mEvaluationValue;
264 }
265
267 {
268  cleanup();
269
270  if (!COptMethod::initialize()) return false;
271
272  mTemperature = * getValue("Temperature").pUDOUBLE;
273  mIterations = * getValue("Iterations").pUINT;
274  mpRandom =
275  CRandom::createGenerator(* (CRandom::Type *) getValue("Random Number Generator").pUINT,
276  * getValue("Seed").pUINT);
277
278  mCurrentIteration = 0;
279
280  if (mpCallBack)
281  mhIterations =
284
285  mBestValue = std::numeric_limits<C_FLOAT64>::infinity();
286  mContinue = true;
287
288  mVariableSize = mpOptItem->size();
289
293
294  return true;
295 }
unsigned C_INT32 mIterations
virtual C_INT32 checkConstraint() const
Definition: COptItem.cpp:401
Definition: COptMethod.h:58
virtual bool initialize()
Definition: COptMethod.cpp:189
#define K
COptMethodCoranaWalk(const COptMethodCoranaWalk &src, const CCopasiContainer *pParent=NULL)
unsigned C_INT32 mCurrentIteration
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
COptProblem * mpOptProblem
Definition: COptMethod.h:56
virtual void output(const COutputInterface::Activity &activity)
#define C_INT32
Definition: copasi.h:90
virtual bool progressItem(const size_t &handle)
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
virtual bool calculate()
virtual C_FLOAT64 getRandomCC()
Definition: CRandom.cpp:235
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
const std::vector< UpdateMethod * > * mpSetCalculateVariable
Definition: COptMethod.h:65
const Value & getValue() const
const C_FLOAT64 * getLowerBoundValue() const
Definition: COptItem.h:191
CVector< size_t > mAccepted
const CVector< C_FLOAT64 > & getSolutionVariables() const
unsigned C_INT32 * pUINT
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
virtual bool finishItem(const size_t &handle)
virtual bool checkFunctionalConstraints()
const C_FLOAT64 & getStartValue() const
Definition: COptItem.cpp:199
const std::vector< COptItem * > * mpOptItem
Definition: COptMethod.h:70
virtual C_FLOAT64 getRandomCO()
Definition: CRandom.cpp:245
#define C_FLOAT64
Definition: copasi.h:92
CVector< C_FLOAT64 > mCurrent
const C_FLOAT64 * getUpperBoundValue() const
Definition: COptItem.h:198
CVector< C_FLOAT64 > mStep
bool addParameter(const CCopasiParameter &parameter)
const C_FLOAT64 & evaluate()
virtual bool checkParametricConstraints()
#define NS
const C_FLOAT64 & getCalculateValue() const
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
CProcessReport * mpCallBack
#define max(a, b)
Definition: f2c.h:176