Orfeo Toolbox  4.0
itkLBFGSBOptimizer.cxx
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef __itkLBFGSBOptimizer_hxx
19 #define __itkLBFGSBOptimizer_hxx
20 
21 #include "itkLBFGSBOptimizer.h"
22 #include "vnl/algo/vnl_lbfgsb.h"
23 
24 extern "C" {
25 extern double v3p_netlib_dpmeps_();
26 }
27 
28 namespace itk
29 {
37  public vnl_lbfgsb
38 {
39 public:
41  typedef vnl_lbfgsb Superclass;
42 
44  LBFGSBOptimizerHelper(vnl_cost_function & f,
45  LBFGSBOptimizer * const itkObj);
46 
48  virtual bool report_iter();
49 
50 private:
52 };
53 
59  m_Trace(false),
60  m_OptimizerInitialized(false),
61  m_CostFunctionConvergenceFactor(1e+7),
62  m_ProjectedGradientTolerance(1e-5),
63  m_MaximumNumberOfIterations(500),
64  m_MaximumNumberOfEvaluations(500),
65  m_MaximumNumberOfCorrections(5),
66  m_CurrentIteration(0),
67  m_InfinityNormOfProjectedGradient(0.0),
68  m_VnlOptimizer(0)
69 {
70  m_LowerBound = InternalBoundValueType(0);
71  m_UpperBound = InternalBoundValueType(0);
72  m_BoundSelection = InternalBoundSelectionType(0);
73 }
74 
80 {
81  delete m_VnlOptimizer;
82 }
83 
87 void
89 ::PrintSelf(std::ostream & os, Indent indent) const
90 {
91  Superclass::PrintSelf(os, indent);
92  os << indent << "Trace: ";
93  if ( m_Trace )
94  {
95  os << "On";
96  }
97  else
98  {
99  os << "Off";
100  }
101  os << std::endl;
102 
103  os << indent << "LowerBound: " << m_LowerBound << std::endl;
104  os << indent << "UpperBound: " << m_UpperBound << std::endl;
105  os << indent << "BoundSelection: " << m_BoundSelection << std::endl;
106 
107  os << indent << "CostFunctionConvergenceFactor: "
108  << m_CostFunctionConvergenceFactor << std::endl;
109 
110  os << indent << "ProjectedGradientTolerance: "
111  << m_ProjectedGradientTolerance << std::endl;
112 
113  os << indent << "MaximumNumberOfIterations: "
114  << m_MaximumNumberOfIterations << std::endl;
115 
116  os << indent << "MaximumNumberOfEvaluations: "
117  << m_MaximumNumberOfEvaluations << std::endl;
118 
119  os << indent << "MaximumNumberOfCorrections: "
120  << m_MaximumNumberOfCorrections << std::endl;
121 
122  os << indent << "CurrentIteration: "
123  << m_CurrentIteration << std::endl;
124 
125  os << indent << "Value: "
126  << this->GetValue() << std::endl;
127 
128  os << indent << "InfinityNormOfProjectedGradient: "
129  << m_InfinityNormOfProjectedGradient << std::endl;
130 
131  if ( this->m_VnlOptimizer )
132  {
133  os << indent << "Vnl LBFGSB Failure Code: "
134  << this->m_VnlOptimizer->get_failure_code() << std::endl;
135  }
136 }
137 
141 void
143 ::SetTrace(bool flag)
144 {
145  if ( flag == m_Trace )
146  {
147  return;
148  }
149 
150  m_Trace = flag;
151  if ( m_OptimizerInitialized )
152  {
153  m_VnlOptimizer->set_trace(m_Trace);
154  }
155 
156  this->Modified();
157 }
158 
162 void
165  const BoundValueType & value)
166 {
167  this->m_LowerBound = value;
168  if ( m_OptimizerInitialized )
169  {
170  m_VnlOptimizer->set_lower_bound(m_LowerBound);
171  }
172  this->Modified();
173 }
174 
178 void
181  const BoundValueType & value)
182 {
183  this->m_UpperBound = value;
184  if ( m_OptimizerInitialized )
185  {
186  m_VnlOptimizer->set_upper_bound(m_UpperBound);
187  }
188  this->Modified();
189 }
190 
196 ::GetValue() const
197 {
198  return this->GetCachedValue();
199 }
200 
204 void
207  const BoundSelectionType & value)
208 {
209  m_BoundSelection = value;
210  if ( m_OptimizerInitialized )
211  {
212  m_VnlOptimizer->set_bound_selection(m_BoundSelection);
213  }
214  this->Modified();
215 }
216 
223 void
226 {
227  if ( value < 0.0 )
228  {
229  itkExceptionMacro("Value " << value
230  << " is too small for SetCostFunctionConvergenceFactor()"
231  << "a typical range would be from 0.0 to 1e+12");
232  }
233  m_CostFunctionConvergenceFactor = value;
234  if ( m_OptimizerInitialized )
235  {
236  m_VnlOptimizer->set_cost_function_convergence_factor(
237  m_CostFunctionConvergenceFactor);
238  }
239  this->Modified();
240 }
241 
246 void
249 {
250  m_ProjectedGradientTolerance = value;
251  if ( m_OptimizerInitialized )
252  {
253  m_VnlOptimizer->set_projected_gradient_tolerance(
254  m_ProjectedGradientTolerance);
255  }
256  this->Modified();
257 }
258 
260 void
263 {
264  m_MaximumNumberOfIterations = value;
265  this->Modified();
266 }
267 
269 void
272 {
273  m_MaximumNumberOfEvaluations = value;
274  if ( m_OptimizerInitialized )
275  {
276  m_VnlOptimizer->set_max_function_evals(m_MaximumNumberOfEvaluations);
277  }
278  this->Modified();
279 }
280 
282 void
285 {
286  m_MaximumNumberOfCorrections = value;
287  if ( m_OptimizerInitialized )
288  {
289  m_VnlOptimizer->set_max_variable_metric_corrections(
290  m_MaximumNumberOfCorrections);
291  }
292  this->Modified();
293 }
294 
298 void
300 {
301  m_CostFunction = costFunction;
302 
303  const unsigned int numberOfParameters =
304  costFunction->GetNumberOfParameters();
305 
306  CostFunctionAdaptorType *adaptor =
307  new CostFunctionAdaptorType(numberOfParameters);
308 
309  adaptor->SetCostFunction(costFunction);
310 
312  {
313  delete m_VnlOptimizer;
314  }
315 
316  this->SetCostFunctionAdaptor(adaptor);
317 
318  m_VnlOptimizer = new InternalOptimizerType(*adaptor, this);
319 
320  // set the optimizer parameters
321  m_VnlOptimizer->set_lower_bound(m_LowerBound);
322  m_VnlOptimizer->set_upper_bound(m_UpperBound);
323  m_VnlOptimizer->set_bound_selection(m_BoundSelection);
324  m_VnlOptimizer->set_cost_function_convergence_factor(
326  m_VnlOptimizer->set_projected_gradient_tolerance(
328  m_VnlOptimizer->set_max_function_evals(m_MaximumNumberOfEvaluations);
329  m_VnlOptimizer->set_max_variable_metric_corrections(
331 
332  m_OptimizerInitialized = true;
333 
334  this->Modified();
335 }
336 
340 void
343 {
344  // Check if all the bounds parameters are the same size as the initial
345  // parameters.
346  unsigned int numberOfParameters = m_CostFunction->GetNumberOfParameters();
347 
348  if ( this->GetInitialPosition().Size() < numberOfParameters )
349  {
350  itkExceptionMacro(<< "InitialPosition array does not have sufficient number of elements");
351  }
352 
353  if ( m_LowerBound.size() < numberOfParameters )
354  {
355  itkExceptionMacro(<< "LowerBound array does not have sufficient number of elements");
356  }
357 
358  if ( m_UpperBound.size() < numberOfParameters )
359  {
360  itkExceptionMacro(<< "UppperBound array does not have sufficient number of elements");
361  }
362 
363  if ( m_BoundSelection.size() < numberOfParameters )
364  {
365  itkExceptionMacro(<< "BoundSelection array does not have sufficient number of elements");
366  }
367 
368  if ( this->GetMaximize() )
369  {
370  this->GetNonConstCostFunctionAdaptor()->NegateCostFunctionOn();
371  }
372  if(this->m_CostFunctionConvergenceFactor == 0.0 && this->m_ProjectedGradientTolerance == 0.0)
373  {
374  itkExceptionMacro("LBFGSB Optimizer cannot function if both CostFunctionConvergenceFactor"
375  " and ProjectedGradienctTolerance are zero.");
376  }
377  this->SetCurrentPosition( this->GetInitialPosition() );
378 
379  ParametersType parameters( this->GetInitialPosition() );
380 
381  this->InvokeEvent( StartEvent() );
382 
383  // vnl optimizers return the solution by reference
384  // in the variable provided as initial position
385  m_VnlOptimizer->minimize(parameters);
386 
387  if ( parameters.size() != this->GetInitialPosition().Size() )
388  {
389  // set current position to initial position and throw an exception
390  this->SetCurrentPosition( this->GetInitialPosition() );
391  itkExceptionMacro(<< "Error occurred in optimization");
392  }
393 
394  this->SetCurrentPosition(parameters);
395 
396  this->InvokeEvent( EndEvent() );
397 }
398 
399 /*-------------------------------------------------------------------------
400  * helper class
401  *-------------------------------------------------------------------------
402  */
403 
406 ::LBFGSBOptimizerHelper(vnl_cost_function & f,
407  LBFGSBOptimizer * const itkObj):
408  vnl_lbfgsb(f),
409  m_ItkObj(itkObj)
410 {}
411 
413 bool
416 {
417  Superclass::report_iter();
418 
419  m_ItkObj->m_InfinityNormOfProjectedGradient = this->get_inf_norm_projected_gradient();
420  m_ItkObj->InvokeEvent( IterationEvent() );
421  m_ItkObj->m_CurrentIteration = this->num_iterations_;
422 
423  // Return true to terminate the optimization loop.
424  if ( this->num_iterations_ > m_ItkObj->m_MaximumNumberOfIterations )
425  {
426  return true;
427  }
428  else
429  {
430  return false;
431  }
432 }
433 
434 const std::string
436 {
437  std::ostringstream stopConditionDescription;
439  {
440  stopConditionDescription << "Too many iterations. Iterations = "
441  << this->m_CurrentIteration << std::endl;
442  }
443  if ( m_VnlOptimizer )
444  {
445  stopConditionDescription << this->GetNameOfClass() << ": ";
446  switch ( m_VnlOptimizer->get_failure_code() )
447  {
448  case vnl_nonlinear_minimizer::ERROR_FAILURE:
449  stopConditionDescription << "Failure";
450  break;
451  case vnl_nonlinear_minimizer::ERROR_DODGY_INPUT:
452  stopConditionDescription << "Dodgy input";
453  break;
454  case vnl_nonlinear_minimizer::CONVERGED_FTOL:
455  stopConditionDescription << "Function tolerance reached after "
457  << " iterations. "
458  << "The relative reduction of the cost function <= "
460  << " = CostFunctionConvergenceFactor ("
462  << ") * machine precision ("
463  << v3p_netlib_dpmeps_()
464  << ").";
465 
466  break;
467  case vnl_nonlinear_minimizer::CONVERGED_XTOL:
468  stopConditionDescription << "Solution tolerance reached";
469  break;
470  case vnl_nonlinear_minimizer::CONVERGED_XFTOL:
471  stopConditionDescription << "Solution and Function tolerance both reached";
472  break;
473  case vnl_nonlinear_minimizer::CONVERGED_GTOL:
474  stopConditionDescription << "Gradient tolerance reached. "
475  << "Projected gradient tolerance is "
477  break;
478  case vnl_nonlinear_minimizer::FAILED_TOO_MANY_ITERATIONS:
479  stopConditionDescription << "Too many iterations. Iterations = "
481  break;
482  case vnl_nonlinear_minimizer::FAILED_FTOL_TOO_SMALL:
483  stopConditionDescription << "Function tolerance too small";
484  break;
485  case vnl_nonlinear_minimizer::FAILED_XTOL_TOO_SMALL:
486  stopConditionDescription << "Solution tolerance too small";
487  break;
488  case vnl_nonlinear_minimizer::FAILED_GTOL_TOO_SMALL:
489  stopConditionDescription << "Gradient tolerance too small";
490  break;
491  case vnl_nonlinear_minimizer::FAILED_USER_REQUEST:
492  break;
493  }
494  }
495  return stopConditionDescription.str();
496 }
497 } // end namespace itk
498 
499 #endif

Generated at Sat Mar 8 2014 15:06:52 for Orfeo Toolbox with doxygen 1.8.3.1