SPSAOptimizer

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