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://oncora1.lyon.fnclcc.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
  19 /*=========================================================================
  20
  21   Program:   Insight Segmentation & Registration Toolkit
  22   Module:    $RCSfile: clitkLBFGSBOptimizer.cxx,v $
  23   Language:  C++
  24   Date:      $Date: 2010/06/14 17:32:06 $
  25   Version:   $Revision: 1.1 $
  26
  27   Copyright (c) Insight Software Consortium. All rights reserved.
  28   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
  29
  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.
  33
  34 =========================================================================*/
  35 #ifndef __clitkLBFGSBOptimizer_cxx
  36 #define __clitkLBFGSBOptimizer_cxx
  37
  38 #include "clitkLBFGSBOptimizer.h"
  39 #include "vnl/algo/vnl_lbfgsb.h"
  40 #include <vnl/vnl_math.h>
  41
  42
  43 namespace clitk
  44 {
  45
  46
  47 /** \class LBFGSBOptimizerHelper
  48  * \brief Wrapper helper around vnl_lbfgsb.
  49  *
  50  * This class is used to translate iteration events, etc, from
  51  * vnl_lbfgsb into iteration events in ITK.
  52  */
  53 class ITK_EXPORT LBFGSBOptimizerHelper :
  54   public vnl_lbfgsb
  55 {
  56 public:
  57   typedef LBFGSBOptimizerHelper               Self;
  58   typedef vnl_lbfgsb                          Superclass;
  59
  60   /** Create with a reference to the ITK object */
  61   LBFGSBOptimizerHelper( vnl_cost_function& f,
  62                          LBFGSBOptimizer* itkObj );
  63
  64   /** Handle new iteration event */
  65   virtual bool report_iter();
  66
  67 private:
  68   LBFGSBOptimizer* m_ItkObj;
  69 };
  70
  71
  72
  73
  74 /**
  75  * Constructor
  76  */
  77 LBFGSBOptimizer
  78 ::LBFGSBOptimizer()
  79 {
  80   m_OptimizerInitialized = false;
  81   m_VnlOptimizer         = 0;
  82
  83   m_LowerBound       = InternalBoundValueType(0);
  84   m_UpperBound       = InternalBoundValueType(0);
  85   m_BoundSelection   = InternalBoundSelectionType(0);
  86
  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;
  93   m_Value                           = 0.0;
  94   m_InfinityNormOfProjectedGradient = 0.0;
  95
  96 }
  97
  98
  99 /**
 100  * Destructor
 101  */
 102 LBFGSBOptimizer
 103 ::~LBFGSBOptimizer()
 104 {
 105   delete m_VnlOptimizer;
 106 }
 107
 108 /**
 109  * PrintSelf
 110  */
 111 void
 112 LBFGSBOptimizer
 113 ::PrintSelf( std::ostream& os, itk::Indent indent) const
 114 {
 115   Superclass::PrintSelf(os, indent);
 116
 117   os << indent << "LowerBound: " << m_LowerBound << std::endl;
 118   os << indent << "UpperBound: " << m_UpperBound << std::endl;
 119   os << indent << "BoundSelection: " << m_BoundSelection << std::endl;
 120
 121   os << indent << "CostFunctionConvergenceFactor: " <<
 122      m_CostFunctionConvergenceFactor << std::endl;
 123
 124   os << indent << "ProjectedGradientTolerance: " <<
 125      m_ProjectedGradientTolerance << std::endl;
 126
 127   os << indent << "MaximumNumberOfIterations: " <<
 128      m_MaximumNumberOfIterations << std::endl;
 129
 130   os << indent << "MaximumNumberOfEvaluations: " <<
 131      m_MaximumNumberOfEvaluations << std::endl;
 132
 133   os << indent << "MaximumNumberOfCorrections: " <<
 134      m_MaximumNumberOfCorrections << std::endl;
 135
 136   os << indent << "CurrentIteration: " <<
 137      m_CurrentIteration << std::endl;
 138
 139   os << indent << "Value: " <<
 140      m_Value << std::endl;
 141
 142   os << indent << "InfinityNormOfProjectedGradient: " <<
 143      m_InfinityNormOfProjectedGradient << std::endl;
 144
 145   if( this->m_VnlOptimizer ) {
 146     os << indent << "Vnl LBFGSB Failure Code: " <<
 147        this->m_VnlOptimizer->get_failure_code() << std::endl;
 148   }
 149 }
 150
 151 /**
 152  * Set lower bound
 153  */
 154 void
 155 LBFGSBOptimizer
 156 ::SetLowerBound(
 157   const BoundValueType& value )
 158 {
 159   m_LowerBound = value;
 160   if( m_OptimizerInitialized ) {
 161     m_VnlOptimizer->set_lower_bound( m_LowerBound );
 162   }
 163
 164   this->Modified();
 165 }
 166
 167 /**
 168  * Get lower bound
 169  */
 170 const
 171 LBFGSBOptimizer
 172 ::BoundValueType &
 173 LBFGSBOptimizer
 174 ::GetLowerBound()
 175 {
 176   return m_LowerBound;
 177 }
 178
 179 /**
 180  * Set upper bound
 181  */
 182 void
 183 LBFGSBOptimizer
 184 ::SetUpperBound(
 185   const BoundValueType& value )
 186 {
 187   m_UpperBound = value;
 188   if( m_OptimizerInitialized ) {
 189     m_VnlOptimizer->set_upper_bound( m_UpperBound );
 190   }
 191
 192   this->Modified();
 193 }
 194
 195 /**
 196  * Get upper bound
 197  */
 198 const
 199 LBFGSBOptimizer
 200 ::BoundValueType &
 201 LBFGSBOptimizer
 202 ::GetUpperBound()
 203 {
 204   return m_UpperBound;
 205 }
 206
 207
 208 /**
 209  * Set bound selection array
 210  */
 211 void
 212 LBFGSBOptimizer
 213 ::SetBoundSelection(
 214   const BoundSelectionType& value )
 215 {
 216   m_BoundSelection = value;
 217   if( m_OptimizerInitialized ) {
 218     m_VnlOptimizer->set_bound_selection( m_BoundSelection );
 219   }
 220   this->Modified();
 221 }
 222
 223 /**
 224  * Get bound selection array
 225  */
 226 const
 227 LBFGSBOptimizer
 228 ::BoundSelectionType &
 229 LBFGSBOptimizer
 230 ::GetBoundSelection()
 231 {
 232   return m_BoundSelection;
 233 }
 234
 235
 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.
 241  */
 242 void
 243 LBFGSBOptimizer
 244 ::SetCostFunctionConvergenceFactor( double value )
 245 {
 246   if( value < 1.0 ) {
 247     itkExceptionMacro("Value " << value
 248                       << " is too small for SetCostFunctionConvergenceFactor()"
 249                       << "a typical range would be from 1.0 to 1e+12");
 250   }
 251   m_CostFunctionConvergenceFactor = value;
 252   if( m_OptimizerInitialized ) {
 253     m_VnlOptimizer->set_cost_function_convergence_factor(
 254       m_CostFunctionConvergenceFactor );
 255   }
 256   this->Modified();
 257 }
 258
 259
 260 /** Set/Get the ProjectedGradientTolerance. Algorithm terminates
 261  * when the project gradient is below the tolerance. Default value
 262  * is 1e-5.
 263  */
 264 void
 265 LBFGSBOptimizer
 266 ::SetProjectedGradientTolerance( double value )
 267 {
 268   m_ProjectedGradientTolerance = value;
 269   if( m_OptimizerInitialized ) {
 270     m_VnlOptimizer->set_projected_gradient_tolerance(
 271       m_ProjectedGradientTolerance );
 272   }
 273   this->Modified();
 274 }
 275
 276
 277 /** Set/Get the MaximumNumberOfIterations. Default is 500 */
 278 void
 279 LBFGSBOptimizer
 280 ::SetMaximumNumberOfIterations( unsigned int value )
 281 {
 282   m_MaximumNumberOfIterations = value;
 283   this->Modified();
 284 }
 285
 286
 287 /** Set/Get the MaximumNumberOfEvaluations. Default is 500 */
 288 void
 289 LBFGSBOptimizer
 290 ::SetMaximumNumberOfEvaluations( unsigned int value )
 291 {
 292   m_MaximumNumberOfEvaluations = value;
 293   if( m_OptimizerInitialized ) {
 294     m_VnlOptimizer->set_max_function_evals( m_MaximumNumberOfEvaluations );
 295   }
 296   this->Modified();
 297 }
 298
 299
 300 /** Set/Get the MaximumNumberOfCorrections. Default is 5 */
 301 void
 302 LBFGSBOptimizer
 303 ::SetMaximumNumberOfCorrections( unsigned int value )
 304 {
 305   m_MaximumNumberOfCorrections = value;
 306   if( m_OptimizerInitialized ) {
 307     m_VnlOptimizer->set_max_variable_metric_corrections(
 308       m_MaximumNumberOfCorrections );
 309   }
 310   this->Modified();
 311 }
 312
 313
 314 /**
 315  * Connect a Cost Function
 316  */
 317 void
 318 LBFGSBOptimizer::
 319 SetCostFunction( itk::SingleValuedCostFunction * costFunction )
 320 {
 321   m_CostFunction = costFunction;
 322
 323   const unsigned int numberOfParameters =
 324     costFunction->GetNumberOfParameters();
 325
 326   CostFunctionAdaptorType * adaptor =
 327     new CostFunctionAdaptorType( numberOfParameters );
 328
 329   adaptor->SetCostFunction( costFunction );
 330
 331   if( m_OptimizerInitialized ) {
 332     delete m_VnlOptimizer;
 333   }
 334
 335   this->SetCostFunctionAdaptor( adaptor );
 336
 337   m_VnlOptimizer = new InternalOptimizerType( *adaptor, this );
 338
 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 );
 350
 351   m_OptimizerInitialized = true;
 352
 353   this->Modified();
 354 }
 355
 356
 357 /**
 358  * Start the optimization
 359  */
 360 void
 361 LBFGSBOptimizer
 362 ::StartOptimization( void )
 363 {
 364
 365   // Check if all the bounds parameters are the same size as the initial parameters.
 366   unsigned int numberOfParameters = m_CostFunction->GetNumberOfParameters();
 367
 368   if ( this->GetInitialPosition().Size() < numberOfParameters ) {
 369     itkExceptionMacro( << "InitialPosition array does not have sufficient number of elements" );
 370   }
 371
 372   if ( m_LowerBound.size() < numberOfParameters ) {
 373     itkExceptionMacro( << "LowerBound array does not have sufficient number of elements" );
 374   }
 375
 376   if ( m_UpperBound.size() < numberOfParameters ) {
 377     itkExceptionMacro( << "UppperBound array does not have sufficient number of elements" );
 378   }
 379
 380   if ( m_BoundSelection.size() < numberOfParameters ) {
 381     itkExceptionMacro( << "BoundSelection array does not have sufficient number of elements" );
 382   }
 383
 384   if( this->GetMaximize() ) {
 385     this->GetNonConstCostFunctionAdaptor()->NegateCostFunctionOn();
 386   }
 387
 388   this->SetCurrentPosition( this->GetInitialPosition() );
 389
 390   ParametersType parameters( this->GetInitialPosition() );
 391
 392   // vnl optimizers return the solution by reference
 393   // in the variable provided as initial position
 394   m_VnlOptimizer->minimize( parameters );
 395
 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" );
 400   }
 401
 402   this->SetCurrentPosition( parameters );
 403
 404   this->InvokeEvent( itk::EndEvent() );
 405
 406 }
 407
 408
 409 /*-------------------------------------------------------------------------
 410  * helper class
 411  *-------------------------------------------------------------------------
 412  */
 413
 414 /** Create with a reference to the ITK object */
 415 LBFGSBOptimizerHelper
 416 ::LBFGSBOptimizerHelper( vnl_cost_function& f,
 417                          LBFGSBOptimizer* itkObj )
 418   : vnl_lbfgsb( f ),
 419     m_ItkObj( itkObj )
 420 {
 421 }
 422
 423
 424 /** Handle new iteration event */
 425 bool
 426 LBFGSBOptimizerHelper
 427 ::report_iter()
 428 {
 429   Superclass::report_iter();
 430
 431   m_ItkObj->m_InfinityNormOfProjectedGradient =
 432     this->get_inf_norm_projected_gradient();
 433
 434   m_ItkObj->InvokeEvent( itk::IterationEvent() );
 435
 436   m_ItkObj->m_CurrentIteration = this->num_iterations_;
 437
 438   //JV
 439   m_ItkObj->m_Value = m_ItkObj->GetCachedValue();
 440
 441
 442   // Return true to terminate the optimization loop.
 443   if( this->num_iterations_ > m_ItkObj->m_MaximumNumberOfIterations ) {
 444     return true;
 445   } else {
 446     return false;
 447   }
 448 }
 449
 450
 451 } // end namespace clitk
 452
 453 #endif