registration/clitkGenericOptimizer.h

   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 #ifndef __clitkGenericOptimizer_h
  20 #define __clitkGenericOptimizer_h
  21
  22 //clitk include
  23 #include "clitkLBFGSBOptimizer.h"
  24
  25 //itk include
  26 #include "itkAmoebaOptimizer.h"
  27 #include "itkPowellOptimizer.h"
  28 #include "itkFRPROptimizer.h"
  29 #include "itkRegularStepGradientDescentOptimizer.h"
  30 #include "itkVersorRigid3DTransformOptimizer.h"
  31 #include "itkConjugateGradientOptimizer.h"
  32 #include "itkLBFGSOptimizer.h"
  33 #include "itkSPSAOptimizer.h"
  34
  35 // general
  36 #include <iomanip>
  37
  38 /*Requires at least this section in the associated *.ggo file: (modify defaults to fit application)
  39
  40 option "optimizer"      -       "0=Simplex, 1=Powell, 2=FRPR, 3=Regular Step GD, 4=VersorRigid3D, 5=Conjugated Gradient, 6=L-BFGS, 7=L-BFGS-B" int no default="0"
  41 option "delta"          -       "0: Initial delta, otherwise automatic"                                                 double  no
  42 option "step"           -       "1,2,3,4: Initial stepsize (to be multiplied with the gradient)"                        double  no      default="2.0"
  43 option "relax"          -       "3,4: Relaxation of the stepsize (multiplied each time the gradient changes sign)"      double  no      default="0.7"
  44 option "valueTol"       -       "0,1,2: Tolerance on the function"                                                      double  no      default="0.01"
  45 option "stepTol"        -       "0,1,3,4: Tolerance on the step size"                                                   double  no      default="0.1"
  46 option "gradTol"        -       "3,4,6,7: Tolerance on the (projected) gradient magnitude (7: 1=low->1e-10=high precision)"     double  no      default="1e-5"
  47 option "lineAcc"        -       "6: Line accuracy (eg: high=0.1, low=0.9)"      double  no      default="0.9"
  48 option "convFactor"     -       "7: Convergence factor: terminate if factor*machine_precision>reduction in cost (1e+12 low -> 1e+1 high precision) "    double          no      default="1e+7"
  49 option "maxIt"          -       "0-7: Maximum number of iterations"                                     int     no      default="500"
  50 option "maxLineIt"      -       "1,2: Maximum number of line iterations"                                        int     no      default="50"
  51 option "maxEval"        -       "6,7: Maximum number of evaluations"                                            int     no      default="500"
  52 option "maxCorr"        -       "7: Maximum number of corrections"                                              int     no      default="5"
  53 option "selectBound"    -       "7: Select the type of bound: 0=none, 1=u, 2=u&l, 3=l"                  int     no      default="0"
  54 option "lowerBound"     -       "7: The lower bound"                                                    double  no
  55 option "upperBound"     -       "7: The upper bound"                                                    double  no
  56
  57 */
  58
  59 namespace clitk
  60 {
  61 template<class args_info_type>
  62 class GenericOptimizer : public itk::LightObject
  63 {
  64 public:
  65   //==============================================
  66   typedef GenericOptimizer     Self;
  67   typedef itk::LightObject     Superclass;
  68   typedef itk::SmartPointer<Self>            Pointer;
  69   typedef itk::SmartPointer<const Self>      ConstPointer;
  70
  71   /** Method for creation through the object factory. */
  72   itkNewMacro(Self);
  73
  74
  75   typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
  76   typedef OptimizerType::Pointer OptimizerPointer;
  77
  78   //==============================================
  79   //Set members
  80   void SetArgsInfo(args_info_type args_info) {
  81     m_ArgsInfo= args_info;
  82     m_Verbose=m_ArgsInfo.verbose_flag;
  83   }
  84   void SetMaximize(bool f) {
  85     m_Maximize=f;
  86     m_MaximizeGiven=true;
  87   }
  88
  89   void SetNumberOfParameters(unsigned int n) {
  90     m_NumberOfParameters=n;
  91   }
  92
  93   //==============================================
  94   //Get members
  95   OptimizerPointer GetOptimizerPointer(void) {
  96
  97     if (!m_MaximizeGiven)std::cerr<<"Warning: Maximize/Minimize was not given to the optimizer!!"<<std::endl;
  98     else if(m_Verbose) {
  99       if (m_Maximize)  std::cout<<"Maximizing similarity measure..."<<std::endl;
 100       else std::cout<<"Minimizing similarity measure..."<<std::endl;
 101     }
 102
 103
 104     //============================================================================
 105     //switch on  the  metric type chosen and set metric specific members
 106     switch (m_ArgsInfo.optimizer_arg) {
 107
 108     case 0: {
 109       itk::AmoebaOptimizer::Pointer m =
 110         itk::AmoebaOptimizer::New();
 111       m_OptimizerIsVNL=true;
 112
 113       //Set Parameters for this optimizer
 114       m->SetMaximize(m_Maximize);
 115       if (m_ArgsInfo.delta_given) {
 116         itk::Array<double> delta(m_NumberOfParameters);
 117         for (unsigned int i =0 ; i < m_NumberOfParameters; i ++)
 118           delta[i]=m_ArgsInfo.delta_arg;
 119         m->SetInitialSimplexDelta(delta);
 120         m->SetAutomaticInitialSimplex(false);
 121       } else m->SetAutomaticInitialSimplex(true);
 122       m->SetMaximumNumberOfIterations(m_ArgsInfo.maxIt_arg);
 123       m->SetFunctionConvergenceTolerance(m_ArgsInfo.valueTol_arg);
 124       m->SetParametersConvergenceTolerance(m_ArgsInfo.stepTol_arg);
 125       if (m_Verbose) std::cout<<"Using Simplex Optimizer..."<<std::endl;
 126       m_Optimizer=m;
 127       break;
 128     }
 129
 130     case 1: {
 131       itk::PowellOptimizer::Pointer m =
 132         itk::PowellOptimizer::New();
 133       m_OptimizerIsVNL=false;
 134
 135       //Set Parameters for this optimizer
 136       m->SetMaximize(m_Maximize);
 137       m->SetStepLength(m_ArgsInfo.step_arg);
 138       m->SetStepTolerance(m_ArgsInfo.stepTol_arg);
 139       m->SetValueTolerance (m_ArgsInfo.valueTol_arg);
 140       m->SetMaximumIteration(m_ArgsInfo.maxIt_arg);
 141       m->SetMaximumLineIteration(m_ArgsInfo.maxLineIt_arg);
 142       if (m_Verbose) std::cout<<"Using Powell Optimizer..."<<std::endl;
 143       m_Optimizer=m;
 144       break;
 145     }
 146
 147     case 2: {
 148       itk::FRPROptimizer::Pointer m =
 149         itk::FRPROptimizer::New();
 150       m_OptimizerIsVNL=false;
 151
 152       //Set Parameters for this optimizer
 153       m->SetMaximize(m_Maximize);
 154       m->SetStepLength(m_ArgsInfo.step_arg);
 155       m->SetStepTolerance(m_ArgsInfo.stepTol_arg);
 156       m->SetValueTolerance (m_ArgsInfo.valueTol_arg);
 157       m->SetMaximumIteration(m_ArgsInfo.maxIt_arg);
 158       m->SetMaximumLineIteration(m_ArgsInfo.maxLineIt_arg);
 159       if (m_Verbose) std::cout<<"Using Powell Optimizer..."<<std::endl;
 160       m_Optimizer=m;
 161       break;
 162     }
 163
 164     case 3: {
 165       itk::RegularStepGradientDescentOptimizer::Pointer m =
 166         itk::RegularStepGradientDescentOptimizer::New();
 167       m_OptimizerIsVNL=false;
 168
 169       //Set Parameters for this optimizer
 170       m->SetMaximize(m_Maximize);
 171       m->SetMinimumStepLength( m_ArgsInfo.stepTol_arg );
 172       m->SetRelaxationFactor(  m_ArgsInfo.relax_arg );
 173       m->SetMaximumStepLength( m_ArgsInfo.step_arg );
 174       m->SetNumberOfIterations(m_ArgsInfo.maxIt_arg  );
 175       m->SetGradientMagnitudeTolerance(m_ArgsInfo.gradTol_arg  );
 176       if (m_Verbose) std::cout<<"Using Regular Step Gradient Descent Optimizer..."<<std::endl;
 177       m_Optimizer=m;
 178       break;
 179     }
 180
 181     case 4: {
 182       itk::VersorRigid3DTransformOptimizer::Pointer m
 183       = itk::VersorRigid3DTransformOptimizer::New();
 184       m_OptimizerIsVNL=false;
 185
 186       //Set Parameters for this optimizer
 187       m->SetMaximize(m_Maximize);
 188       m->SetMinimumStepLength( m_ArgsInfo.stepTol_arg );
 189       m->SetRelaxationFactor(  m_ArgsInfo.relax_arg );
 190       m->SetMaximumStepLength( m_ArgsInfo.step_arg );
 191       m->SetNumberOfIterations(m_ArgsInfo.maxIt_arg  );
 192       m->SetGradientMagnitudeTolerance(m_ArgsInfo.gradTol_arg  );
 193       if (m_Verbose) std::cout<<"Using VersorRigid3DTransform Optimizer..."<<std::endl;
 194       m_Optimizer=m;
 195       break;
 196     }
 197
 198     case 5: {
 199
 200       itk::ConjugateGradientOptimizer::Pointer m
 201       = itk::ConjugateGradientOptimizer::New();
 202       m_OptimizerIsVNL=true;
 203
 204       //Set Parameters for this optimizer
 205       m->SetMaximize(m_Maximize);
 206       if (m_Verbose) std::cout<<"Using Conjugated Gradient Optimizer..."<<std::endl;
 207       m_Optimizer=m;
 208       break;
 209     }
 210
 211     case 6: {
 212       itk::LBFGSOptimizer::Pointer m
 213       = itk::LBFGSOptimizer::New();
 214       m_OptimizerIsVNL=true;
 215
 216       //Set Parameters for this optimizer
 217       m->SetMaximize(m_Maximize);
 218       m->SetGradientConvergenceTolerance(m_ArgsInfo.gradTol_arg);
 219       m->SetLineSearchAccuracy(m_ArgsInfo.lineAcc_arg);
 220       m->SetMaximumNumberOfFunctionEvaluations(m_ArgsInfo.maxEval_arg);
 221       m->SetTrace(true);
 222       if (m_Verbose) std::cout<<"Using L-BFGS Optimizer..."<<std::endl;
 223       m_Optimizer=m;
 224       break;
 225     }
 226
 227     case 7: {
 228       LBFGSBOptimizer::Pointer m
 229       = LBFGSBOptimizer::New();
 230       m_OptimizerIsVNL=true;
 231
 232       //Set Parameters for this optimizer
 233       m->SetMaximize(m_Maximize);
 234       m->SetMaximumNumberOfEvaluations(m_ArgsInfo.maxEval_arg);
 235       m->SetMaximumNumberOfIterations(m_ArgsInfo.maxIt_arg);
 236       m->SetMaximumNumberOfCorrections(m_ArgsInfo.maxCorr_arg);
 237       m->SetCostFunctionConvergenceFactor( m_ArgsInfo.convFactor_arg);
 238       m->SetProjectedGradientTolerance( m_ArgsInfo.gradTol_arg);
 239
 240       // Bounds
 241       clitk::LBFGSBOptimizer::BoundSelectionType boundSelect( m_NumberOfParameters );
 242       clitk::LBFGSBOptimizer::BoundValueType upperBound( m_NumberOfParameters );
 243       clitk::LBFGSBOptimizer::BoundValueType lowerBound( m_NumberOfParameters );
 244       boundSelect.Fill( m_ArgsInfo.selectBound_arg );
 245       upperBound.Fill( m_ArgsInfo.upperBound_arg );
 246       lowerBound.Fill( m_ArgsInfo.lowerBound_arg );
 247       m->SetBoundSelection( boundSelect );
 248       m->SetUpperBound( upperBound );
 249       m->SetLowerBound( lowerBound );
 250
 251       if (m_Verbose) std::cout<<"Using L-BFGS-B Optimizer..."<<std::endl;
 252       m_Optimizer=m;
 253       break;
 254     }
 255
 256     //     case 8:
 257     //       {
 258     //  itk::SPSAOptimizer::Pointer m =
 259     //    itk::SPSAOptimizerr::Pointer::New();
 260     //      m_OptimizerIsVNL=false;
 261
 262     //  //Set Parameters for this optimizer
 263     //  m->SetMaximize(m_Maximize);
 264     //  m->SetGamma(m_ArgsInfo.gamma_arg);
 265     //  m->SetAlpha(m_ArgsInfo.gamma_arg);
 266     //  m->SetA(m_ArgsInfo.bigA_arg);
 267     //  m->Seta(m_ArgsInfo.littleA_arg);
 268     //  m->Setc(m_ArgsInfo.littleC_arg);
 269     //  m->SetMaximumNumberOfIterations(m_ArgsInfo.maxIt_arg);
 270     //  m->SetMinimumNumberOfIterations(m_ArgsInfo.minIt_arg);
 271     //  m->SetNumberOfPerturbations(m_ArgsInfo.per_arg);
 272     //  m->SetStateOfConvergenceDecayRate(args_inf.conv_arg);
 273     //  if (m_Verbose) std::cout<<"Using SPSA Optimizer..."<<std::endl;
 274     //  m_Optimizer=m;
 275     //  break;
 276     //       }
 277
 278     }
 279
 280
 281     //============================================================================
 282     //return the pointer
 283     return m_Optimizer;
 284   }
 285
 286
 287   void OutputIterationInfo(void) const {
 288     switch (m_ArgsInfo.optimizer_arg) {
 289     case 0: {
 290       itk::AmoebaOptimizer* o=dynamic_cast<itk::AmoebaOptimizer*>(m_Optimizer.GetPointer());
 291       std::cout<<std::setprecision(6);
 292       //if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetIteration()<<std::endl;
 293       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCachedCurrentPosition()<<std::endl;
 294       if (m_OutputValue)  std::cout<<"Value: "<< o->GetCachedValue()<<std::endl;
 295       // if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetCachedDerivative()<<std::endl;
 296       std::cout<<std::endl;
 297       break;
 298     }
 299     case 1: {
 300       itk::PowellOptimizer* o=dynamic_cast<itk::PowellOptimizer*>(m_Optimizer.GetPointer());
 301       std::cout<<std::setprecision(6);
 302       if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
 303       if (m_OutputIteration)  std::cout<<"Line Iteration: "<< o->GetCurrentLineIteration()<<std::endl;
 304       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCurrentPosition()<<std::endl;
 305       if (m_OutputValue)  std::cout<<"Value: "<< o->GetCurrentCost()<<std::endl;
 306       //if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetDerivative()<<std::endl;
 307       std::cout<<std::endl;
 308       break;
 309     }
 310     case 2: {
 311       itk::FRPROptimizer* o=dynamic_cast<itk::FRPROptimizer*>(m_Optimizer.GetPointer());
 312       std::cout<<std::setprecision(6);
 313       if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
 314       if (m_OutputIteration)  std::cout<<"Line Iteration: "<< o->GetCurrentLineIteration()<<std::endl;
 315       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCurrentPosition()<<std::endl;
 316       if (m_OutputValue)  std::cout<<"Value: "<< o->GetCurrentCost()<<std::endl;
 317       //if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetDerivative()<<std::endl;
 318       std::cout<<std::endl;
 319       break;
 320     }
 321     case 3: {
 322       itk::RegularStepGradientDescentOptimizer* o=dynamic_cast<itk::RegularStepGradientDescentOptimizer*>(m_Optimizer.GetPointer());
 323       std::cout<<std::setprecision(6);
 324       if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
 325       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCurrentPosition()<<std::endl;
 326       if (m_OutputPosition)  std::cout<<"Step length: "<< o->GetCurrentStepLength()<<std::endl;
 327       if (m_OutputValue)  std::cout<<"Value: "<< o->GetValue()<<std::endl;
 328       if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetGradient()<<std::endl;
 329       std::cout<<std::endl;
 330       break;
 331     }
 332     case 4: {
 333       itk::VersorRigid3DTransformOptimizer* o=dynamic_cast<itk::VersorRigid3DTransformOptimizer*>(m_Optimizer.GetPointer());
 334       std::cout<<std::setprecision(6);
 335       if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
 336       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCurrentPosition()<<std::endl;
 337       if (m_OutputPosition)  std::cout<<"Step length: "<< o->GetCurrentStepLength()<<std::endl;
 338       if (m_OutputValue)  std::cout<<"Value: "<< o->GetValue()<<std::endl;
 339       if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetGradient()<<std::endl;
 340       std::cout<<std::endl;
 341       break;
 342     }
 343     case 5: {
 344       itk::ConjugateGradientOptimizer* o=dynamic_cast<itk::ConjugateGradientOptimizer*>(m_Optimizer.GetPointer());
 345       std::cout<<std::setprecision(6);
 346       if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
 347       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCachedCurrentPosition()<<std::endl;
 348       if (m_OutputValue)  std::cout<<"Value: "<< o->GetCachedValue()<<std::endl;
 349       if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetCachedDerivative()<<std::endl;
 350       std::cout<<std::endl;
 351       break;
 352     }
 353     case 6: {
 354       itk::LBFGSOptimizer* o=dynamic_cast<itk::LBFGSOptimizer*>(m_Optimizer.GetPointer());
 355       std::cout<<std::setprecision(6);
 356       //if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
 357       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCachedCurrentPosition()<<std::endl;
 358       if (m_OutputValue)  std::cout<<"Value: "<< o->GetCachedValue()<<std::endl;
 359       //          if (m_OutputValue)  std::cout<<"Norm Projected Gradient: "<< o->GetInfinityNormOfProjectedGradient()<<std::endl;
 360       if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetCachedDerivative()<<std::endl;
 361       std::cout<<std::endl;
 362       break;
 363     }
 364     case 7: {
 365       clitk::LBFGSBOptimizer* o=dynamic_cast<clitk::LBFGSBOptimizer*>(m_Optimizer.GetPointer());
 366       std::cout<<std::setprecision(6);
 367       if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
 368       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCachedCurrentPosition()<<std::endl;
 369       if (m_OutputValue)  std::cout<<"Value: "<< o->GetValue()<<std::endl;
 370       if (m_OutputValue)  std::cout<<"Norm Projected Gradient:"<< o->GetInfinityNormOfProjectedGradient()<<std::endl;
 371       if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetCachedDerivative()<<std::endl;
 372       std::cout<<std::endl;
 373       break;
 374     }
 375     //       case 8:
 376     //  {
 377     //    itk::VersorRigid3DTransformOptimizer* o=dynamic_cast<itk::VersorRigid3DTransformOptimizer*>(m_Optimizer.GetPointer());
 378     //    std::cout<<std::setprecision(6);
 379     //    if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetIteration()<<std::endl;
 380     //    if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCurrentPosition()<<std::endl;
 381     //    if (m_OutputValue)  std::cout<<"Value: "<< o->GetCurrentCost()<<std::endl;
 382     //    if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetDerivative()<<std::endl;
 383     //    std::cout<<std::endl;
 384     //  }
 385     }
 386   }
 387
 388   void SetOutputIteration(bool o) {
 389     m_OutputIteration=o;
 390   }
 391   void SetOutputPosition(bool o) {
 392     m_OutputPosition=o;
 393   }
 394   void SetOutputValue(bool o) {
 395     m_OutputValue=o;
 396   }
 397   void SetOutputGradient(bool o) {
 398     m_OutputGradient=o;
 399   }
 400
 401
 402   //==============================================
 403 protected:
 404   GenericOptimizer() {
 405     m_Optimizer=NULL;
 406     m_Maximize=false;
 407     m_Verbose=false;
 408     m_MaximizeGiven=false;
 409     m_NumberOfParameters=1;
 410     m_OptimizerIsVNL=false;
 411     m_OutputIteration=true;
 412     m_OutputPosition=false;
 413     m_OutputValue=true;
 414     m_OutputGradient=false;
 415
 416   }
 417   ~GenericOptimizer() {};
 418
 419 private:
 420   args_info_type m_ArgsInfo;
 421   OptimizerPointer m_Optimizer;
 422   bool m_Maximize;
 423   bool m_Verbose;
 424   bool m_MaximizeGiven;
 425   unsigned int m_NumberOfParameters;
 426   bool m_OptimizerIsVNL;
 427   bool m_OutputIteration;
 428   bool m_OutputPosition;
 429   bool m_OutputValue;
 430   bool m_OutputGradient;
 431 };
 432
 433 } // end namespace clitk
 434
 435 #endif // #define __clitkGenericOptimizer_h