gradientFunction.h

Go to the documentation of this file.

#ifndef _DLR_GRADIENTFUNCTION_H_
#define _DLR_GRADIENTFUNCTION_H_

#include <functional>

namespace dlr {

  namespace optimization {

    template <class Functor>
    class GradientFunction
      : public std::unary_function<typename Functor::argument_type,
                                   typename Functor::result_type>
    {
    public:
      GradientFunction(const Functor& functor, double epsilon=1.0e-6) :
        m_functor(functor), m_epsilon(epsilon) {
        if(epsilon == 0.0) {
          DLR_THROW3(ValueException, "GradientFunction::GradientFunction()",
                     "Invalid value (0.0) for argument epsilon.");
        }
      }

      virtual ~GradientFunction() {}

      typename Functor::argument_type
      gradient(const typename Functor::argument_type& theta);

      typename Functor::result_type
      operator()(const typename Functor::argument_type& theta) {
        return m_functor(theta);
      }

    private:
      Functor m_functor;
      double m_epsilon;

    }; // class GradientFunction

  } // namespace optimization

} // namespace dlr


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

namespace dlr {

  using optimization::GradientFunction;

} // namespace dlr


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

namespace dlr {

  namespace optimization {

    // Numerically approximate the gradient of this->operator() by
    // divided differences.
    template <class Functor>
    typename Functor::argument_type
    GradientFunction<Functor>::
    gradient(const typename Functor::argument_type& theta)
    {
      // Create some vectors to use as input to operator()().
      typename Functor::argument_type thetaMinus(theta.size());
      typename Functor::argument_type thetaPlus(theta.size());
      // Return value must be a vector, too, so use argument_type.
      typename Functor::argument_type result(theta.size());
      // Initialize arguments.
      for(size_t index = 0; index < theta.size(); ++index) {
        thetaMinus[index] = theta[index];
        thetaPlus[index] = theta[index];
      }
      // Now compute each partial derivative.
      for(size_t index = 0; index < theta.size(); ++index) {
        // Set up the difference.
        thetaMinus[index] = theta[index] - m_epsilon;
        thetaPlus[index] = theta[index] + m_epsilon;
        // Compute divided difference.
        typename Functor::result_type valueMinus =
          m_functor.operator()(thetaMinus);
        typename Functor::result_type valuePlus =
          m_functor.operator()(thetaPlus);
        result[index] = (valuePlus - valueMinus) / (2 * m_epsilon);
        // Reset argument values.
        thetaMinus[index] = theta[index];
        thetaPlus[index] = theta[index];
      }
      return result;
    }

  } // namespace optimization

} // namespace dlr

#endif // #ifndef _DLR_GRADIENTFUNCTION_H_

---