freAmoebaSOOptimizer.cxx

Go to the documentation of this file.

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

  Program:   F.R.E.E. - flexible registration evaluation engine
  Version:   v.1.0.0
  Date:      $Date: 2006/09/01 12:00:00 $
  Module:    $RCSfile: freAmoebaSOOptimizer.cxx,v $
  Language:  C++


  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.

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

#include "freAmoebaSOOptimizer.h"
#include "itkCommand.h"
#include "itkEventObject.h"

namespace FREE
{

AmoebaSOOptimizer
::AmoebaSOOptimizer()
  : m_InitialSimplexDelta(1)  // initial size
{
  m_OptimizerInitialized           = false;
  m_VnlOptimizer                   = 0;
  m_MaximumNumberOfIterations      = 500;
  m_ParametersConvergenceTolerance = 1e-8;
  m_FunctionConvergenceTolerance   = 1e-4;
  m_AutomaticInitialSimplex        = true;
        m_InitialSimplexDelta.Fill(itk::NumericTraits<ParametersType::ValueType>::One);
}

AmoebaSOOptimizer
::~AmoebaSOOptimizer()
{
  if (m_VnlOptimizer) delete  m_VnlOptimizer;
}

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

  os << indent << "MaximumNumberOfIterations: " 
     << m_MaximumNumberOfIterations << std::endl;
  os << indent << "ParametersConvergenceTolerance: "
     << m_ParametersConvergenceTolerance << std::endl;
  os << indent << "FunctionConvergenceTolerance: "
     << m_FunctionConvergenceTolerance << std::endl;
  os << indent << "AutomaticInitialSimplex: "
     << (m_AutomaticInitialSimplex ? "On" : "Off") << std::endl;
  os << indent << "InitialSimplexDelta: "
     << m_InitialSimplexDelta << std::endl;
}

AmoebaSOOptimizer::MeasureType
AmoebaSOOptimizer
::GetValue (const ParametersType &parameters)
{
  ParametersType scaledParameters = parameters;
  if(m_ScalesInitialized)
  {
    const ScalesType scales = this->GetScales();
    for(unsigned int i=0;i<parameters.size();i++)
    {
      scaledParameters[i] *= scales[i]; 
    }
  }
  return this->GetNonConstCostFunctionAdaptor()->f( scaledParameters );
};

AmoebaSOOptimizer::DecomposedMeasureType
AmoebaSOOptimizer
::GetDecomposedValue (const ParametersType &parameters)
{
        this->GetValue(parameters);
        return this->GetNonConstCostFunctionAdaptor()->GetCachedDecomposedValue();
};

void
AmoebaSOOptimizer
::SetMaximumNumberOfIterations( unsigned int n )
{
  if ( n == m_MaximumNumberOfIterations )
    {
    return;
    }

  m_MaximumNumberOfIterations = n;
  if ( m_OptimizerInitialized )
    {
    m_VnlOptimizer->set_max_iterations( static_cast<int>( n ) );
    }
  this->Modified();
}

void
AmoebaSOOptimizer
::SetParametersConvergenceTolerance( double tol )
{
  if ( tol == m_ParametersConvergenceTolerance )
    {
    return;
    }

  m_ParametersConvergenceTolerance = tol;
  if ( m_OptimizerInitialized )
    {
    m_VnlOptimizer->set_x_tolerance( tol );
    }

  this->Modified();
}

void
AmoebaSOOptimizer
::SetFunctionConvergenceTolerance( double tol )
{
  if ( tol == m_FunctionConvergenceTolerance )
    {
    return;
    }

  m_FunctionConvergenceTolerance = tol;
  if ( m_OptimizerInitialized )
    {
    m_VnlOptimizer->set_f_tolerance( tol );
    }

  this->Modified();
}

void
AmoebaSOOptimizer::
SetCostFunction( CostFunctionType * costFunction )
{
        const unsigned int numberOfParameters = 
                costFunction->GetNumberOfParameters();

        CostFunctionAdaptorType * adaptor = 
                new CostFunctionAdaptorType( numberOfParameters );

        this->m_CostFunction = costFunction;
        adaptor->SetCostFunction( costFunction );

        if( m_OptimizerInitialized )
        { 
                delete m_VnlOptimizer;
        }

        this->SetCostFunctionAdaptor( adaptor );

        m_VnlOptimizer = new vnl_amoeba( *adaptor );

        // set up optimizer parameters
        m_VnlOptimizer->set_max_iterations( static_cast<int>( m_MaximumNumberOfIterations ) );
        m_VnlOptimizer->set_x_tolerance( m_ParametersConvergenceTolerance );
        m_VnlOptimizer->set_f_tolerance( m_FunctionConvergenceTolerance );

        m_OptimizerInitialized = true;
}

void
AmoebaSOOptimizer
::StartOptimization( void )
{
    
        this->InvokeEvent( itk::StartEvent() );

  if( this->GetMaximize() )
    {
    this->GetNonConstCostFunctionAdaptor()->NegateCostFunctionOn();
    }

  ParametersType initialPosition = this->GetInitialPosition();

  ParametersType parameters( initialPosition );

  // If the user provides the scales then we set otherwise we don't
  // for computation speed.
  // We also scale the initial parameters up if scales are defined.
  // This compensates for later scaling them down in the cost function adaptor
  // and at the end of this function.  
  if(m_ScalesInitialized)
    {
    ScalesType scales = this->GetScales();
    this->GetNonConstCostFunctionAdaptor()->SetScales(scales);
    for(unsigned int i=0;i<parameters.size();i++)
      {
      parameters[i] *= scales[i]; 
      }
    }

  // vnl optimizers return the solution by reference 
  // in the variable provided as initial position
  if (m_AutomaticInitialSimplex)
    {
    m_VnlOptimizer->minimize( parameters );
    }
  else
    {
    InternalParametersType delta( m_InitialSimplexDelta );
    // m_VnlOptimizer->verbose = 1;
    m_VnlOptimizer->minimize( parameters, delta );
    }
  
   // we scale the parameters down if scales are defined
  if(m_ScalesInitialized)
    {
    ScalesType scales = this->GetScales();
    for(unsigned int i=0;i<parameters.size();i++)
      {
      parameters[i] /= scales[i]; 
      }
    }

  this->SetCurrentPosition( parameters );

        m_BestPosition = parameters;
        m_BestValue     = this->GetCachedValue();
    
        this->InvokeEvent( itk::EndEvent() );
}

vnl_amoeba * 
AmoebaSOOptimizer
::GetOptimizer()
{
  return m_VnlOptimizer;
}

} // end namespace FREE

---