Evaluable.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2009 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
// psrchive/More/MEAL/MEAL/Evaluable.h
 
#ifndef __MEAL_Evaluable_H
#define __MEAL_Evaluable_H
 
#include "MEAL/Function.h"
#include "MEAL/EvaluationPolicy.h"
 
#include "complex_math.h"
#include "Estimate.h"

template<typename T>
struct EstimateTraits {
  typedef Estimate<T> type;
};
 
template<typename T>
void add_variance (typename EstimateTraits<T>::type& result,
                   double var, const T& grad)
{
  result.var += var * grad * grad;
}
 
namespace MEAL {

  template<typename T>
  class Evaluable : public Function {
 
  public:

    typedef T Result;

    static const char* Name;

    Evaluable ()
    {
      evaluation_policy = default_evaluation_policy (this);
    }

    Evaluable (const Evaluable& copy) : Function (copy)
    {
      evaluation_policy = default_evaluation_policy (this);
    }

    Evaluable& operator = (const Evaluable& copy)
    {
      Function::operator = (copy);
      return *this;
    }

    T evaluate (std::vector<T>* grad=0) const
    { 
      T result = evaluation_policy->evaluate (grad);
      if (!true_math::finite(result))
        throw Error (InvalidState, "MEAL::Evaluable::evaluate for " + this->get_name(),
                     "non-finite result");
      return result;
    }

    typename EstimateTraits<T>::type estimate () const
    {
      std::vector< T > gradient;
      typename EstimateTraits<T>::type result = evaluate (&gradient);
 
      const unsigned nparam = get_nparam();
 
      if (gradient.size() != nparam)
        throw Error (InvalidState, "MEAL::Evaluable::evaluate",
                     "gradient.size=%d != nparam=%d", gradient.size(), nparam);
 
      for (unsigned iparam=0; iparam<nparam; iparam++)
        add_variance( result, get_variance(iparam), gradient[iparam]);
 
      return result;
    }
 
  protected:
 
    template<class P> friend class CalculatePolicy;

    virtual void calculate (T& result, std::vector<T>*) = 0;

    Reference::To< EvaluationPolicy< T > > evaluation_policy;
 
    void copy_evaluation_policy (const Evaluable* other)
    {
      evaluation_policy = other->evaluation_policy;
    }

    void calculate (Evaluable* eval, T& result, std::vector<T>* grad)
    {
      eval->calculate (result, grad);
    }
 
  };
 
}
 
template<typename T>
const char* MEAL::Evaluable<T>::Name = "Evaluable<T>";
 
#endif