SPSAOptimizer
From KitwarePublic
Jump to navigationJump to search
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
