SPSAOptimizer

From KitwarePublic
Revision as of 13:36, 17 March 2005 by Blezek (talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigationJump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.

Dan Blezek's implemetation of the SPSA optimizer. Because I couldn't upload the files as media, here they are:

itkSPSAOptimizer.h

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkSPSAOptimizer.h,v $
  Language:  C++
  Date:      $Date: 2005/03/10 13:32:40 $
  Version:   $Revision: 1.1 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkSPSAOptimizer_h
#define __itkSPSAOptimizer_h

#include "itkSingleValuedNonLinearOptimizer.h"
#include <log4cxx/logger.h>

namespace itk
{
  
/** \class SPSAOptimizer
 * \brief Implements the Simultaseous Perturbation Stochastic Approximation optimizer
 * \ingroup Numerics Optimizers
 */  
class ITK_EXPORT SPSAOptimizer : 
    public SingleValuedNonLinearOptimizer
{
public:
  /** Standard class typedefs. */
  typedef SPSAOptimizer          Self;
  typedef SingleValuedNonLinearOptimizer    Superclass;
  typedef SmartPointer<Self>                Pointer;
  typedef SmartPointer<const Self>          ConstPointer;

  /**  ParametersType typedef.
   *  It defines a position in the optimization search space. */
  typedef Superclass::ParametersType         ParametersType;
  
  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods). */
  itkTypeMacro( SPSAOptimizer, SingleValuedNonLinearOptimizer );

  /** Methods to configure the cost function. */
  itkGetConstMacro( Maximize, bool );
  itkSetMacro( Maximize, bool );
  itkBooleanMacro( Maximize );
  bool GetMinimize( ) const
  { return !m_Maximize; }
  void SetMinimize(bool v)
  { this->SetMaximize(!v); }
  void MinimizeOn()
  { this->MaximizeOff(); }
  void MinimizeOff()
  { this->MaximizeOn(); }
  
  /** Advance one step following the gradient direction. */
  virtual void AdvanceOneStep( void );

  /** Start optimization. */
  void    StartOptimization( void );

  /** Resume previously stopped optimization with current parameters
   * \sa StopOptimization. */
  void    ResumeOptimization( void );

  /** Stop optimization.
   * \sa ResumeOptimization */
  void    StopOptimization( void );

  /** Set the number of iterations. */
  itkSetMacro( NumberOfIterations, unsigned long );

  /** Get the number of iterations. */
  itkGetConstMacro( NumberOfIterations, unsigned long );

  /** Get the current iteration number. */
  itkGetConstMacro( CurrentIteration, unsigned int );

  /** Get the current value. */
  itkGetConstMacro( Value, double );
  
  /** Get the current value. */
  itkSetMacro( Alpha, double );
  itkGetConstMacro( Alpha, double );

  /** Get the current value. */
  itkSetMacro( Gamma, double );
  itkGetConstMacro( Gamma, double );

  /** Get the current value. */
  itkSetMacro( a, double );
  itkGetConstMacro( a, double );

  /** Get the current value. */
  itkSetMacro( A, double );
  itkGetConstMacro( A, double );

  /** Get the current value. */
  itkSetMacro( c, double );
  itkGetConstMacro( c, double );

protected:
  SPSAOptimizer();
  virtual ~SPSAOptimizer() {};
  void PrintSelf(std::ostream& os, Indent indent) const;

  // made protected so subclass can access
  bool                          m_Maximize;

 private:
  static log4cxx::LoggerPtr logger;

  SPSAOptimizer(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented
  
  bool                          m_Stop;
  double                        m_Value;
  unsigned long                 m_NumberOfIterations;
  unsigned long                 m_CurrentIteration;
  double m_Alpha, m_Gamma, m_a, m_A, m_c;
  ParametersType m_CurrentPosition;
  
};

} // end namespace itk


#endif



and itkSPSAOPtimizer.cxx

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkSPSAOptimizer.cxx,v $
  Language:  C++
  Date:      $Date: 2005/03/10 13:32:40 $
  Version:   $Revision: 1.1 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef _itkSPSAOptimizer_txx
#define _itkSPSAOptimizer_txx

#include "itkSPSAOptimizer.h"
#include "itkCommand.h"
#include "itkEventObject.h"
#include "itkExceptionObject.h"
#include "vnl/vnl_math.h"
#include "vnl/vnl_sample.h"


namespace itk
{
  log4cxx::LoggerPtr SPSAOptimizer::logger;

/**
 * Constructor
 */
SPSAOptimizer
::SPSAOptimizer()
{
  if ( logger == NULL )
    {
    logger = log4cxx::Logger::getLogger ( "org.itk.Code.Numerics.SPSAOptimizer" );
    }
  logger->debug ( "Constructor" );
  itkDebugMacro("Constructor");

  m_NumberOfIterations = 100;
  m_CurrentIteration = 0;
  m_Maximize = false;
  m_Value = 0.0;
  m_Alpha = 0.602;
  m_Gamma = 0.101;
  m_a = 0.1;
  m_A = 0.1;
  m_c = 0.1;
}



void
SPSAOptimizer
::PrintSelf(std::ostream& os, Indent indent) const
{
  Superclass::PrintSelf(os,indent);

  os << indent << "NunberOfIterations: "
     << m_NumberOfIterations << std::endl;
  os << indent << "Maximize: "
     << m_Maximize << std::endl;
  os << indent << "CurrentIteration: "
     << m_CurrentIteration;
  os << indent << "Value: "
     << m_Value;
  if (m_CostFunction)
    {
    os << indent << "CostFunction: "
       << m_CostFunction;
    }
  os << indent << "Alpha: " << m_Alpha << std::endl;
  os << indent << "Gamma: " << m_Gamma << std::endl;
  os << indent << "a: " << m_a << std::endl;
  os << indent << "A: " << m_A << std::endl;
  os << indent << "c: " << m_c << std::endl;
  os << std::endl;
}


/**
 * Start the optimization
 */
void
SPSAOptimizer
::StartOptimization( void )
{

  itkDebugMacro("StartOptimization");
  logger->debug ( "StartOptimization" );   
  m_CurrentIteration   = 0;

  this->SetCurrentPosition( this->GetInitialPosition() );
  this->ResumeOptimization();

}



/**
 * Resume the optimization
 */
void
SPSAOptimizer
::ResumeOptimization( void )
{
  
  itkDebugMacro("ResumeOptimization");

  m_Stop = false;

  InvokeEvent( StartEvent() );
  while( !m_Stop ) 
    {

    try
      {
      m_Value = m_CostFunction->GetValue( this->GetCurrentPosition() );
      }
    catch( ExceptionObject& err )
      {
      // An exception has occurred. 
      // Terminate immediately.
      StopOptimization();

      // Pass exception to caller
      throw err;
      }


    if( m_Stop )
      {
      break;
      }
  
    AdvanceOneStep();

    m_CurrentIteration++;

    if( m_CurrentIteration >= m_NumberOfIterations )
      {
      StopOptimization();
      break;
      }
    
    }
    

}


/**
 * Stop optimization
 */
void
SPSAOptimizer
::StopOptimization( void )
{

  itkDebugMacro("StopOptimization");

  m_Stop = true;
  InvokeEvent( EndEvent() );
}





/**
 * Advance one Step following the gradient direction
 */
void
SPSAOptimizer
::AdvanceOneStep( void )
{ 
  char buffer[1024];
  itkDebugMacro("AdvanceOneStep");
  logger->debug ( "AdvanceOneStep" );
  double direction;
  if( this->m_Maximize ) 
    {
    direction = 1.0;
    }
  else 
    {
    direction = -1.0;
    }

  const unsigned int spaceDimension = 
    m_CostFunction->GetNumberOfParameters();

  const ParametersType & currentPosition = this->GetCurrentPosition();

  ScalesType scales = this->GetScales();

  // Make sure the scales have been set properly
  if (scales.size() != spaceDimension)
    {
    itkExceptionMacro(<< "The size of Scales is "
                      << scales.size()
                      << ", but the NumberOfParameters for the CostFunction is "
                      << spaceDimension
                      << ".");
    }


  double ak, ck, yplus, yminus;
  double a = m_a, c = m_c, A = m_A;

  ak = a / pow (A + m_CurrentIteration, m_Alpha);
  ck = c / pow ( m_CurrentIteration+1, m_Gamma );
  ParametersType delta ( spaceDimension );
  ParametersType theta ( spaceDimension );
  
  sprintf ( buffer, "ak: %f ck: %f", ak, ck );
  logger->debug ( buffer );

  for ( unsigned int i = 0; i < spaceDimension; i++ ) {
    delta[i] = 2 * vnl_math_rnd ( vnl_sample_uniform ( 0.0f, 1.0f ) ) - 1;
    delta[i] = delta[i] / scales[i];
  }

  for ( unsigned int i = 0; i < spaceDimension; i++ ) {
    theta[i] = currentPosition[i] + ck * delta[i];
  }

  yplus = m_CostFunction->GetValue ( theta );
  
  for ( unsigned int i = 0; i < spaceDimension; i++ ) {
    theta[i] = currentPosition[i] - ck * delta[i];
  }

  yminus = m_CostFunction->GetValue ( theta );

  sprintf ( buffer, "New Position: " );
  double ydiff = direction * ( yplus - yminus );
  for ( unsigned int i = 0; i < spaceDimension; i++ ) {
    theta[i] = currentPosition[i] - ak * ydiff / ( 2.0 * ck * delta[i] );
    sprintf ( buffer, "%s %f, ", buffer, theta[i] );
  }
  logger->debug ( buffer ); 
  this->SetCurrentPosition( theta );

  this->InvokeEvent( IterationEvent() );

}



} // end namespace itk

#endif