1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
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.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
19 /*=========================================================================
21 Program: Insight Segmentation & Registration Toolkit
22 Module: $RCSfile: clitkLBFGSBOptimizer.cxx,v $
24 Date: $Date: 2010/06/14 17:32:06 $
25 Version: $Revision: 1.1 $
27 Copyright (c) Insight Software Consortium. All rights reserved.
28 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
30 This software is distributed WITHOUT ANY WARRANTY; without even
31 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
32 PURPOSE. See the above copyright notices for more information.
34 =========================================================================*/
35 #ifndef __clitkLBFGSBOptimizer_cxx
36 #define __clitkLBFGSBOptimizer_cxx
38 #include "clitkLBFGSBOptimizer.h"
39 #include "vnl/algo/vnl_lbfgsb.h"
40 #include <vnl/vnl_math.h>
47 /** \class LBFGSBOptimizerHelper
48 * \brief Wrapper helper around vnl_lbfgsb.
50 * This class is used to translate iteration events, etc, from
51 * vnl_lbfgsb into iteration events in ITK.
53 class ITK_EXPORT LBFGSBOptimizerHelper :
57 typedef LBFGSBOptimizerHelper Self;
58 typedef vnl_lbfgsb Superclass;
60 /** Create with a reference to the ITK object */
61 LBFGSBOptimizerHelper( vnl_cost_function& f,
62 LBFGSBOptimizer* itkObj );
64 /** Handle new iteration event */
65 virtual bool report_iter();
68 LBFGSBOptimizer* m_ItkObj;
80 m_OptimizerInitialized = false;
83 m_LowerBound = InternalBoundValueType(0);
84 m_UpperBound = InternalBoundValueType(0);
85 m_BoundSelection = InternalBoundSelectionType(0);
87 m_CostFunctionConvergenceFactor = 1e+7;
88 m_ProjectedGradientTolerance = 1e-5;
89 m_MaximumNumberOfIterations = 500;
90 m_MaximumNumberOfEvaluations = 500;
91 m_MaximumNumberOfCorrections = 5;
92 m_CurrentIteration = 0;
94 m_InfinityNormOfProjectedGradient = 0.0;
105 delete m_VnlOptimizer;
113 ::PrintSelf( std::ostream& os, itk::Indent indent) const
115 Superclass::PrintSelf(os, indent);
117 os << indent << "LowerBound: " << m_LowerBound << std::endl;
118 os << indent << "UpperBound: " << m_UpperBound << std::endl;
119 os << indent << "BoundSelection: " << m_BoundSelection << std::endl;
121 os << indent << "CostFunctionConvergenceFactor: " <<
122 m_CostFunctionConvergenceFactor << std::endl;
124 os << indent << "ProjectedGradientTolerance: " <<
125 m_ProjectedGradientTolerance << std::endl;
127 os << indent << "MaximumNumberOfIterations: " <<
128 m_MaximumNumberOfIterations << std::endl;
130 os << indent << "MaximumNumberOfEvaluations: " <<
131 m_MaximumNumberOfEvaluations << std::endl;
133 os << indent << "MaximumNumberOfCorrections: " <<
134 m_MaximumNumberOfCorrections << std::endl;
136 os << indent << "CurrentIteration: " <<
137 m_CurrentIteration << std::endl;
139 os << indent << "Value: " <<
140 m_Value << std::endl;
142 os << indent << "InfinityNormOfProjectedGradient: " <<
143 m_InfinityNormOfProjectedGradient << std::endl;
145 if( this->m_VnlOptimizer ) {
146 os << indent << "Vnl LBFGSB Failure Code: " <<
147 this->m_VnlOptimizer->get_failure_code() << std::endl;
157 const BoundValueType& value )
159 m_LowerBound = value;
160 if( m_OptimizerInitialized ) {
161 m_VnlOptimizer->set_lower_bound( m_LowerBound );
185 const BoundValueType& value )
187 m_UpperBound = value;
188 if( m_OptimizerInitialized ) {
189 m_VnlOptimizer->set_upper_bound( m_UpperBound );
209 * Set bound selection array
214 const BoundSelectionType& value )
216 m_BoundSelection = value;
217 if( m_OptimizerInitialized ) {
218 m_VnlOptimizer->set_bound_selection( m_BoundSelection );
224 * Get bound selection array
228 ::BoundSelectionType &
230 ::GetBoundSelection()
232 return m_BoundSelection;
236 /** Set/Get the CostFunctionConvergenceFactor. Algorithm terminates
237 * when the reduction in cost function is less than factor * epsmcj
238 * where epsmch is the machine precision.
239 * Typical values for factor: 1e+12 for low accuracy;
240 * 1e+7 for moderate accuracy and 1e+1 for extremely high accuracy.
244 ::SetCostFunctionConvergenceFactor( double value )
247 itkExceptionMacro("Value " << value
248 << " is too small for SetCostFunctionConvergenceFactor()"
249 << "a typical range would be from 1.0 to 1e+12");
251 m_CostFunctionConvergenceFactor = value;
252 if( m_OptimizerInitialized ) {
253 m_VnlOptimizer->set_cost_function_convergence_factor(
254 m_CostFunctionConvergenceFactor );
260 /** Set/Get the ProjectedGradientTolerance. Algorithm terminates
261 * when the project gradient is below the tolerance. Default value
266 ::SetProjectedGradientTolerance( double value )
268 m_ProjectedGradientTolerance = value;
269 if( m_OptimizerInitialized ) {
270 m_VnlOptimizer->set_projected_gradient_tolerance(
271 m_ProjectedGradientTolerance );
277 /** Set/Get the MaximumNumberOfIterations. Default is 500 */
280 ::SetMaximumNumberOfIterations( unsigned int value )
282 m_MaximumNumberOfIterations = value;
287 /** Set/Get the MaximumNumberOfEvaluations. Default is 500 */
290 ::SetMaximumNumberOfEvaluations( unsigned int value )
292 m_MaximumNumberOfEvaluations = value;
293 if( m_OptimizerInitialized ) {
294 m_VnlOptimizer->set_max_function_evals( m_MaximumNumberOfEvaluations );
300 /** Set/Get the MaximumNumberOfCorrections. Default is 5 */
303 ::SetMaximumNumberOfCorrections( unsigned int value )
305 m_MaximumNumberOfCorrections = value;
306 if( m_OptimizerInitialized ) {
307 m_VnlOptimizer->set_max_variable_metric_corrections(
308 m_MaximumNumberOfCorrections );
315 * Connect a Cost Function
319 SetCostFunction( itk::SingleValuedCostFunction * costFunction )
321 m_CostFunction = costFunction;
323 const unsigned int numberOfParameters =
324 costFunction->GetNumberOfParameters();
326 CostFunctionAdaptorType * adaptor =
327 new CostFunctionAdaptorType( numberOfParameters );
329 adaptor->SetCostFunction( costFunction );
331 if( m_OptimizerInitialized ) {
332 delete m_VnlOptimizer;
335 this->SetCostFunctionAdaptor( adaptor );
337 m_VnlOptimizer = new InternalOptimizerType( *adaptor, this );
339 // set the optimizer parameters
340 m_VnlOptimizer->set_lower_bound( m_LowerBound );
341 m_VnlOptimizer->set_upper_bound( m_UpperBound );
342 m_VnlOptimizer->set_bound_selection( m_BoundSelection );
343 m_VnlOptimizer->set_cost_function_convergence_factor(
344 m_CostFunctionConvergenceFactor );
345 m_VnlOptimizer->set_projected_gradient_tolerance(
346 m_ProjectedGradientTolerance );
347 m_VnlOptimizer->set_max_function_evals( m_MaximumNumberOfEvaluations );
348 m_VnlOptimizer->set_max_variable_metric_corrections(
349 m_MaximumNumberOfCorrections );
351 m_OptimizerInitialized = true;
358 * Start the optimization
362 ::StartOptimization( void )
365 // Check if all the bounds parameters are the same size as the initial parameters.
366 unsigned int numberOfParameters = m_CostFunction->GetNumberOfParameters();
368 if ( this->GetInitialPosition().Size() < numberOfParameters ) {
369 itkExceptionMacro( << "InitialPosition array does not have sufficient number of elements" );
372 if ( m_LowerBound.size() < numberOfParameters ) {
373 itkExceptionMacro( << "LowerBound array does not have sufficient number of elements" );
376 if ( m_UpperBound.size() < numberOfParameters ) {
377 itkExceptionMacro( << "UppperBound array does not have sufficient number of elements" );
380 if ( m_BoundSelection.size() < numberOfParameters ) {
381 itkExceptionMacro( << "BoundSelection array does not have sufficient number of elements" );
384 if( this->GetMaximize() ) {
385 this->GetNonConstCostFunctionAdaptor()->NegateCostFunctionOn();
388 this->SetCurrentPosition( this->GetInitialPosition() );
390 ParametersType parameters( this->GetInitialPosition() );
392 // vnl optimizers return the solution by reference
393 // in the variable provided as initial position
394 m_VnlOptimizer->minimize( parameters );
396 if ( parameters.size() != this->GetInitialPosition().Size() ) {
397 // set current position to initial position and throw an exception
398 this->SetCurrentPosition( this->GetInitialPosition() );
399 itkExceptionMacro( << "Error occured in optimization" );
402 this->SetCurrentPosition( parameters );
404 this->InvokeEvent( itk::EndEvent() );
409 /*-------------------------------------------------------------------------
411 *-------------------------------------------------------------------------
414 /** Create with a reference to the ITK object */
415 LBFGSBOptimizerHelper
416 ::LBFGSBOptimizerHelper( vnl_cost_function& f,
417 LBFGSBOptimizer* itkObj )
424 /** Handle new iteration event */
426 LBFGSBOptimizerHelper
429 Superclass::report_iter();
431 m_ItkObj->m_InfinityNormOfProjectedGradient =
432 this->get_inf_norm_projected_gradient();
434 m_ItkObj->InvokeEvent( itk::IterationEvent() );
436 m_ItkObj->m_CurrentIteration = this->num_iterations_;
439 m_ItkObj->m_Value = m_ItkObj->GetCachedValue();
442 // Return true to terminate the optimization loop.
443 if( this->num_iterations_ > m_ItkObj->m_MaximumNumberOfIterations ) {
451 } // end namespace clitk