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.
4 // All rights reserved.
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"
17 #include "COptTask.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 =
282  mpCallBack->addItem("Iterations",
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
COptTask * mpParentTask
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