optimizerNelderMead.h

Go to the documentation of this file.

#ifndef _DLR_OPTIMIZERNELDERMEAD_H_
#define _DLR_OPTIMIZERNELDERMEAD_H_

#include <functional>
#include <vector>
#include <dlrOptimization/optimizer.h>
#include <dlrCommon/exception.h>

namespace dlr {

  namespace optimization {

    // template <std::unary_function Functor>
    template <class Functor>
    class OptimizerNelderMead
      : public Optimizer<Functor> {
    public:

      // Typedefs for convenience
      typedef typename Functor::argument_type argument_type;
      typedef typename Functor::result_type result_type;

    
      OptimizerNelderMead();

    
      explicit OptimizerNelderMead(const Functor& functor);

    
      OptimizerNelderMead(const OptimizerNelderMead& source);

    
      ~OptimizerNelderMead();

    
      virtual std::vector<size_t>
      getNumberOfFunctionCalls();


      void
      setDelta(const argument_type& delta);

    
      void
      setNumberOfRestarts(size_t numberOfRestarts) {
        m_numberOfRestarts = numberOfRestarts;
      }

    
      void
      setParameters(argument_type delta,
                    size_t functionCallLimit=5000,
                    size_t numberOfRestarts=1,
                    double alpha=1.0,
                    double beta=0.5,
                    double gamma=2.0,
                    double minimumSimplexValueSpan=0.0001,
                    size_t verbosity=0);

    
      virtual void
      setStartPoint(argument_type startPoint);

      virtual void
      setVerbosity(int verbosity) {m_verbosity = verbosity;}

      OptimizerNelderMead&
      operator=(const OptimizerNelderMead& source);

    protected:

      void
      computeAxisSums(
        const std::vector<argument_type>& currentPoints,
        argument_type& axisSums);


      void
      doNelderMead(std::vector<argument_type>& currentPoints,
                   std::vector<result_type>& currentValues,
                   size_t& numberOfFunctionCalls);


      result_type
      evaluateMove(std::vector<argument_type>& currentPoints,
                   std::vector<result_type>& currentValues,
                   const argument_type& axisSums,
                   double factor);


      std::pair<typename Functor::argument_type, typename Functor::result_type>
      run();

      argument_type m_delta;
      size_t m_functionCallLimit;
      size_t m_numberOfRestarts;
      double m_alpha;
      double m_beta;
      double m_gamma;
      double m_minimumSimplexValueSpan;
      bool m_deltaValueHack;

      std::vector<size_t> m_functionCallCount;
      argument_type m_theta0;
      size_t m_verbosity;

    }; // class OptimizerNelderMead

  } // namespace optimization

} // namespace dlr


/* ======= Declarations to maintain compatibility with legacy code. ======= */

namespace dlr {

  using optimization::OptimizerNelderMead;

} // namespace dlr


/*******************************************************************
 * Member function definitions follow.  This would be a .cpp file
 * if it weren't templated.
 *******************************************************************/

#include <cmath>

namespace dlr {

  namespace optimization {

    namespace privateCode {

      template <class Type0, class Type1>
      class extractSecond
        : public std::unary_function<std::pair<Type0, Type1>, Type1>
      {
      public:
        Type1 operator()(const std::pair<Type0, Type1>& inputPair) {
          return inputPair.second;
        }
      };

      template <class Type>
      std::vector<size_t>
      argsort(const std::vector<Type>& inputVector)
      {
        std::vector< std::pair<Type, size_t> >
          sortVector(inputVector.size());
        std::vector<size_t> resultVector(inputVector.size());

        for(size_t index = 0; index < inputVector.size(); ++index) {
          sortVector[index] =
            std::pair<Type, size_t>(inputVector[index], index);
        }
        std::sort(sortVector.begin(), sortVector.end());
        std::transform(sortVector.begin(), sortVector.end(), resultVector.begin(),
                       extractSecond<Type, size_t>());
        return resultVector;    
      }

