COPASI API  4.16.103
COptMethodSA.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 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 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2003 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 // adapted by Stefan Hoops, August 2006, from the original Gepasi file
16 // nelmea.cpp :
17 
18 #include <cmath>
19 
20 #include "copasi.h"
21 
22 #include "COptMethodSA.h"
23 #include "COptProblem.h"
24 #include "COptItem.h"
25 #include "COptTask.h"
26 
30 
31 #define STORED 2
32 #define NS 10
33 #define K 1.0
34 
36  COptMethod(CCopasiTask::optimization, CCopasiMethod::SimulatedAnnealing, pParent)
37 {
38  addParameter("Start Temperature", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0);
39  addParameter("Cooling Factor", CCopasiParameter::UDOUBLE, (C_FLOAT64) 0.85);
40  addParameter("Tolerance", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.e-006);
41  addParameter("Random Number Generator", CCopasiParameter::UINT, (unsigned C_INT32) CRandom::mt19937);
42  addParameter("Seed", CCopasiParameter::UINT, (unsigned C_INT32) 0);
43 
44  initObjects();
45 }
46 
48  const CCopasiContainer * pParent):
49  COptMethod(src, pParent)
50 {initObjects();}
51 
53 {cleanup();}
54 
56 {
58 }
59 
61 {
62  if (!initialize())
63  {
64  if (mpCallBack)
66 
67  return false;
68  }
69 
70  size_t i, j, k, m;
71 
72  size_t h, a;
73  C_FLOAT64 xc, p, c, nt, New;
74  C_FLOAT64 fk[STORED];
75  bool ready;
76 
77  // initial point is first guess but we have to make sure that we
78  // are within the parameter domain
79  for (i = 0; i < mVariableSize; i++)
80  {
81  const COptItem & OptItem = *(*mpOptItem)[i];
82 
83  switch (OptItem.checkConstraint(OptItem.getStartValue()))
84  {
85  case - 1:
86  mCurrent[i] = *OptItem.getLowerBoundValue();
87  break;
88 
89  case 1:
90  mCurrent[i] = *OptItem.getUpperBoundValue();
91  break;
92 
93  case 0:
94  mCurrent[i] = OptItem.getStartValue();
95  break;
96  }
97 
98  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
99 
100  // The step must not contain any zeroes
101  mStep[i] = std::max(fabs(mCurrent[i]), 1.0);
102  }
103 
105 
106  if (!isnan(mEvaluationValue))
107  {
108  // and store that value
111 
112  // We found a new best value lets report it.
114  }
115 
116  // store this in all positions
117  for (a = 0; a < STORED; a++)
118  fk[a] = mCurrentValue;
119 
120  mAccepted = 0;
121 
122  // set the number of steps at one single temperature
123  nt = (C_FLOAT64)(5 * mVariableSize);
124 
125  if (nt < 100) nt = 100;
126 
127  // no temperature reductions yet
128  k = 0;
129 
130  do // number of internal cycles: max(5*mVariableSize, 100) * NS * mVariableSize
131  {
132  for (m = 0; m < nt && mContinue; m++) // step adjustments
133  {
134  for (j = 0; j < NS && mContinue; j++) // adjustment in all directions
135  {
136  for (h = 0; h < mVariableSize && mContinue; h++) // adjustment in one direction
137  {
138  // Calculate the step
139  xc = (2.0 * mpRandom->getRandomCC() - 1) * mStep[h];
140  New = mCurrent[h] + xc;
141 
142  // Set the new parameter value
143  (*(*mpSetCalculateVariable)[h])(New);
144 
145  // Check all parametric constraints
147  {
148  // Undo since not accepted
149  (*(*mpSetCalculateVariable)[h])(mCurrent[h]);
150  continue;
151  }
152 
153  // evaluate the function
154  evaluate();
155 
156  // Check all functional constraints
158  {
159  // Undo since not accepted
160  (*(*mpSetCalculateVariable)[h])(mCurrent[h]);
161  continue;
162  }
163 
164  // here we have a completely feasible point
165  // keep if energy is reduced
167  {
168  // only one value has changed...
169  mCurrent[h] = New;
171  i++; // a new point
172  mAccepted[h]++; // a new point in this coordinate
173 
175  {
176  // and store that value
179 
180  // We found a new best value lets report it.
182  }
183  }
184  else
185  {
186  // keep with probability p, if energy is increased
187  p = exp((mCurrentValue - mEvaluationValue) / (K * mTemperature));
188 
189  if (p > mpRandom->getRandomCO())
190  {
191  // only one value has changed...
192  mCurrent[h] = New;
194  i++; // a new point
195  mAccepted[h]++; // a new point in this coordinate
196  }
197  else
198  // Undo since not accepted
199  (*(*mpSetCalculateVariable)[h])(mCurrent[h]);
200  }
201  }
202  }
203 
204  // update the step sizes
205  for (a = 0; a < mVariableSize; a++)
206  {
207  c = (C_FLOAT64) mAccepted[a] / (C_FLOAT64) NS;
208 
209  if (c > 0.6)
210  mStep[a] *= 1 + 5 * (c - 0.6);
211  else if (c < 0.4)
212  mStep[a] /= 1 + 5 * (0.4 - c);
213 
214  mAccepted[a] = 0;
215  }
216  }
217 
218  // the equilibrium energy
219 
220  k++;
221 
222  // if this is the first cycle ignore the convergence tests
223  if (k == 1)
224  ready = false;
225  else
226  {
227  ready = true;
228 
229  // check termination criterion of not much change since last STORED times
230  for (a = 0; a < STORED; a++)
231  if (fabs(fk[a] - mCurrentValue) > mTolerance)
232  {
233  ready = false;
234  break;
235  }
236 
237  if (!ready)
238  {
239  for (a = 0; a < STORED - 1; a++)
240  fk[a] = fk[a + 1];
241 
242  fk[STORED - 1] = mCurrentValue;
243  }
244  // check the termination criterion of not much larger than last optimal
245  else if (fabs(mCurrentValue - mBestValue) > mTolerance)
246  ready = false;
247  }
248 
249  if (!ready)
250  {
251  i++;
252 
254 
255  for (a = 0; a < mVariableSize; a++)
256  (*(*mpSetCalculateVariable)[a])(mCurrent[a]);
257 
259  }
260 
261  // update the temperature
263 
264  if (mpCallBack)
265  mContinue &= mpCallBack->progressItem(mhTemperature);
266  }
267  while (!ready && mContinue);
268 
269  if (mpCallBack)
271 
272  return true;
273 }
274 
276 {
277  return true;
278 }
279 
281 {
282  // We do not need to check whether the parametric constraints are fulfilled
283  // since the parameters are created within the bounds.
284 
287 
288  // When we leave the either functional domain
289  // we set the objective value +Inf
291  mEvaluationValue = std::numeric_limits<C_FLOAT64>::infinity();
292 
293  return mEvaluationValue;
294 }
295 
297 {
298  cleanup();
299 
300  if (!COptMethod::initialize()) return false;
301 
302  mTemperature = * getValue("Start Temperature").pUDOUBLE;
303  mCoolingFactor = * getValue("Cooling Factor").pUDOUBLE;
304  mTolerance = * getValue("Tolerance").pUDOUBLE;
305  mpRandom =
306  CRandom::createGenerator(* (CRandom::Type *) getValue("Random Number Generator").pUINT,
307  * getValue("Seed").pUINT);
308 
309  if (mpCallBack)
310  mhTemperature =
311  mpCallBack->addItem("Current Temperature",
312  mTemperature);
313 
314  mBestValue = std::numeric_limits<C_FLOAT64>::infinity();
315  mContinue = true;
316 
317  mVariableSize = mpOptItem->size();
318 
322 
323  return true;
324 }
C_FLOAT64 mTolerance
Definition: COptMethodSA.h:109
virtual C_INT32 checkConstraint() const
Definition: COptItem.cpp:401
COptMethodSA(const COptMethodSA &src, const CCopasiContainer *pParent=NULL)
COptTask * mpParentTask
Definition: COptMethod.h:58
virtual bool cleanup()
virtual bool initialize()
Definition: COptMethod.cpp:189
CVector< C_FLOAT64 > mCurrent
Definition: COptMethodSA.h:139
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)
virtual bool initialize()
#define C_INT32
Definition: copasi.h:90
virtual bool progressItem(const size_t &handle)
size_t mVariableSize
Definition: COptMethodSA.h:119
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
C_FLOAT64 mCoolingFactor
Definition: COptMethodSA.h:104
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
C_FLOAT64 mTemperature
Definition: COptMethodSA.h:94
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
Definition: COptMethodSA.h:154
const CVector< C_FLOAT64 > & getSolutionVariables() const
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
virtual bool finishItem(const size_t &handle)
virtual bool checkFunctionalConstraints()
#define NS
CVector< C_FLOAT64 > mStep
Definition: COptMethodSA.h:149
size_t mhTemperature
Definition: COptMethodSA.h:99
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
C_FLOAT64 mCurrentValue
Definition: COptMethodSA.h:144
#define STORED
const C_FLOAT64 * getUpperBoundValue() const
Definition: COptItem.h:198
bool addParameter(const CCopasiParameter &parameter)
virtual bool checkParametricConstraints()
C_FLOAT64 mEvaluationValue
Definition: COptMethodSA.h:129
const C_FLOAT64 & evaluate()
const C_FLOAT64 & getCalculateValue() const
void initObjects()
virtual ~COptMethodSA()
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
CProcessReport * mpCallBack
virtual bool optimise()
C_FLOAT64 mBestValue
Definition: COptMethodSA.h:124
CRandom * mpRandom
Definition: COptMethodSA.h:114
#define K
#define max(a, b)
Definition: f2c.h:176