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