ANIMA  4.0
animaNLOPTOptimizers.cxx
Go to the documentation of this file.
1 /*
2  * ITK Wrapper for the NLOPT library
3  *
4  * This code was graciously provided to us by the CRL at Boston. Please do NOT distribute
5  * it in any form or to anyone.
6  *
7  * Copyright (c) 2010-2011 Children's Hospital Boston.
8  * Benoit Scherrer, CRL (Computational Radiology Laboratory), Harvard Medical School
9  *
10  * This software is licensed by the copyright holder under the terms of the
11  * Open Software License version 3.0.
12  * http://www.opensource.org/licenses/osl-3.0.php
13  *
14  * Attribution Notice.
15  *
16  * This research was carried out in the Computational Radiology Laboratory of
17  * Children's Hospital, Boston and Harvard Medical School.
18  * http://www.crl.med.harvard.edu
19  * For more information contact: simon.warfield@childrens.harvard.edu
20  *
21  * This research work was made possible by Grant Number R01 RR021885 (Principal
22  * Investigator: Simon K. Warfield, Ph.D.) to Children's Hospital, Boston
23  * from the National Center for Research Resources (NCRR), a component of the
24  * National Institutes of Health (NIH).
25 */
26 #include "animaNLOPTOptimizers.h"
27 
28 #include <math.h>
29 #include <vnl/vnl_math.h>
30 
31 #include <itkTimeProbe.h>
32 
33 using namespace std;
34 
35 namespace anima
36 {
37  /**********************************************************************************************/
45  NLOPTOptimizers::NLOPTOptimizers()
46  {
47  m_VerboseLevel = 0;
48  m_Maximize = false;
49  m_ErrorCode = NLOPT_FAILURE;
50  m_Algorithm = NLOPT_LN_BOBYQA;
51  m_LocalOptimizer = NLOPT_LN_BOBYQA;
52 
53  m_StopValSet = false;
54  m_StopVal = 0;
55  m_FTolRel = -1;
56  m_FTolAbs = -1;
57  m_XTolRel = -1;
58  m_XTolAbs = -1;
59  m_MaxEval = -1;
60  m_MaxTime = -1;
61  m_PopulationSize = -1;
62  m_VectorStorageSize = -1;
63  m_CurrentCost = 0;
64  }
65 
66  /**********************************************************************************************/
74  NLOPTOptimizers::~NLOPTOptimizers()
75  {
76  }
77 
78  /**********************************************************************************************/
95  double NLOPTOptimizers::NloptFunctionWrapper(unsigned n, const double *x, double *grad, void *data)
96  {
97  NLOPTOptimizers *optimizer = static_cast<NLOPTOptimizers *>(data);
98 
99  //-----------------------------------------
100  // Copy the nlopt position to a itk position
101  // (takes into account the itk scale)
102  //-----------------------------------------
103  NLOPTOptimizers::ParametersType itkCurrentPosition(n);
104  for ( unsigned int i=0; i<n ; i++ )
105  itkCurrentPosition[i] = x[i]/optimizer->GetScales()[i];
106 
107  double f = (optimizer->GetCostFunction()->GetValue(itkCurrentPosition));
108 
109  //-----------------------------------------
110  // If needed, also compute the derivative and copy
111  // back in *grad
112  //-----------------------------------------
113  if ( grad!=NULL )
114  {
115  DerivativeType derivative;
116  optimizer->GetCostFunction()->GetDerivative (itkCurrentPosition, derivative);
117  for ( unsigned int i=0; i<n ; i++ )
118  grad[i] = derivative[i];
119  }
120 
121  return f;
122  }
123 
124  /**********************************************************************************************/
134  void NLOPTOptimizers::StartOptimization()
135  {
136  if( m_CostFunction.IsNull() )
137  throw itk::ExceptionObject(__FILE__, __LINE__, "Error. The cost function has not been set", "NLOPTOptimizers");
138 
139  //---------------------------------------------
140  // n = number of parameters
141  //---------------------------------------------
142  unsigned int n = m_CostFunction->GetNumberOfParameters();
143 
144  //---------------------------------------------
145  // Set x = initial position
146  // Take into account the ITK scales
147  //---------------------------------------------
148  if ( this->GetInitialPosition().GetSize() != n )
149  throw itk::ExceptionObject(__FILE__, __LINE__, "Invalid initial position parameter. Its size should be equal to the number of parameters", "NLOPTOptimizers");
150  m_CurrentPosition.set_size(n);
151 
152  double *x = new double[n];
153  const double *in_x = GetInitialPosition().data_block();
154  for ( unsigned int i=0; i<n ; i++ )
155  {
156  x[i] = in_x[i] * GetScales()[i];
157  m_CurrentPosition[i] = in_x[i];
158  }
159 
160  //---------------------------------------------
161  // lb/ub = lower bound/upper bound
162  // Take into account the ITK scales
163  //---------------------------------------------
164  double *lb=NULL;
165  double *ub=NULL;
166  if ( m_LowerBoundParameters.GetSize()!=0 )
167  {
168  if ( m_LowerBoundParameters.GetSize()!=n )
169  throw itk::ExceptionObject(__FILE__, __LINE__, "Invalid lower bound parameter. Its size should be equal to the number of parameters", "NLOPTOptimizers");
170 
171  lb = new double[n];
172  for ( unsigned int i=0; i<n; i++ )
173  lb[i] = m_LowerBoundParameters[i]*GetScales()[i];
174  }
175  if ( m_UpperBoundParameters.GetSize()!=0 )
176  {
177  if ( m_UpperBoundParameters.GetSize()!=n )
178  throw itk::ExceptionObject(__FILE__, __LINE__, "Invalid upper bound parameter. Its size should be equal to the number of parameters", "NLOPTOptimizers");
179 
180  ub = new double[n];
181  for ( unsigned int i=0; i<n; i++ )
182  ub[i] = m_UpperBoundParameters[i]*GetScales()[i];
183  }
184 
185  //---------------------------------------------
186  // Creates the NLOPT structure and fill it
187  //---------------------------------------------
188  m_NloptOptions = nlopt_create((::nlopt_algorithm)(int)m_Algorithm, n);
189  nlopt_set_lower_bounds(m_NloptOptions, lb);
190  nlopt_set_upper_bounds(m_NloptOptions, ub);
191  if ( m_Maximize )
192  nlopt_set_max_objective(m_NloptOptions, (nlopt_func)this->NloptFunctionWrapper, (void *)this);
193  else
194  nlopt_set_min_objective(m_NloptOptions, (nlopt_func)this->NloptFunctionWrapper, (void *)this);
195 
196  if ( m_StopValSet ) nlopt_set_stopval(m_NloptOptions, m_StopVal);
197 
198  nlopt_set_ftol_rel(m_NloptOptions, m_FTolRel);
199  nlopt_set_ftol_abs(m_NloptOptions, m_FTolAbs);
200  nlopt_set_xtol_rel(m_NloptOptions, m_XTolRel);
201  nlopt_set_xtol_abs1(m_NloptOptions, m_XTolAbs);
202  nlopt_set_maxtime(m_NloptOptions, m_MaxTime);
203  nlopt_set_vector_storage(m_NloptOptions, m_VectorStorageSize);
204  nlopt_set_maxeval(m_NloptOptions, m_MaxEval);
205  nlopt_set_population(m_NloptOptions, m_PopulationSize);
206 
207  for (unsigned int i = 0;i < m_InequalityConstraints.size();++i)
208  nlopt_add_inequality_constraint(m_NloptOptions, (nlopt_func)ConstraintsFunctionType::GetConstraintValue, m_InequalityConstraints[i]->GetAdditionalData(), m_InequalityConstraints[i]->GetTolerance());
209 
210  for (unsigned int i = 0;i < m_EqualityConstraints.size();++i)
211  nlopt_add_equality_constraint(m_NloptOptions, (nlopt_func)ConstraintsFunctionType::GetConstraintValue, m_EqualityConstraints[i]->GetAdditionalData(), m_EqualityConstraints[i]->GetTolerance());
212 
213  //----------------------------------------
214  // Setup local optimizer for algorithms
215  // using sequences of local optimizations
216  //----------------------------------------
217  m_NloptLocalOptions = nlopt_create((::nlopt_algorithm)(int)m_LocalOptimizer, n);
218 
219  if ( m_StopValSet ) nlopt_set_stopval(m_NloptLocalOptions, m_StopVal);
220 
221  nlopt_set_ftol_rel(m_NloptLocalOptions, m_FTolRel);
222  nlopt_set_ftol_abs(m_NloptLocalOptions, m_FTolAbs);
223  nlopt_set_xtol_rel(m_NloptLocalOptions, m_XTolRel);
224  nlopt_set_xtol_abs1(m_NloptLocalOptions, m_XTolAbs);
225  nlopt_set_maxtime(m_NloptLocalOptions, m_MaxTime);
226  nlopt_set_vector_storage(m_NloptLocalOptions, m_VectorStorageSize);
227  nlopt_set_maxeval(m_NloptLocalOptions, m_MaxEval);
228 
229  //--------------------------------------
230  // Plug local optimizer into global one
231  //--------------------------------------
232  nlopt_set_local_optimizer(m_NloptOptions, m_NloptLocalOptions);
233 
234  this->InvokeEvent( itk::StartEvent() );
235 
236  //----------------------------------------
237  // Run the NLOPT optimizer !
238  //----------------------------------------
239  double valf;
240  m_ErrorCode = static_cast<nlopt_result>((int)nlopt_optimize (m_NloptOptions, x, &valf));
241  SetCurrentCost(valf);
242 
243  //----------------------------------------
244  // Converts back using the scales
245  //----------------------------------------
247  for ( unsigned int i=0; i<n ; i++ )
248  p[i] = x[i]/GetScales()[i];
249  this->SetCurrentPosition(p);
250 
251  //----------------------------------------
252  // Free memory
253  //----------------------------------------
254  delete[] x;
255  if ( ub!=NULL ) delete[] ub;
256  if ( lb!=NULL ) delete[] lb;
257 
258  nlopt_destroy(m_NloptOptions);
259  nlopt_destroy(m_NloptLocalOptions);
260 
261  this->InvokeEvent( itk::EndEvent() );
262  }
263 
264  /**********************************************************************************************/
275  bool NLOPTOptimizers::isSuccessful() const
276  {
277  if (m_ErrorCode < 0
278  && m_ErrorCode != NLOPT_ROUNDOFF_LIMITED
279  && m_ErrorCode != NLOPT_FORCED_STOP)
280  return false;
281  else
282  return true;
283  }
284 
285  /**********************************************************************************************/
295  std::string NLOPTOptimizers::GetErrorCodeDescription() const
296  {
297  switch ( m_ErrorCode )
298  {
299  case NLOPT_SUCCESS:
300  return "Success";
301 
302  case NLOPT_STOPVAL_REACHED:
303  return "Optimization stopped because 'StopVal' was reached.";
304  case NLOPT_FTOL_REACHED:
305  return "Optimization stopped because 'FTolRel' or 'FTolAbs' was reached.";
306  case NLOPT_XTOL_REACHED:
307  return "Optimization stopped because 'XTolRel' or 'XTolAbs' was reached.";
308  case NLOPT_MAXEVAL_REACHED:
309  return "Optimization stopped because 'MaxEval' was reached.";
310  case NLOPT_MAXTIME_REACHED:
311  return "Optimization stopped because 'MaxTime' was reached.";
312 
313 
314  case NLOPT_FAILURE:
315  return "Generic failure code.";
316  case NLOPT_INVALID_ARGS:
317  return "Invalid arguments (e.g. lower bounds are bigger than upper bounds, an unknown algorithm was specified, etcetera).";
318  case NLOPT_OUT_OF_MEMORY:
319  return "Out of memory.";
320  case NLOPT_ROUNDOFF_LIMITED:
321  return "Halted because roundoff errors limited progress.";
322  case NLOPT_FORCED_STOP:
323  return "Halted because of a forced termination.";
324 
325  default:
326  break;
327  }
328  return "Unkown error code";
329  }
330 
331  /**********************************************************************************************/
346  void NLOPTOptimizers::AddInequalityConstraint(ConstraintsFunctionType *constraint)
347  {
348  m_InequalityConstraints.push_back(constraint);
349  }
350 
351  /**********************************************************************************************/
359  void NLOPTOptimizers::ClearInequalityConstraints()
360  {
361  m_InequalityConstraints.clear();
362  }
363 
364  /**********************************************************************************************/
381  void NLOPTOptimizers::AddEqualityConstraint(ConstraintsFunctionType *constraint)
382  {
383  m_EqualityConstraints.push_back(constraint);
384  }
385 
386  /**********************************************************************************************/
394  void NLOPTOptimizers::ClearEqualityConstraints()
395  {
396  m_EqualityConstraints.clear();
397  }
398 
399  /**********************************************************************************************/
409  void NLOPTOptimizers::SetLowerBoundParameters( const ParametersType& p )
410  {
411  m_LowerBoundParameters.SetSize(p.GetSize());
412  for ( unsigned int i=0; i<p.GetSize(); i++ )
413  m_LowerBoundParameters[i] = p[i];
414  }
415 
416  /**********************************************************************************************/
426  void NLOPTOptimizers::SetUpperBoundParameters( const ParametersType& p )
427  {
428  m_UpperBoundParameters.SetSize(p.GetSize());
429  for ( unsigned int i=0; i<p.GetSize(); i++ )
430  m_UpperBoundParameters[i] = p[i];
431  }
432 
436  void NLOPTOptimizers::PrintSelf(std::ostream& os, itk::Indent indent) const
437  {
438  Superclass::PrintSelf(os,indent);
439 
440  int major, minor, bugfix;
441  nlopt_version(&major, &minor, &bugfix);
442  os << indent << "NLOPT v"<<major<<"."<<minor<<"."<<bugfix<<std::endl;
443 
444  os << indent << "Algorithm: " << ::nlopt_algorithm_name((::nlopt_algorithm)(int)m_Algorithm) << std::endl;
445  os << indent << "Maximize On/Off " << m_Maximize << std::endl;
446 
447  os << indent << "StopVal: "<<m_StopVal<<std::endl;
448  os << indent << "FTolRel: "<<m_FTolRel<<std::endl;
449  os << indent << "FTolAbs: "<<m_FTolAbs<<std::endl;
450  os << indent << "XTolRel: "<<m_XTolRel<<std::endl;
451  os << indent << "MaxEval: "<<m_MaxEval<<std::endl;
452  os << indent << "MaxTime: "<<m_MaxTime<<std::endl;
453  }
454 
455 } // end of namespace anima
Implements an ITK wrapper for the NLOPT library.
STL namespace.
MeasureType GetValue() const
SingleValuedNonLinearOptimizer::ParametersType ParametersType