]> Creatis software - clitk.git/blob - registration/clitkLBFGSBOptimizer.cxx
Moved from repository clitk to clitk.private/tests_dav
[clitk.git] / registration / clitkLBFGSBOptimizer.cxx
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to: 
5   - University of LYON              http://www.universite-lyon.fr/
6   - Léon Bérard cancer center       http://www.centreleonberard.fr
7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
8
9   This software is distributed WITHOUT ANY WARRANTY; without even
10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11   PURPOSE.  See the copyright notices for more information.
12
13   It is distributed under dual licence
14
15   - BSD        See included LICENSE.txt file
16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
18 #ifndef __clitkLBFGSBOptimizer_cxx
19 #define __clitkLBFGSBOptimizer_cxx
20 #include "clitkLBFGSBOptimizer.h"
21 #include "vnl/algo/vnl_lbfgsb.h"
22 #include <vnl/vnl_math.h>
23
24 extern "C" {
25   extern double v3p_netlib_dpmeps_();
26 }
27
28 namespace clitk
29 {
30
31
32 /** \class LBFGSBOptimizerHelper
33  * \brief Wrapper helper around vnl_lbfgsb.
34  *
35  * This class is used to translate iteration events, etc, from
36  * vnl_lbfgsb into iteration events in ITK.
37  */
38 class ITK_EXPORT LBFGSBOptimizerHelper :
39     public vnl_lbfgsb
40 {
41 public:
42   typedef LBFGSBOptimizerHelper               Self;
43   typedef vnl_lbfgsb                          Superclass;
44
45   /** Create with a reference to the ITK object */
46   LBFGSBOptimizerHelper( vnl_cost_function& f,
47                          LBFGSBOptimizer* itkObj );
48
49   /** Handle new iteration event */
50   virtual bool report_iter();
51
52 private:
53   LBFGSBOptimizer* m_ItkObj;
54 };
55
56   
57
58
59 /**
60  * Constructor
61  */
62 LBFGSBOptimizer
63 ::LBFGSBOptimizer()
64 {
65   m_OptimizerInitialized = false;
66   m_VnlOptimizer         = 0;
67
68   m_LowerBound       = InternalBoundValueType(0);
69   m_UpperBound       = InternalBoundValueType(0); 
70   m_BoundSelection   = InternalBoundSelectionType(0);
71
72   m_CostFunctionConvergenceFactor   = 1e+7;
73   m_ProjectedGradientTolerance      = 1e-5;
74   m_MaximumNumberOfIterations       = 500;
75   m_MaximumNumberOfEvaluations      = 500;
76   m_MaximumNumberOfCorrections      = 5;
77   m_CurrentIteration                = 0;
78   m_Value                           = 0.0;
79   m_InfinityNormOfProjectedGradient = 0.0;
80   m_StopConditionDescription.str("");
81  
82 }
83
84
85 /**
86  * Destructor
87  */
88 LBFGSBOptimizer
89 ::~LBFGSBOptimizer()
90 {
91   delete m_VnlOptimizer;
92 }
93
94 /**
95  * PrintSelf
96  */
97 void
98 LBFGSBOptimizer
99 ::PrintSelf( std::ostream& os, itk::Indent indent) const
100 {  
101   Superclass::PrintSelf(os, indent);
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     m_Value << 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
138 /**
139  * Set lower bound
140  */
141 void
142 LBFGSBOptimizer
143 ::SetLowerBound(
144 const BoundValueType& value )
145 {
146   m_LowerBound = value;
147   if( m_OptimizerInitialized )
148     {
149     m_VnlOptimizer->set_lower_bound( m_LowerBound );
150     }
151
152   this->Modified();
153 }
154
155 /**
156  * Get lower bound
157  */
158 const
159 LBFGSBOptimizer 
160 ::BoundValueType &
161 LBFGSBOptimizer
162 ::GetLowerBound()
163 {
164   return m_LowerBound;
165
166
167 /**
168  * Set upper bound
169  */
170 void
171 LBFGSBOptimizer
172 ::SetUpperBound(
173 const BoundValueType& value )
174 {
175   m_UpperBound = value;
176   if( m_OptimizerInitialized )
177     {
178     m_VnlOptimizer->set_upper_bound( m_UpperBound );
179     }
180
181   this->Modified();
182 }
183
184 /**
185  * Get upper bound
186  */
187 const
188 LBFGSBOptimizer 
189 ::BoundValueType &
190 LBFGSBOptimizer
191 ::GetUpperBound()
192 {
193   return m_UpperBound;
194
195
196
197 /**
198  * Set bound selection array
199  */
200 void
201 LBFGSBOptimizer
202 ::SetBoundSelection(
203 const BoundSelectionType& value )
204 {
205   m_BoundSelection = value;
206   if( m_OptimizerInitialized )
207     {
208     m_VnlOptimizer->set_bound_selection( m_BoundSelection );
209     }
210   this->Modified();
211 }
212
213 /**
214  * Get bound selection array
215  */
216 const
217 LBFGSBOptimizer 
218 ::BoundSelectionType &
219 LBFGSBOptimizer
220 ::GetBoundSelection()
221 {
222   return m_BoundSelection;
223
224
225
226 /** Set/Get the CostFunctionConvergenceFactor. Algorithm terminates
227  * when the reduction in cost function is less than factor * epsmcj
228  * where epsmch is the machine precision.
229  * Typical values for factor: 1e+12 for low accuracy; 
230  * 1e+7 for moderate accuracy and 1e+1 for extremely high accuracy.
231  */
232 void
233 LBFGSBOptimizer
234 ::SetCostFunctionConvergenceFactor( double value )
235 {
236   if( value < 1.0 )
237     {
238     itkExceptionMacro("Value " << value 
239       << " is too small for SetCostFunctionConvergenceFactor()"
240       << "a typical range would be from 1.0 to 1e+12");
241     }
242   m_CostFunctionConvergenceFactor = value;
243   if( m_OptimizerInitialized )
244     {
245     m_VnlOptimizer->set_cost_function_convergence_factor(
246       m_CostFunctionConvergenceFactor );
247     }
248   this->Modified();
249 }
250
251
252 /** Set/Get the ProjectedGradientTolerance. Algorithm terminates
253  * when the project gradient is below the tolerance. Default value
254  * is 1e-5.
255  */
256 void
257 LBFGSBOptimizer
258 ::SetProjectedGradientTolerance( double value )
259 {
260   m_ProjectedGradientTolerance = value;
261   if( m_OptimizerInitialized )
262     {
263     m_VnlOptimizer->set_projected_gradient_tolerance(
264       m_ProjectedGradientTolerance );
265     }
266   this->Modified();
267 }
268
269
270 /** Set/Get the MaximumNumberOfIterations. Default is 500 */
271 void
272 LBFGSBOptimizer
273 ::SetMaximumNumberOfIterations( unsigned int value )
274 {
275   m_MaximumNumberOfIterations = value;
276   this->Modified();
277 }
278
279
280 /** Set/Get the MaximumNumberOfEvaluations. Default is 500 */
281 void
282 LBFGSBOptimizer
283 ::SetMaximumNumberOfEvaluations( unsigned int value )
284 {
285   m_MaximumNumberOfEvaluations = value;
286   if( m_OptimizerInitialized )
287     {
288     m_VnlOptimizer->set_max_function_evals( m_MaximumNumberOfEvaluations );
289     }
290   this->Modified();
291 }
292
293
294 /** Set/Get the MaximumNumberOfCorrections. Default is 5 */
295 void
296 LBFGSBOptimizer
297 ::SetMaximumNumberOfCorrections( unsigned int value )
298 {
299   m_MaximumNumberOfCorrections = value;
300   if( m_OptimizerInitialized )
301     {
302     m_VnlOptimizer->set_max_variable_metric_corrections(
303       m_MaximumNumberOfCorrections );
304     }
305   this->Modified();
306 }
307
308
309 /**
310  * Connect a Cost Function
311  */
312 void
313 LBFGSBOptimizer::
314 SetCostFunction( itk::SingleValuedCostFunction * costFunction )
315 {
316   m_CostFunction = costFunction;
317
318   const unsigned int numberOfParameters = 
319     costFunction->GetNumberOfParameters();
320
321   CostFunctionAdaptorType * adaptor = 
322     new CostFunctionAdaptorType( numberOfParameters );
323        
324   adaptor->SetCostFunction( costFunction );
325
326   if( m_OptimizerInitialized )
327     { 
328     delete m_VnlOptimizer;
329     }
330     
331   this->SetCostFunctionAdaptor( adaptor );
332
333   m_VnlOptimizer = new InternalOptimizerType( *adaptor, this );
334
335   // set the optimizer parameters
336   m_VnlOptimizer->set_lower_bound( m_LowerBound );
337   m_VnlOptimizer->set_upper_bound( m_UpperBound );
338   m_VnlOptimizer->set_bound_selection( m_BoundSelection );
339   m_VnlOptimizer->set_cost_function_convergence_factor(
340     m_CostFunctionConvergenceFactor );
341   m_VnlOptimizer->set_projected_gradient_tolerance( 
342     m_ProjectedGradientTolerance );
343   m_VnlOptimizer->set_max_function_evals( m_MaximumNumberOfEvaluations );
344   m_VnlOptimizer->set_max_variable_metric_corrections(
345     m_MaximumNumberOfCorrections );
346
347   m_OptimizerInitialized = true;
348
349   this->Modified();
350 }
351
352
353 /**
354  * Start the optimization
355  */
356 void
357 LBFGSBOptimizer
358 ::StartOptimization( void )
359 {
360   
361   // Check if all the bounds parameters are the same size as the initial parameters.
362   unsigned int numberOfParameters = m_CostFunction->GetNumberOfParameters();
363
364   if ( this->GetInitialPosition().Size() < numberOfParameters )
365     {
366     itkExceptionMacro( << "InitialPosition array does not have sufficient number of elements" );
367     }
368
369   if ( m_LowerBound.size() < numberOfParameters )
370     {
371     itkExceptionMacro( << "LowerBound array does not have sufficient number of elements" );
372     }
373
374   if ( m_UpperBound.size() < numberOfParameters )
375     {
376     itkExceptionMacro( << "UppperBound array does not have sufficient number of elements" );
377     }
378
379   if ( m_BoundSelection.size() < numberOfParameters )
380     {
381     itkExceptionMacro( << "BoundSelection array does not have sufficient number of elements" );
382     }
383
384   if( this->GetMaximize() )
385     {
386     this->GetNonConstCostFunctionAdaptor()->NegateCostFunctionOn();
387     }
388
389   this->SetCurrentPosition( this->GetInitialPosition() );
390
391   ParametersType parameters( this->GetInitialPosition() );
392
393   // Clear the description
394   m_StopConditionDescription.str("");
395
396   // vnl optimizers return the solution by reference 
397   // in the variable provided as initial position
398   m_VnlOptimizer->minimize( parameters );
399
400   if ( parameters.size() != this->GetInitialPosition().Size() )
401     {
402     // set current position to initial position and throw an exception
403     this->SetCurrentPosition( this->GetInitialPosition() );
404     itkExceptionMacro( << "Error occured in optimization" );
405     }
406
407   this->SetCurrentPosition( parameters );
408
409   this->InvokeEvent( itk::EndEvent() );
410
411 }
412
413
414 /*-------------------------------------------------------------------------
415  * helper class
416  *-------------------------------------------------------------------------
417  */
418
419 /** Create with a reference to the ITK object */
420 LBFGSBOptimizerHelper
421 ::LBFGSBOptimizerHelper( vnl_cost_function& f,
422                          LBFGSBOptimizer* itkObj )
423   : vnl_lbfgsb( f ),
424     m_ItkObj( itkObj )
425 {
426 }
427
428
429 /** Handle new iteration event */
430 bool
431 LBFGSBOptimizerHelper
432 ::report_iter()
433 {
434   Superclass::report_iter();
435
436   m_ItkObj->m_InfinityNormOfProjectedGradient =
437     this->get_inf_norm_projected_gradient();
438
439   m_ItkObj->InvokeEvent( itk::IterationEvent() );
440
441   m_ItkObj->m_CurrentIteration = this->num_iterations_;
442
443   //JV
444   m_ItkObj->m_Value = m_ItkObj->GetCachedValue();
445
446
447   // Return true to terminate the optimization loop.
448   if( this->num_iterations_ > m_ItkObj->m_MaximumNumberOfIterations )
449     {
450     m_ItkObj->m_StopConditionDescription <<
451       "Too many iterations. Iterations = " <<
452       m_ItkObj->m_MaximumNumberOfIterations;
453     return true;
454     }
455   else
456     {
457     return false;
458     }
459 }
460
461 const std::string
462 LBFGSBOptimizer::
463 GetStopConditionDescription() const
464 {
465   if (m_VnlOptimizer)
466     {
467     m_StopConditionDescription.str("");
468     m_StopConditionDescription << this->GetNameOfClass() << ": ";
469     switch (m_VnlOptimizer->get_failure_code())
470       {
471       case vnl_nonlinear_minimizer::ERROR_FAILURE:
472         m_StopConditionDescription << "Failure";
473         break;
474       case vnl_nonlinear_minimizer::ERROR_DODGY_INPUT:
475         m_StopConditionDescription << "Dodgy input";
476         break;
477       case vnl_nonlinear_minimizer::CONVERGED_FTOL:
478         m_StopConditionDescription << "Function tolerance reached after "
479                                    << m_CurrentIteration
480                                    << " iterations. "
481                                    << "The relative reduction of the cost function <= "
482                                    << m_CostFunctionConvergenceFactor * v3p_netlib_dpmeps_()
483                                    << " = CostFunctionConvergenceFactor ("
484                                    << m_CostFunctionConvergenceFactor
485                                    << ") * machine precision ("
486                                    << v3p_netlib_dpmeps_()
487                                    << ").";
488         break;
489       case vnl_nonlinear_minimizer::CONVERGED_XTOL:
490         m_StopConditionDescription << "Solution tolerance reached";
491         break;
492       case vnl_nonlinear_minimizer::CONVERGED_XFTOL:
493         m_StopConditionDescription << "Solution and Function tolerance both reached";
494         break;
495       case vnl_nonlinear_minimizer::CONVERGED_GTOL:
496         m_StopConditionDescription << "Gradient tolerance reached. "
497                                    << "Projected gradient tolerance is "
498                                    << m_ProjectedGradientTolerance;
499         break;
500       case vnl_nonlinear_minimizer::FAILED_TOO_MANY_ITERATIONS:
501         m_StopConditionDescription << "Too many iterations. Iterations = "
502                                    << m_MaximumNumberOfEvaluations;
503         break;
504       case vnl_nonlinear_minimizer::FAILED_FTOL_TOO_SMALL:
505         m_StopConditionDescription << "Function tolerance too small";
506         break;
507       case vnl_nonlinear_minimizer::FAILED_XTOL_TOO_SMALL:
508         m_StopConditionDescription << "Solution tolerance too small";
509         break;
510       case vnl_nonlinear_minimizer::FAILED_GTOL_TOO_SMALL:
511         m_StopConditionDescription << "Gradient tolerance too small";
512         break;
513       case vnl_nonlinear_minimizer::FAILED_USER_REQUEST:
514         break;
515       }
516     return m_StopConditionDescription.str();
517     }
518   else
519     {
520     return std::string("");
521     }
522 }
523
524 } // end namespace clitk
525
526 #endif