      template <class Type>
      std::vector<Type>
      take(const std::vector<Type>& inputVector,
            const std::vector<size_t>& indexVector)
      {
        if(inputVector.size() != indexVector.size()) {
          std::ostringstream message;
          message << "Incorrectly Sized Arguments:\n"
                  << " First Argument: " <<inputVector.size()
                  << ", Second Argument: " << indexVector.size() << ".";
          DLR_THROW(dlr::ValueException, "take()",
                    message.str().c_str());
        }
        std::vector<Type> resultVector(inputVector.size());
        for(size_t index = 0; index < resultVector.size(); ++index) {
          if(indexVector[index] >= inputVector.size()) {
            std::ostringstream message;
            message << "Out of Range Index:"
                    << indexVector[index] << ".";
            DLR_THROW(dlr::ValueException, "take()",
                      message.str().c_str());
          }
          resultVector[index] = inputVector[indexVector[index]];
        }
        return resultVector;
      }

    } // namespace privateCode

  } // namespace optimization
  
} // namespace dlr


namespace dlr {

  namespace optimization {

    // Sets parameters to reasonable values for functions which take
    // values and arguments in the "normal" range of 0 to 100 or so.
    // template <std::unary_function Functor>
    template <class Functor>
    OptimizerNelderMead<Functor>::
    OptimizerNelderMead()
      : Optimizer<Functor>(),
        m_functionCallCount(),
        m_verbosity(0)
    {
      this->setParameters(argument_type());
      m_deltaValueHack = true;
    }


    // Constructor which specifies the specific Functor instance to
    // use.
    // template <std::unary_function Functor>
    template <class Functor>
    OptimizerNelderMead<Functor>::
    OptimizerNelderMead(const Functor& functor)
      : Optimizer<Functor>(functor),
        m_functionCallCount(),
        m_verbosity(0)
    {
      this->setParameters(argument_type());
      m_deltaValueHack = true;
    }

    // Copy constructor.
    // template <std::unary_function Functor>
    template <class Functor>
    OptimizerNelderMead<Functor>::
    OptimizerNelderMead(const OptimizerNelderMead& source)
      : Optimizer<Functor>(source),
        m_functionCallLimit(source.m_functionCallLimit),
        m_numberOfRestarts(source.m_numberOfRestarts),
        m_alpha(source.m_alpha),
        m_beta(source.m_beta),
        m_gamma(source.m_gamma),
        m_minimumSimplexValueSpan(source.m_minimumSimplexValueSpan),
        m_deltaValueHack(source.m_deltaValueHack),
        m_functionCallCount(source.m_functionCallCount),
        m_verbosity(source.m_verbosity)
    {
      copyArgumentType(source.m_delta, m_delta);
    }

    // Destructor.
    // template <std::unary_function Functor>
    template <class Functor>
    OptimizerNelderMead<Functor>::
    ~OptimizerNelderMead()
    {
      // Empty
    }

    // Queries the number of functionCalls required to complete the
    // previous minimization.
    // template <std::unary_function Functor>
    template <class Functor>
    std::vector<size_t>
    OptimizerNelderMead<Functor>::
    getNumberOfFunctionCalls()
    {
      return m_functionCallCount;
    }


    // This method sets the spacing of the initial points used in the
    // nonlinear optimization, without affecting any other optimization
    // parameters.
    template <class Functor>
    void
    OptimizerNelderMead<Functor>::
    setDelta(const typename OptimizerNelderMead<Functor>::argument_type& delta)
    {
      copyArgumentType(delta, m_delta);
    }
  
  
    // Sets minimization parameters.
    // template <std::unary_function Functor>
    template <class Functor>
    void
    OptimizerNelderMead<Functor>::
    setParameters(typename OptimizerNelderMead<Functor>::argument_type delta,
                  size_t functionCallLimit,
                  size_t numberOfRestarts,
                  double alpha,
                  double beta,
                  double gamma,
                  double minimumSimplexValueSpan,
                  size_t verbosity)
    {
      copyArgumentType(delta, m_delta);
      m_functionCallLimit = functionCallLimit;
      m_numberOfRestarts = numberOfRestarts;
      m_alpha = alpha;
      m_beta = beta;
      m_gamma = gamma;
      m_minimumSimplexValueSpan = minimumSimplexValueSpan;
      m_verbosity = verbosity;
      m_deltaValueHack = false;

      // Inherited member
      Optimizer<Functor>::m_needsOptimization = true;
    }

