COPASI API  4.16.103
COptMethodSteepestDescent.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) 2005 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include "copasi.h"
16 
18 #include "COptProblem.h"
19 #include "COptItem.h"
20 #include "COptTask.h"
21 
22 #include "FminBrent.h"
23 
25 
27  COptMethod(CCopasiTask::optimization, CCopasiMethod::SteepestDescent, pParent),
28  mIterations(100),
29  mTolerance(1e-6),
30  mContinue(true),
31  mBestValue(std::numeric_limits< C_FLOAT64 >::infinity()),
32  mValue(0.0),
33  mVariableSize(0),
34  mIndividual(0),
35  mGradient(0),
36  mpDescent(new FDescentTemplate<COptMethodSteepestDescent>(this, &COptMethodSteepestDescent::descentLine)),
37  mCurrentIteration(0)
38 {
39  addParameter("Iteration Limit", CCopasiParameter::UINT, (unsigned C_INT32) 100);
40  addParameter("Tolerance", CCopasiParameter::DOUBLE, (C_FLOAT64) 1e-6);
41 }
42 
44  const CCopasiContainer * pParent): COptMethod(src, pParent),
45  mIterations(src.mIterations),
46  mTolerance(src.mTolerance),
47  mContinue(src.mContinue),
48  mBestValue(src.mBestValue),
49  mValue(src.mValue),
50  mVariableSize(src.mVariableSize),
51  mIndividual(src.mIndividual),
52  mGradient(src.mGradient),
53  mpDescent(new FDescentTemplate<COptMethodSteepestDescent>(this, &COptMethodSteepestDescent::descentLine)),
54  mCurrentIteration(src.mCurrentIteration)
55 {}
56 
58 {
60 
61  cleanup();
62 }
63 
65 {
66  if (!initialize()) return false;
67 
68  size_t i, k;
69  C_FLOAT64 tmp, x0, alpha, mn, mx, fmn, fmx;
70  bool calc_grad;
71 
72  // initial point is first guess but we have to make sure that we
73  // are within the parameter domain
74  for (i = 0; i < mVariableSize; i++)
75  {
76  const COptItem & OptItem = *(*mpOptItem)[i];
77 
78  switch (OptItem.checkConstraint(OptItem.getStartValue()))
79  {
80  case - 1:
81  mIndividual[i] = *OptItem.getLowerBoundValue();
82  break;
83 
84  case 1:
85  mIndividual[i] = *OptItem.getUpperBoundValue();
86  break;
87 
88  case 0:
89  mIndividual[i] = OptItem.getStartValue();
90  break;
91  }
92 
93  (*(*mpSetCalculateVariable)[i])(mIndividual[i]);
94  }
95 
96  fmx = mBestValue = evaluate();
97 
99 
100  // We found a new best value lets report it.
101  //if (mpReport) mpReport->printBody();
103 
104  bool SolutionFound = false;
105 
106  for (mCurrentIteration = 0; mCurrentIteration < mIterations && mContinue && !SolutionFound; mCurrentIteration++)
107  {
108  // calculate the direction of steepest descent
109  // by central finite differences
110  gradient();
111  // minimise the function in this direction
112  // find brackets for the minimisation
113  // mn = 0.0; md = 1.0;
114  // mnbrak(&mn, &md, &mx, &fmn, &fmd, &fmx, descent_line);
115  // make sure that no parameter will exceed its bounds
117 
118  for (i = 0; i < mVariableSize; i++)
119  {
120  if (fabs(mGradient[i]) > std::numeric_limits< C_FLOAT64 >::epsilon())
121  {
122  if (mGradient[i] > 0)
123  {
124  tmp = *(*mpOptItem)[i]->getUpperBoundValue();
125  }
126 
127  else
128  {
129  tmp = *(*mpOptItem)[i]->getLowerBoundValue();
130  }
131 
132  // calculate the size of the largest jump
133  tmp = (tmp - mIndividual[i]) / mGradient[i];
134 
135  // keep it if it is the smallest
136  if (tmp < x0) x0 = tmp;
137  }
138  else mGradient[i] = 0.0;
139  }
140 
141  if (x0 < mTolerance) x0 = mTolerance;
142 
143  // we will move at a rate of 1/10 this size
144  mn = mx = alpha = 0.0;
145 
146  for (k = 0, calc_grad = false; (k < 9) && !calc_grad && !SolutionFound; k++)
147  {
148  // set the minimum to the last successful step
149  mn = mx;
150  fmn = fmx;
151  // increment alpha
152  alpha += 0.1 * x0;
153  // set the maximum to it
154  mx = alpha;
155 
156  // take one step in that direction
157  fmx = descentLine(alpha);
158 
159  fmx = evaluate();
160 
161  // if this was an upward step find the minimum
162  if (fmx > fmn)
163  {
164  //md = mn + (mx-mn)/2;
165  //Brent(mn, md, mx, descent_line, &alpha, &tmp, 1e-6, 50);
166 
167  FminBrent(mn, mx, mpDescent, &alpha, &tmp, mTolerance, 5);
168 
169  // take one step in that direction
170  fmx = descentLine(alpha);
171 
172  // time to evaluate the new steepest descent direction
173  calc_grad = true;
174  }
175 
176  if (fabs(fmx - mBestValue) < mTolerance)
177  SolutionFound = true;
178  }
179 
180  for (i = 0; i < mVariableSize; i++)
181  mIndividual[i] = *(*mpOptItem)[i]->getObjectValue();
182 
183  if (fmx < mBestValue)
184  {
185  mBestValue = fmx;
186 
188 
189  // We found a new best value lets report it.
190  //if (mpReport) mpReport->printBody();
192  }
193  }
194 
195  return true;
196 }
197 
199 {
200  // pdelete all used variables
201  return true;
202 }
203 
205 {
206  cleanup();
207 
208  if (!COptMethod::initialize()) return false;
209 
210  mIterations = * getValue("Iteration Limit").pUINT;
211  mTolerance = * getValue("Tolerance").pDOUBLE;
212 
213  mContinue = true;
214  mVariableSize = mpOptItem->size();
217 
218  mBestValue = std::numeric_limits<C_FLOAT64>::infinity();
219 
220  return true;
221 }
222 
224 {
225  size_t i;
226 
227  C_FLOAT64 y;
228  C_FLOAT64 x;
229 
230  y = evaluate();
231 
232  for (i = 0; i < mVariableSize && mContinue; i++)
233  {
234  if ((x = *(*mpOptItem)[i]->getObjectValue()) != 0.0)
235  {
236  (*(*mpSetCalculateVariable)[i])(x * 1.001);
237  mGradient[i] = (y - evaluate()) / (x * 0.001);
238  }
239 
240  else
241  {
242  (*(*mpSetCalculateVariable)[i])(1e-7);
243  mGradient[i] = (y - evaluate()) / 1e-7;
244  }
245 
246  (*(*mpSetCalculateVariable)[i])(x);
247  }
248 }
249 
251 {
252  for (size_t i = 0; i < mVariableSize; i++)
253  (*(*mpSetCalculateVariable)[i])(mIndividual[i] + x * mGradient[i]);
254 
255  return evaluate();
256 }
257 
258 // evaluate the fitness of one individual
260 {
261  // evaluate the fitness
263 
265 
266  // when we leave the either the parameter or functional domain
267  // we penalize the objective value by forcing it to be larger
268  // than the best value recorded so far.
269  if (mValue < mBestValue &&
272  mValue = mBestValue + fabs(mBestValue - mValue);
273 
274  return mValue;
275 }
276 
278 {
280 }
virtual C_INT32 checkConstraint() const
Definition: COptItem.cpp:401
COptTask * mpParentTask
Definition: COptMethod.h:58
#define pdelete(p)
Definition: copasi.h:215
virtual bool initialize()
Definition: COptMethod.cpp:189
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
COptMethodSteepestDescent(const CCopasiContainer *pParent=NULL)
virtual bool calculate()
const std::vector< UpdateMethod * > * mpSetCalculateVariable
Definition: COptMethod.h:65
const Value & getValue() const
const C_FLOAT64 * getLowerBoundValue() const
Definition: COptItem.h:191
unsigned C_INT32 * pUINT
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
int FminBrent(double a, double b, FDescent *pF, double *min, double *fmin, double tol, int maxiter)
Definition: FminBrent.cpp:30
virtual bool checkFunctionalConstraints()
const C_FLOAT64 & getStartValue() const
Definition: COptItem.cpp:199
const std::vector< COptItem * > * mpOptItem
Definition: COptMethod.h:70
#define C_FLOAT64
Definition: copasi.h:92
const C_FLOAT64 * getUpperBoundValue() const
Definition: COptItem.h:198
bool addParameter(const CCopasiParameter &parameter)
C_FLOAT64 descentLine(const C_FLOAT64 &x)
virtual bool checkParametricConstraints()
const C_FLOAT64 & getCalculateValue() const
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
#define max(a, b)
Definition: f2c.h:176