COPASI API  4.16.103
COptMethodTruncatedNewton.cpp
Go to the documentation of this file.
1 // Begin CVS Header
2 // $Source: /Volumes/Home/Users/shoops/cvs/copasi_dev/copasi/optimization/COptMethodTruncatedNewton.cpp,v $
3 // $Revision: 1.11 $
4 // $Name: $
5 // $Author: shoops $
6 // $Date: 2012/06/20 21:16:37 $
7 // End CVS Header
8 
9 // Copyright (C) 2012 - 2010 by Pedro Mendes, Virginia Tech Intellectual
10 // Properties, Inc., University of Heidelberg, and The University
11 // of Manchester.
12 // All rights reserved.
13 
14 // Copyright (C) 2008 by Pedro Mendes, Virginia Tech Intellectual
15 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
16 // and The University of Manchester.
17 // All rights reserved.
18 
19 // Copyright (C) 2001 - 2007 by Pedro Mendes, Virginia Tech Intellectual
20 // Properties, Inc. and EML Research, gGmbH.
21 // All rights reserved.
22 
23 #include "copasi.h"
24 
26 #include "COptProblem.h"
27 #include "COptItem.h"
28 #include "COptTask.h"
29 
32 
34  COptMethod(CCopasiTask::optimization, CCopasiMethod::TruncatedNewton, pParent),
36  mpCTruncatedNewton(new CTruncatedNewton())
37 {initObjects();}
38 
40  const CCopasiContainer * pParent):
41  COptMethod(src, pParent),
43  mpCTruncatedNewton(new CTruncatedNewton())
44 {initObjects();}
45 
47 {
50  cleanup();
51 }
52 
54 {
56 }
57 
59 {
60  if (!initialize()) return false;
61 
62  C_FLOAT64 fest;
63  C_INT lw, ierror = 0;
64  lw = 14 * mVariableSize;
65 
69  CVector< C_FLOAT64 > dwork(lw);
70 
73 
74  // initial point is the first guess but we have to make sure that
75  // we are within the parameter domain
76  C_INT i, repeat;
77 
78  for (i = 0; i < mVariableSize; i++)
79  {
80  const COptItem & OptItem = *(*mpOptItem)[i];
81 
82  // :TODO: In COPASI the bounds are not necessarry fixed.
83  // Since evaluate checks for boundaries and constraints this is not
84  // needed. The question remaining is how does tnbc_ handle unconstraint problems?
85  // low[i] = *OptItem.getLowerBoundValue();
86  // up[i] = *OptItem.getUpperBoundValue();
87 
88  mCurrent[i] = OptItem.getStartValue();
89 
90  switch (OptItem.checkConstraint(mCurrent[i]))
91  {
92  case - 1:
93  mCurrent[i] = *OptItem.getLowerBoundValue();
94  break;
95 
96  case 1:
97  mCurrent[i] = *OptItem.getUpperBoundValue();
98  break;
99 
100  case 0:
101  break;
102  }
103 
104  // set the value
105  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
106  }
107 
108  // Report the first value as the current best
109  mBestValue = evaluate();
110  mBest = mCurrent;
112 
113  repeat = 0;
114 
115  while (repeat < 10 && mContinue)
116  {
117  repeat++;
118 
119  // estimate minimum is 1/10 initial function value
120  fest = (1 - pow(0.9, (C_FLOAT64) repeat)) * mEvaluationValue;
121  ierror = 0;
122 
123  // minimise
124  try
125  {
126  mpCTruncatedNewton->tnbc_(&ierror, &mVariableSize, mCurrent.array(), &fest, mGradient.array(), dwork.array(),
127  &lw, mpTruncatedNewton, low.array(), up.array(), iPivot.array());
128  mEvaluationValue = fest;
129  }
130 
131  // This signals that the user opted to interupt
132  catch (bool)
133  {
134  break;
135  }
136 
137  if (ierror < 0)
138  fatalError(); // Invalid parameter values.
139 
140  // The way the method is currently implemented may lead to parameters just outside the boundaries.
141  // We need to check whether the current value is within the boundaries or whether the corrected
142  // leads to an improved solution.
143 
144  bool withinBounds = true;
145 
146  for (i = 0; i < mVariableSize; i++)
147  {
148  const COptItem & OptItem = *(*mpOptItem)[i];
149 
150  //force it to be within the bounds
151  switch (OptItem.checkConstraint(mCurrent[i]))
152  {
153  case - 1:
154  withinBounds = false;
155  mCurrent[i] = *OptItem.getLowerBoundValue();
156  break;
157 
158  case 1:
159  withinBounds = false;
160  mCurrent[i] = *OptItem.getUpperBoundValue();
161  break;
162 
163  case 0:
164  break;
165  }
166 
167  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
168  }
169 
170  evaluate();
171 
172  // Is the corrected value better than solution?
174  {
175  // We found a new best value lets report it.
176  // and store that value
177  mBest = mCurrent;
179 
181 
182  // We found a new best value lets report it.
184  }
185 
186  // We found a solution
187  if (withinBounds)
188  break;
189 
190  // Choosing another starting point will be left to the user
191 #ifdef XXXX
192 
193  // Try another starting point
194  for (i = 0; i < mVariableSize; i++)
195  {
196  mCurrent[i] *= 1.2;
197  const COptItem & OptItem = *(*mpOptItem)[i];
198 
199  //force it to be within the bounds
200  switch (OptItem.checkConstraint(mCurrent[i]))
201  {
202  case - 1:
203  mCurrent[i] = *OptItem.getLowerBoundValue();
204  break;
205 
206  case 1:
207  mCurrent[i] = *OptItem.getUpperBoundValue();
208  break;
209 
210  case 0:
211  break;
212  }
213 
214  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
215  }
216 
217  evaluate();
218 
219  // Check whether we improved
221  {
222  // We found a new best value lets report it.
223  // and store that value
224  mBest = mCurrent;
226 
228 
229  // We found a new best value lets report it.
231  }
232 
233 #endif // XXXX
234  }
235 
236  return true;
237 }
238 
240 {
241  cleanup();
242 
243  if (!COptMethod::initialize()) return false;
244 
245  mVariableSize = (C_INT) mpOptItem->size();
248 
249  mContinue = true;
250  mBestValue = std::numeric_limits<C_FLOAT64>::infinity();
252 
253  return true;
254 }
255 
257 {
258  return true;
259 }
260 
261 // callback function, evaluate the value of the objective function and its gradient
262 //(by finite differences), translated by f2c, edited by Pedro and then modified for COPASI by Joseph
264 {
265  C_INT i;
266 
267  // set the parameter values
268  for (i = 0; i < *n; i++)
269  (*(*mpSetCalculateVariable)[i])(x[i]);
270 
271  //carry out the function evaluation
272  *f = evaluate();
273 
274  // Check whether we improved
276  {
277  // We found a new best value lets report it.
278  // and store that value
279  for (i = 0; i < *n; i++)
280  mBest[i] = x[i];
281 
284 
285  // We found a new best value lets report it.
287  }
288 
289  // Calculate the gradient
290  for (i = 0; i < *n && mContinue; i++)
291  {
292  if (x[i] != 0.0)
293  {
294  (*(*mpSetCalculateVariable)[i])(x[i] * 1.001);
295  g[i] = (evaluate() - *f) / (x[i] * 0.001);
296  }
297 
298  else
299  {
300  (*(*mpSetCalculateVariable)[i])(1e-7);
301  g[i] = (evaluate() - *f) / 1e-7;
302  }
303 
304  (*(*mpSetCalculateVariable)[i])(x[i]);
305  }
306 
307  if (!mContinue)
308  throw bool(mContinue);
309 
310  return 0;
311 }
312 
314 {
315  // We do not need to check whether the parametric constraints are fulfilled
316  // since the parameters are created within the bounds.
319 
320  // when we leave the either the parameter or functional domain
321  // we penalize the objective value by forcing it to be larger
322  // than the best value recorded so far.
327 
328  return mEvaluationValue;
329 }
#define C_INT
Definition: copasi.h:115
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
#define fatalError()
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 calculate()
const std::vector< UpdateMethod * > * mpSetCalculateVariable
Definition: COptMethod.h:65
const C_FLOAT64 * getLowerBoundValue() const
Definition: COptItem.h:191
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
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
CType * array()
Definition: CVector.h:139
const C_FLOAT64 * getUpperBoundValue() const
Definition: COptItem.h:198
C_INT sFun(C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
virtual bool checkParametricConstraints()
COptMethodTruncatedNewton(const COptMethodTruncatedNewton &src, const CCopasiContainer *pParent=NULL)
const C_FLOAT64 & getCalculateValue() const
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
int tnbc_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *, C_INT *)
#define max(a, b)
Definition: f2c.h:176