    // Sets the initial conditions for the minimization.
    // template <std::unary_function Functor>
    template <class Functor>
    void
    OptimizerNelderMead<Functor>::
    setStartPoint(argument_type startPoint)
    {
      m_theta0 = startPoint;

      // Inherited member
      Optimizer<Functor>::m_needsOptimization = true;
    }
                  
    // Assignment operator.
    // template <std::unary_function Functor>
    template <class Functor>
    OptimizerNelderMead<Functor>&
    OptimizerNelderMead<Functor>::
    operator=(const OptimizerNelderMead<Functor>& source)
    {
      Optimizer<Functor>::operator=(source);
      m_delta = source.m_delta;
      m_functionCallLimit = source.m_functionCallLimit;
      m_numberOfRestarts = source.m_numberOfRestarts;
      m_alpha = source.m_alpha;
      m_beta = source.m_beta;
      m_gamma = source.m_gamma;
      m_minimumSimplexValueSpan = source.m_minimumSimplexValueSpan;
      m_deltaValueHack = source.m_deltaValueHack;
      m_functionCallCount = source.m_functionCallCount;
    }

    // template <std::unary_function Functor>
    template <class Functor>
    void
    OptimizerNelderMead<Functor>::
    computeAxisSums(
      const std::vector<typename OptimizerNelderMead<Functor>::argument_type>&
      currentPoints,
      typename OptimizerNelderMead<Functor>::argument_type& axisSums)
    {
      if(currentPoints.size() == 0) {
        DLR_THROW(LogicException, "OptimizerNelderMead::computeAxisSums()",
                  "Vector Size of Current Points is 0... "
                  "something's not right.");
        return;
      }
      copyArgumentType(currentPoints[0], axisSums);
      for(size_t index = 1; index < currentPoints.size(); ++index) {
        axisSums += currentPoints[index];
      }
    }
  
    // Run the actual simplex search.  Modifies all arguments.
    // template <std::unary_function Functor>
    template <class Functor>
    void
    OptimizerNelderMead<Functor>::
    doNelderMead(
      std::vector<typename OptimizerNelderMead<Functor>::argument_type>&
      currentPoints,
      std::vector<typename OptimizerNelderMead<Functor>::result_type>&
      currentValues,
      size_t& numberOfFunctionCalls)
    {
      size_t dimension = currentValues.size() - 1;
      argument_type axisSums;
      this->computeAxisSums(currentPoints, axisSums);
      while(1) {
        std::vector<size_t> indices = privateCode::argsort(currentValues);
        currentValues = privateCode::take(currentValues, indices);
        currentPoints = privateCode::take(currentPoints, indices);

        if(m_verbosity >= 3) {
          std::cout << "\rCurrent best: " << currentValues[0] << std::flush;
        }

        double simplexValueSpan = 0;
        if((std::fabs(currentValues[dimension])
            + std::fabs(currentValues[0])) > 0) {
          simplexValueSpan =
            (2.0 * std::fabs(currentValues[dimension] - currentValues[0])
             / (std::fabs(currentValues[dimension])
                + std::fabs(currentValues[0])));
        }
        if(m_verbosity >= 3) {
          std::cout << "       \t       simplexValueSpan: "
                    << simplexValueSpan << "     "
                    << std::flush;
        }
        if((simplexValueSpan < m_minimumSimplexValueSpan)
           || (numberOfFunctionCalls >= m_functionCallLimit)) {
          if(m_verbosity >= 3) {
            std::cout << "\n";
          }
          if(m_verbosity >= 1) {
            std::cout << "Terminating with:\n"
                      << "  simplexValueSpan = " << simplexValueSpan
                      << " (" << m_minimumSimplexValueSpan << ")\n";
            std::cout << "  numberOfFunctionCalls = " << numberOfFunctionCalls
                      << " (" << m_functionCallLimit << ")" << std::endl;
          }
          break;
        }
        result_type newValue =
          this->evaluateMove(currentPoints, currentValues, axisSums, 
                             (-1.0 * m_alpha));
        numberOfFunctionCalls += 1;
        this->computeAxisSums(currentPoints, axisSums);
        if(newValue <= currentValues[0]) {
          newValue = evaluateMove(currentPoints, currentValues, axisSums, 
                                  m_gamma);
          numberOfFunctionCalls += 1;
          this->computeAxisSums(currentPoints, axisSums);
        } else if (newValue >= currentValues[dimension - 1]) {
          result_type oldMaxValue = currentValues[dimension];
          newValue = evaluateMove(currentPoints, currentValues, axisSums, 
                                  m_beta);
          numberOfFunctionCalls += 1;
          if(newValue >= oldMaxValue) {
            for(size_t i = 1; i < dimension + 1; ++i) {
              currentPoints[i] = 0.5 * (currentPoints[i] + currentPoints[0]);
              currentValues[i] = m_functor(currentPoints[i]);
              numberOfFunctionCalls += 1;
            }
          }
          this->computeAxisSums(currentPoints, axisSums);
        }
      }
    }


    // Evaluate, and maybe accept, a new point.  You'll want to believe
    // the following math before reading this code.
    // 
    // xMean + factor*(xMax - xMean)
    // = xMean + factor*xMax - factor*xMean
    // = (1 - factor)*xMean + factor*xMax
    // = ((1 - factor)/n)*sum(xSubi, xSubi != xMax) + factor*xMax
    // = ((1 - factor)/n)*sum(xSubi) - ((1 - factor)/n)*xMax + factor*xMax
    // = ((1 - factor)/n)*sum(xSubi) + (factor - ((1 - factor)/n))*xMax
    // template <std::unary_function Functor>
    template <class Functor>
    typename OptimizerNelderMead<Functor>::result_type
    OptimizerNelderMead<Functor>::
    evaluateMove(std::vector<argument_type>& currentPoints,
                 std::vector<result_type>& currentValues,
                 const argument_type& axisSums,
                 double factor)
    {
      size_t dimension = currentPoints.size() - 1;
      double centroidFactor = (1.0 - factor) / dimension;
      double extrapolationFactor = factor - centroidFactor;
      argument_type newPoint =
        ((centroidFactor * axisSums)
         + (extrapolationFactor * currentPoints[dimension]));
      result_type newValue = m_functor(newPoint);
      if(newValue < currentValues[dimension]) {
        currentValues[dimension] = newValue;
        currentPoints[dimension] = newPoint;
      }
      return newValue;
    }

    // Perform the minimization (top level).
    // template <std::unary_function Functor>
    template <class Functor>
    std::pair<typename Functor::argument_type, typename Functor::result_type>
    OptimizerNelderMead<Functor>::
    run()
    {
      size_t dimension = m_theta0.size();
      if(dimension == 0) {
        DLR_THROW(ValueException, "OptimizerNelderMead::run()",
                  "invalid starting point has zero size.");
      }
      if(m_deltaValueHack) {
        // Initialize delta using only operations that we expect to be
        // supported by argtype, and which we don't expect to be
        // defined in surprising ways.
        m_delta = (m_theta0 * 0.0) + 1.0;
        m_deltaValueHack = false;
      }
    
      // The algorithm requires dimension + 1 initial points and values.
      std::vector<argument_type> currentPoints(dimension + 1);
      std::vector<result_type> currentValues(dimension + 1);
      // First point is the one specified by the user.
      copyArgumentType(m_theta0, currentPoints[0]);
      currentValues[0] = m_functor(currentPoints[0]);

      // We'll repeat the whole minimization several times.
      m_functionCallCount.clear();
      for(size_t i = 0; i < m_numberOfRestarts + 1; ++i) {
        // The remaining points will be built using m_delta
        for(size_t j = 0; j < dimension; ++j) {
          copyArgumentType(currentPoints[0], currentPoints[j + 1]);
          (currentPoints[j + 1])[j] += m_delta[j];
          currentValues[j + 1] = m_functor(currentPoints[j + 1]);
        }
        size_t numberOfFunctionCalls = (i == 0) ? 1 : 0;
        numberOfFunctionCalls += dimension;
        // Run the minimization algorithm.
        this->doNelderMead(currentPoints, currentValues, numberOfFunctionCalls);
        // Have to add dimension to numberOfFunctionCalls because of the
        // initialization steps above.
        m_functionCallCount.push_back(numberOfFunctionCalls + dimension);
      }
      return std::make_pair(currentPoints[0], currentValues[0]);
    }

  } // namespace optimization

} // namespace dlr


#endif // #ifndef _DLR_OPTIMIZERNELDERMEAD_H_

---