c2_function.hh

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// $Id: c2_function.hh 97045 2016-05-23 09:49:44Z gcosmo $

#ifndef __has_c2_function_hh
#define __has_c2_function_hh 1

// MSVC does not automatically define numerical constants such as M_PI 
// this came from the msdn website, so it should be right...
#ifdef _MSC_VER
#define _USE_MATH_DEFINES
#define c2_isnan _isnan
#define c2_isfinite _finite
#else
#define c2_isnan std::isnan
#define c2_isfinite std::isfinite
#endif

#include <cmath>
#include <vector>
#include <utility>
#include <string>
#include <stdexcept>
#include <typeinfo>
#include <sstream>
#include <limits> // fails under gcc-4.3 without this here, was ok 

class c2_exception : public std::exception {
public:
    c2_exception(const char msgcode[]) : info(msgcode) { }
    virtual ~c2_exception() throw() { }
    virtual const char* what() const throw() { return info.c_str(); }
private:
    std::string info;
};

// put these forward references here, and with a bogus typename to make swig
template <typename float_type> class c2_composed_function_p;
template <typename float_type> class c2_sum_p;
template <typename float_type> class c2_diff_p;
template <typename float_type> class c2_product_p;
template <typename float_type> class c2_ratio_p;
template <typename float_type> class c2_piecewise_function_p;
template <typename float_type> class c2_quadratic_p;
template <typename float_type> class c2_ptr;
template <typename float_type> class c2_fblock 
{        
public:
        float_type x;
        float_type y;
        float_type yp;
        float_type ypp; 
        bool ypbad;
        bool yppbad; 
};

template <typename float_type=double> class c2_function {
public:    
        const std::string cvs_header_vers() const { return 
                "c2_function.hh 490 2012-04-10 19:05:40Z marcus ";
        }
        
        const std::string cvs_file_vers() const ;
        
public:
  virtual ~c2_function() { 
    if(sampling_grid && !no_overwrite_grid) delete sampling_grid;         
    if(root_info) delete root_info; 
    if(owner_count) {
      std::ostringstream outstr;
      outstr << "attempt to delete an object with non-zero ownership in class";
      outstr << typeid(*this).name() << std::endl;
      //throw c2_exception(outstr.str().c_str());
    }
  }
        
  virtual float_type value_with_derivatives(float_type x, float_type *yprime, 
    float_type *yprime2) const throw(c2_exception) =0 ; 
        
        inline float_type operator () (float_type x) const throw(c2_exception) 
        { return value_with_derivatives(x, (float_type *)0, (float_type *)0); } 

    inline float_type operator () (float_type x, float_type *yprime, 
                            float_type *yprime2) const throw(c2_exception) 
        { return value_with_derivatives(x, yprime, yprime2); } 
        
  float_type find_root(float_type lower_bracket, float_type upper_bracket, 
                       float_type start, 
                       float_type value, int *error=0, 
                       float_type *final_yprime=0, 
                       float_type *final_yprime2=0 ) const throw(c2_exception); 
  float_type partial_integrals(std::vector<float_type> xgrid, 
                               std::vector<float_type> *partials = 0,
                               float_type abs_tol=1e-12, 
                               float_type rel_tol=1e-12, 
                               int derivs=2, bool adapt=true, 
                               bool extrapolate=true) 
    const throw(c2_exception);
        
  float_type integral(float_type amin, float_type amax, 
                      std::vector<float_type> *partials = 0,
                      float_type abs_tol=1e-12, float_type rel_tol=1e-12, 
                      int derivs=2, bool adapt=true, bool extrapolate=true) 
    const throw(c2_exception);

  c2_piecewise_function_p<float_type> *adaptively_sample(
                 float_type amin, float_type amax,
                 float_type abs_tol=1e-12, float_type rel_tol=1e-12,
                 int derivs=2, std::vector<float_type> *xvals=0, 
                 std::vector<float_type> *yvals=0) const throw(c2_exception);
        
  inline float_type xmin() const { return fXMin; }
  inline float_type xmax() const { return fXMax; }
  void set_domain(float_type amin, float_type amax) { fXMin=amin; fXMax=amax; }
                
  size_t get_evaluations() const { return evaluations; }
  void reset_evaluations()  const { evaluations=0; } 
  inline void increment_evaluations() const { evaluations++; }
        
  bool check_monotonicity(const std::vector<float_type> &data, 
                          const char message[]) const throw(c2_exception);
        
  virtual void set_sampling_grid(const std::vector<float_type> &grid) 
    throw(c2_exception); 
        
  std::vector<float_type> *get_sampling_grid_pointer() const { 
    return sampling_grid; } 

  virtual void get_sampling_grid(float_type amin, float_type amax, 
                                 std::vector<float_type> &grid) const ;
        
  void preen_sampling_grid(std::vector<float_type> *result) const;
  void refine_sampling_grid(std::vector<float_type> &grid, 
                            size_t refinement) const;
        
  c2_function<float_type> &normalized_function(
      float_type amin, 
      float_type amax, 
      float_type norm=1.0) const throw(c2_exception);
  c2_function<float_type> &square_normalized_function(
      float_type amin, float_type amax, 
      float_type norm=1.0) 
    const throw(c2_exception);
  c2_function<float_type> &square_normalized_function(
      float_type amin, float_type amax, const c2_function<float_type> &weight, 
      float_type norm=1.0)
    const throw(c2_exception);

  c2_sum_p<float_type> &operator + (const c2_function<float_type> &rhs) const 
  { return *new c2_sum_p<float_type>(*this, rhs); }
  c2_diff_p<float_type> &operator - (const c2_function<float_type> &rhs) const 
  { return *new c2_diff_p<float_type>(*this, rhs); }
  c2_product_p<float_type> &operator * 
  (const c2_function<float_type> &rhs) const 
  { return *new c2_product_p<float_type>(*this, rhs); }
  c2_ratio_p<float_type> &operator / (const c2_function<float_type> &rhs) const 
  { return *new c2_ratio_p<float_type>(*this, rhs); }
  c2_composed_function_p<float_type> & operator ()
    (const c2_function<float_type> &inner) const 
  { return *new c2_composed_function_p<float_type>((*this), inner); }

  float_type get_trouble_point() const { return bad_x_point; }
        
  void claim_ownership() const { owner_count++; }
  size_t release_ownership_for_return() const throw(c2_exception) { 
    if(!owner_count) {
      std::ostringstream outstr;
      outstr << "attempt to release ownership of an unowned function in class ";
      outstr << typeid(*this).name() << std::endl;
      throw c2_exception(outstr.str().c_str());
    }
    owner_count--; 
    return owner_count;
  }
  void release_ownership() const throw(c2_exception) { 
    if(!release_ownership_for_return()) delete this;
  }
  size_t count_owners() const { return owner_count; }

protected:
  c2_function(const c2_function<float_type>  &src) 
  : sampling_grid(0),
    no_overwrite_grid(false),
    fXMin(src.fXMin), fXMax(src.fXMax), root_info(0), owner_count(0)
  {} // copy constructor only copies domain, and is only for internal use
  c2_function() : 
    sampling_grid(0), no_overwrite_grid(0), 
    fXMin(-std::numeric_limits<float_type>::max()), 
    fXMax(std::numeric_limits<float_type>::max()), root_info(0), owner_count(0)
  {}

  virtual void set_sampling_grid_pointer(std::vector<float_type> &grid) 
  { 
    if (sampling_grid && !no_overwrite_grid) delete sampling_grid; 
    sampling_grid=&grid; no_overwrite_grid=1; 
  }
        
  std::vector<float_type> * sampling_grid;
  bool no_overwrite_grid;
        
  float_type fXMin, fXMax;
  mutable size_t evaluations;
  mutable float_type bad_x_point;
public: 
  inline void fill_fblock(c2_fblock<float_type> &fb) const throw(c2_exception) 
  {
    fb.y=value_with_derivatives(fb.x, &fb.yp, &fb.ypp);
    fb.ypbad=c2_isnan(fb.yp) || !c2_isfinite(fb.yp);
    fb.yppbad=c2_isnan(fb.ypp) || !c2_isfinite(fb.ypp);
    increment_evaluations();
  }

private:
  struct recur_item { 
    c2_fblock<float_type> f1; size_t depth; 
    float_type previous_estimate, abs_tol, step_sum; 
    bool done; 
    size_t f0index, f2index;
  };
        

  struct c2_integrate_recur { 
    c2_fblock<float_type> *f0, *f1;
    float_type abs_tol, rel_tol, eps_scale, extrap_coef, extrap2, 
      dx_tolerance, abs_tol_min;
    std::vector< recur_item > *rb_stack;
    int  derivs;
    bool adapt, extrapolate, inited;
  };

  struct c2_sample_recur { 
    c2_fblock<float_type> *f0, *f1;
    float_type abs_tol, rel_tol, dx_tolerance, abs_tol_min;
    int derivs;
    c2_piecewise_function_p<float_type> *out;
    std::vector<float_type> *xvals, *yvals;
    std::vector< recur_item > *rb_stack;
    bool inited;
  };

  struct c2_root_info {
    c2_fblock<float_type> lower, upper;
    bool inited;
  };
        
  float_type integrate_step(struct c2_integrate_recur &rb) 
    const throw(c2_exception);
    
  void sample_step(struct c2_sample_recur &rb) const throw(c2_exception);

  mutable struct  c2_root_info *root_info;
  
  mutable size_t owner_count;
};

template <typename float_type=double> class c2_classic_function_p 
  : public c2_function<float_type> {
public:
  c2_classic_function_p(const float_type (*c_func)(float_type)) 
    : c2_function<float_type>(), func(c_func)  {}

  virtual float_type value_with_derivatives(
    float_type x, float_type *yprime, 
    float_type *yprime2) const throw(c2_exception) 
  {
    if(!func) 
      throw c2_exception("c2_classic_function called with null function");
    if(yprime) *yprime=0;
    if(yprime2) *yprime2=0;
    return func(x); 
  }
  virtual ~c2_classic_function_p() { }
        
protected:
  const float_type (*func)(float_type);
};

template <typename float_type> class c2_const_ptr {
public:
  c2_const_ptr() : func(0)  {}
  c2_const_ptr(const c2_function<float_type> &f) : func(0) 
  { this->set_function(&f); }
  c2_const_ptr(const c2_const_ptr<float_type> &src) : func(0)
  { this->set_function(src.get_ptr()); }
  void set_function(const c2_function<float_type> *f) 
  { 
    if(func) func->release_ownership();
    func=f; 
    if(func) func->claim_ownership();
  }
        
  const c2_const_ptr<float_type> & operator =
  (const c2_const_ptr<float_type> &f) 
  { this->set_function(f.get_ptr()); return f; }
  const c2_function<float_type> & operator =
  (const c2_function<float_type> &f) 
  { this->set_function(&f); return f; }
  void release_for_return() throw(c2_exception)
  {        
    if(func) func->release_ownership_for_return();
    func=0;
  }
  void unset_function(void) { this->set_function(0);  }
  ~c2_const_ptr() { this->set_function(0); }
        
  inline const c2_function<float_type> &get() const throw(c2_exception) 
  { 
    if(!func) throw c2_exception("c2_ptr accessed uninitialized");
    return *func;
  }
  inline const c2_function<float_type> *get_ptr() const { return func; }
  inline const c2_function<float_type> *operator -> () const 
  { return &get(); }
  bool valid() const { return func != 0; }
        
  operator const c2_function<float_type>& () const { return this->get(); }
        

  float_type operator()(float_type x) const throw(c2_exception) 
  { return get()(x); }
  float_type operator()(float_type x, float_type *yprime, 
                        float_type *yprime2) const throw(c2_exception) 
  { return get().value_with_derivatives(x, yprime, yprime2); }
  c2_sum_p<float_type> &operator + (const c2_function<float_type> &rhs) 
    const throw(c2_exception)
  { return *new c2_sum_p<float_type>(get(), rhs); }
  c2_diff_p<float_type> &operator - (const c2_function<float_type> &rhs)  
    const throw(c2_exception)
  { return *new c2_diff_p<float_type>(get(), rhs); }
  c2_product_p<float_type> &operator * (const c2_function<float_type> &rhs) 
    const throw(c2_exception)
  { return *new c2_product_p<float_type>(get(), rhs); }
  c2_ratio_p<float_type> &operator / (const c2_function<float_type> &rhs)  
    const throw(c2_exception)
  { return *new c2_ratio_p<float_type>(get(), rhs); }
  c2_composed_function_p<float_type> & 
  operator ()(const c2_function<float_type> &inner) 
    const throw(c2_exception)
  { return *new c2_composed_function_p<float_type>(get(), inner); }
        
protected:
  const c2_function<float_type> * func;
};

template <typename float_type> class c2_ptr : public c2_const_ptr<float_type >
{
public:
  c2_ptr() : c2_const_ptr<float_type>()  {}
  c2_ptr(c2_function<float_type> &f) : 
    c2_const_ptr<float_type>() { this->set_function(&f); } 
  c2_ptr(const c2_ptr<float_type> &src) : 
    c2_const_ptr<float_type>() { this->set_function(src.get_ptr()); }
  inline c2_function<float_type> &get() const throw(c2_exception) 
  { return *const_cast<c2_function<float_type>*>(
      &c2_const_ptr<float_type>::get()); }
  inline c2_function<float_type> *get_ptr() const 
  { return const_cast<c2_function<float_type>*>(this->func); }
  inline c2_function<float_type> *operator -> () const 
  { return &get(); }
  const c2_ptr<float_type> & operator =(const c2_ptr<float_type> &f) 
  { this->set_function(f.get_ptr()); return f; }
  c2_function<float_type> & operator =(c2_function<float_type> &f) 
  { this->set_function(&f); return f; }
private:
  void operator =(const c2_const_ptr<float_type> &) { }
  void operator =(const c2_function<float_type> &) { }
};

template <typename float_type, template <typename> class c2_class > 
class c2_typed_ptr 
  : public c2_const_ptr<float_type> {
public:
  c2_typed_ptr() : c2_ptr<float_type>()  {}
  c2_typed_ptr(c2_class<float_type> &f)
    : c2_const_ptr<float_type>() { this->set_function(&f); }
  c2_typed_ptr(const c2_typed_ptr<float_type, c2_class> &src) 
    : c2_const_ptr<float_type>() { this->set_function(src.get_ptr()); }
        
  inline c2_class<float_type> &get() const throw(c2_exception) 
  { 
    return *static_cast<c2_class<float_type> *>
      (const_cast<c2_function<float_type>*>(&c2_const_ptr<float_type>::get()));
  }
  inline c2_class<float_type> *operator -> () const 
                { return &get(); }
  inline c2_class<float_type> *get_ptr() const 
  { return static_cast<c2_class<float_type> *>(
       const_cast<c2_function<float_type>*>(this->func)); }

  operator c2_class<float_type>&() const { return get(); }
  void operator =(const c2_typed_ptr<float_type, c2_class> &f) 
  { this->set_function(f.get_ptr()); }
  void operator =(c2_class<float_type> &f) 
  { this->set_function(&f); }
private:
  void operator =(const c2_const_ptr<float_type> &) { }
  void operator =(const c2_function<float_type> &) { }
};

template <typename float_type=double> class c2_plugin_function_p : 
        public c2_function<float_type> {
public:
  c2_plugin_function_p() : c2_function<float_type>(), func()  {}
  c2_plugin_function_p(c2_function<float_type> &f) : 
    c2_function<float_type>(),func(f)  { }

  void set_function(c2_function<float_type> *f) 
  {
    func.set_function(f);
    if(f) this->set_domain(f->xmin(), f->xmax());
  }
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception) 
  {
    if(!func.valid()) 
      throw c2_exception("c2_plugin_function_p called uninitialized");
    return func->value_with_derivatives(x, yprime, yprime2); 
  }
  virtual ~c2_plugin_function_p() { }
        
  void unset_function() { func.unset_function(); }
        
  virtual void get_sampling_grid(float_type amin, float_type amax, 
                                 std::vector<float_type> &grid) const 
  {        
    if(!func.valid()) 
      throw c2_exception("c2_plugin_function_p called uninitialized");
    if(this->sampling_grid) 
      c2_function<float_type>::get_sampling_grid(amin, amax, grid);
    else  func->get_sampling_grid(amin, amax, grid);
  }
protected:
        c2_ptr<float_type> func;
};

template <typename float_type=double> class c2_const_plugin_function_p 
  : public c2_plugin_function_p<float_type> {
public:
  c2_const_plugin_function_p() : c2_plugin_function_p<float_type>()  {}
  c2_const_plugin_function_p(const c2_function<float_type> &f) : 
    c2_plugin_function_p<float_type>()  { this->set_function(&f); }
  void set_function(const c2_function<float_type> *f) 
  { c2_plugin_function_p<float_type>::set_function(
       const_cast<c2_function<float_type>*>(f)); }
  virtual ~c2_const_plugin_function_p() { }
        
  const c2_function<float_type> &get() const throw(c2_exception) 
  { return this->func.get(); }
};

template <typename float_type=double> class c2_binary_function 
  : public c2_function<float_type> {
public:
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw (c2_exception) 
  { 
    if(stub) 
      throw c2_exception("attempt to evaluate a c2_binary_function stub");
    return this->combine(*Left.get_ptr(), *Right.get_ptr(), x, yprime, yprime2);
  }

  virtual ~c2_binary_function() { }

protected:
  c2_binary_function( float_type (*combiner)(
     const c2_function<float_type> &left, 
     const c2_function<float_type> &right, 
     float_type x, float_type *yprime, 
     float_type *yprime2),
                      const c2_function<float_type> &left,  
                      const c2_function<float_type> &right) : 
    c2_function<float_type>(), combine(combiner), Left(left), 
    Right(right), stub(false)
  { 
    this->set_domain(
                     (left.xmin() > right.xmin()) ? left.xmin() : right.xmin(), 
                     (left.xmax() < right.xmax()) ? left.xmax() : right.xmax()
                     );
  } 
  c2_binary_function(float_type (*combiner)(
      const c2_function<float_type> &left, 
      const c2_function<float_type> &right, 
      float_type x, float_type *yprime, float_type *yprime2)
                     ) : c2_function<float_type>(), combine(combiner), 
    Left(), Right(), stub(true) { }
        
public:
  float_type (* const combine)(
      const c2_function<float_type> &left, 
      const c2_function<float_type> &right, 
      float_type x, float_type *yprime, float_type *yprime2);
        
protected:        
  const c2_const_ptr<float_type> Left,  Right;
  bool stub;
};

template <typename float_type=double> class c2_scaled_function_p 
  : public c2_function<float_type> {
public:
  c2_scaled_function_p(const c2_function<float_type> &outer, 
                       float_type scale) : 
    c2_function<float_type>(), func(outer), yscale(scale) { }
        
  void reset(float_type scale) { yscale=scale; }
        
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw (c2_exception) 
  { 
    float_type y=this->func->value_with_derivatives(x, yprime, yprime2); 
    if(yprime) (*yprime)*=yscale; 
    if(yprime2) (*yprime2)*=yscale; 
    return y*yscale; 
  }

protected:
  c2_scaled_function_p<float_type>() : func() {} 
  const c2_const_ptr<float_type> func;
  float_type yscale;
};

template <typename float_type=double> class c2_cached_function_p 
  : public c2_function<float_type> {
public:
  c2_cached_function_p(const c2_function<float_type> &f) 
  : c2_function<float_type>(),
    func(f), init(false)  {}
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception) 
  {
    if(!init || x != x0) {
      y=this->func->value_with_derivatives(x, &yp, &ypp);
      x0=x;
      init=true;
    }
    if(yprime) *yprime=yp;
    if(yprime2) *yprime2=ypp;
    return y; 
  }

protected:
  c2_cached_function_p() : func() {}
  const c2_const_ptr<float_type> func;
  mutable bool init;
  mutable float_type x0, y, yp, ypp;
        
};

template <typename float_type=double> class  c2_composed_function_p 
  : public c2_binary_function<float_type> {
public:
  c2_composed_function_p(const c2_function<float_type> &outer, 
                         const c2_function<float_type> &inner) : 
    c2_binary_function<float_type>(combine, outer, inner) { 
    this->set_domain(inner.xmin(), inner.xmax()); }
  c2_composed_function_p() : c2_binary_function<float_type>(combine) {} 

  static float_type combine(const c2_function<float_type> &left, 
                            const c2_function<float_type> &right, 
                            float_type x, float_type *yprime, 
                            float_type *yprime2) throw(c2_exception)
  {
    float_type y0, y1;
    if(yprime || yprime2) {
      float_type yp0, ypp0, yp1, ypp1;
      y0=right.value_with_derivatives(x, &yp0, &ypp0);
      y1=left.value_with_derivatives(y0, &yp1, &ypp1);
      if(yprime) *yprime=yp1*yp0;
      if(yprime2) *yprime2=ypp0*yp1+yp0*yp0*ypp1;
    } else {
      y0=right(x);
      y1=left(y0);
    }
    return y1;
  }        
};

template <typename float_type=double> class c2_sum_p 
  : public c2_binary_function<float_type> {
public:        
  c2_sum_p(const c2_function<float_type> &left, 
           const c2_function<float_type> &right) 
    : c2_binary_function<float_type>(combine, left, right) {}
  c2_sum_p() : c2_binary_function<float_type>(combine) {} ; 

  static float_type combine(const c2_function<float_type> &left, 
                            const c2_function<float_type> &right, 
                            float_type x, float_type *yprime, 
                            float_type *yprime2) throw(c2_exception)
  {
    float_type y0, y1;
    if(yprime || yprime2) {
      float_type yp0, ypp0, yp1, ypp1;
      y0=left.value_with_derivatives(x, &yp0, &ypp0);
      y1=right.value_with_derivatives(x, &yp1, &ypp1);
      if(yprime) *yprime=yp0+yp1;
      if(yprime2) *yprime2=ypp0+ypp1;
    } else {
      y0=left(x);
      y1=right(x);
    }
    return y0+y1;
  }
};

template <typename float_type=double> class c2_diff_p 
  : public c2_binary_function<float_type> {
public:        
  c2_diff_p(const c2_function<float_type> &left, 
            const c2_function<float_type> &right) 
    : c2_binary_function<float_type>(combine, left, right) {}
  c2_diff_p() : c2_binary_function<float_type>(combine) {} ; 

  static float_type combine(const c2_function<float_type> &left, 
                            const c2_function<float_type> &right, 
                            float_type x, float_type *yprime, 
                            float_type *yprime2) throw(c2_exception)
  {
    float_type y0, y1;
    if(yprime || yprime2) {
      float_type yp0, ypp0, yp1, ypp1;
      y0=left.value_with_derivatives(x, &yp0, &ypp0);
      y1=right.value_with_derivatives(x, &yp1, &ypp1);
      if(yprime) *yprime=yp0-yp1;
      if(yprime2) *yprime2=ypp0-ypp1;
    } else {
      y0=left(x);
      y1=right(x);
    }
    return y0-y1;
  }
};


template <typename float_type=double> class c2_product_p 
  : public c2_binary_function<float_type> {
public:        
  c2_product_p(const c2_function<float_type> &left, 
               const c2_function<float_type> &right) 
    : c2_binary_function<float_type>(combine, left, right) {}
  c2_product_p() : c2_binary_function<float_type>(combine) {} ; 
  
  static float_type combine(const c2_function<float_type> &left, 
                            const c2_function<float_type> &right, 
                            float_type x, float_type *yprime, 
                            float_type *yprime2)  throw(c2_exception)
    {
      float_type y0, y1;
      if(yprime || yprime2) {
        float_type yp0, ypp0, yp1, ypp1;
        y0=left.value_with_derivatives(x, &yp0, &ypp0);
        y1=right.value_with_derivatives(x, &yp1, &ypp1);
        if(yprime) *yprime=y1*yp0+y0*yp1;
        if(yprime2) *yprime2=ypp0*y1+2.0*yp0*yp1+ypp1*y0;
      } else {
        y0=left(x);
        y1=right(x);
      }
      return y0*y1;
    }
};


template <typename float_type=double> class c2_ratio_p 
  : public c2_binary_function<float_type> {
public:        
  c2_ratio_p(const c2_function<float_type> &left, 
             const c2_function<float_type> &right) 
    : c2_binary_function<float_type>(combine, left, right) {}
  c2_ratio_p() : c2_binary_function<float_type>(combine) {} ; 
        
  static float_type combine(const c2_function<float_type> &left, 
                            const c2_function<float_type> &right, 
                            float_type x, float_type *yprime, 
                            float_type *yprime2) throw(c2_exception)
  {
    float_type y0, y1;
    if(yprime || yprime2) {
      float_type yp0, ypp0, yp1, ypp1;
      y0=left.value_with_derivatives(x, &yp0, &ypp0);
      y1=right.value_with_derivatives(x, &yp1, &ypp1);
      if(yprime) *yprime=(yp0*y1-y0*yp1)/(y1*y1); // first deriv of ratio
      if(yprime2) *yprime2=(y1*y1*ypp0+y0*(2*yp1*yp1-y1*ypp1)-2*y1*yp0*yp1)
                    /(y1*y1*y1); 
    } else {
      y0=left(x);
      y1=right(x);
    }
    return y0/y1;
  }
};

template <typename float_type> class c2_constant_p 
  : public c2_function<float_type> {
public:
  c2_constant_p(float_type x) : c2_function<float_type>(), value(x) {}
  void reset(float_type val) { value=val; }
  virtual float_type value_with_derivatives(
          float_type, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception) 
  { if(yprime) *yprime=0; if(yprime2) *yprime2=0; return value; }
        
private:
  float_type value;
};

template <typename float_type> class c2_transformation {
public:
  c2_transformation(bool transformed,
                    float_type (*xin)(float_type), 
                    float_type (*xinp)(float_type), 
                    float_type (*xinpp)(float_type), 
                    float_type (*xout)(float_type)
                    ) :
    fTransformed(transformed), fHasStaticTransforms(true),
    pIn(xin), pInPrime(xinp), pInDPrime(xinpp), pOut(xout) { }

  c2_transformation(bool transformed) :
    fTransformed(transformed), fHasStaticTransforms(false),
    pIn(report_error), pInPrime(report_error), pInDPrime(report_error), 
    pOut(report_error) { }
  virtual ~c2_transformation() { }
  const bool fTransformed;
  const bool fHasStaticTransforms;

  float_type (* const pIn)(float_type);
  float_type (* const pInPrime)(float_type);
  float_type (* const pInDPrime)(float_type);
  float_type (* const pOut)(float_type);
        
  virtual float_type fIn(float_type x) const { return pIn(x); }
  virtual float_type fInPrime(float_type x) const { return pInPrime(x); }
  virtual float_type fInDPrime(float_type x) const { return pInDPrime(x); }
  virtual float_type fOut(float_type x) const { return pOut(x); }

protected:
  static float_type report_error(float_type x)  { 
    throw c2_exception("use of improperly constructed axis transform"); 
    return x; }
  static float_type ident(float_type x)  { return x; }
  static float_type one(float_type)  { return 1; }
  static float_type zero(float_type)  { return 0; }
  static float_type recip(float_type x)  { return 1.0/x; }
  static float_type recip_prime(float_type x)  { return -1/(x*x); }
  static float_type recip_prime2(float_type x)  { return 2/(x*x*x); }
};

template <typename float_type> class c2_transformation_linear 
  : public c2_transformation<float_type> {
public:
  c2_transformation_linear() : c2_transformation<float_type>(
     false, this->ident, 
     this->one, this->zero, 
     this->ident) { }
  virtual ~c2_transformation_linear() { }
};
template <typename float_type> class c2_transformation_log 
  : public c2_transformation<float_type> {
public:
  c2_transformation_log() : c2_transformation<float_type>(
     true, std::log, this->recip, 
     this->recip_prime, std::exp) { }
  virtual ~c2_transformation_log() { }
};
template <typename float_type> class c2_transformation_recip 
  : public c2_transformation<float_type> {
public:
  c2_transformation_recip() : c2_transformation<float_type>(
     true, this->recip, 
     this->recip_prime, 
     this->recip_prime2, this->recip) { }
  virtual ~c2_transformation_recip() { }
};

template <typename float_type> 
class c2_function_transformation {
public:
  c2_function_transformation(
     const c2_transformation<float_type> &xx, 
     const c2_transformation<float_type> &yy) :
    isIdentity(!(xx.fTransformed || yy.fTransformed)), X(xx), Y(yy) { }
  virtual ~c2_function_transformation() { delete &X; delete &Y; }
  virtual float_type evaluate(
          float_type xraw, 
          float_type y, float_type yp0, float_type ypp0,
          float_type *yprime, float_type *yprime2) const;
  const bool isIdentity;
  const c2_transformation<float_type> &X;
  const c2_transformation<float_type> &Y;
};

template <typename float_type> class c2_lin_lin_function_transformation : 
  public c2_function_transformation<float_type> {
public:
  c2_lin_lin_function_transformation() : 
    c2_function_transformation<float_type>(
       *new c2_transformation_linear<float_type>, 
       *new c2_transformation_linear<float_type>
                                           ) { }
  virtual ~c2_lin_lin_function_transformation() { }
};

template <typename float_type> class c2_log_log_function_transformation : 
        public c2_function_transformation<float_type> {
public:
        c2_log_log_function_transformation() : 
                c2_function_transformation<float_type>(
                        *new c2_transformation_log<float_type>, 
                        *new c2_transformation_log<float_type>
                ) { }
        virtual ~c2_log_log_function_transformation() { }
};

template <typename float_type> class c2_lin_log_function_transformation : 
        public c2_function_transformation<float_type> {
public:
        c2_lin_log_function_transformation() : 
                c2_function_transformation<float_type>(
                        *new c2_transformation_linear<float_type>, 
                        *new c2_transformation_log<float_type>
                ) { }
        virtual ~c2_lin_log_function_transformation() { }
};

template <typename float_type> class c2_log_lin_function_transformation : 
        public c2_function_transformation<float_type> {
public:
        c2_log_lin_function_transformation() : 
                c2_function_transformation<float_type>(
                        *new c2_transformation_log<float_type>, 
                        *new c2_transformation_linear<float_type>
                ) { }
        virtual ~c2_log_lin_function_transformation() { }
};

template <typename float_type> class c2_arrhenius_function_transformation : 
        public c2_function_transformation<float_type> {
public:
        c2_arrhenius_function_transformation() : 
                c2_function_transformation<float_type>(
                        *new c2_transformation_recip<float_type>, 
                        *new c2_transformation_log<float_type>
                ) { }
        virtual ~c2_arrhenius_function_transformation() { }
};

template <typename float_type=double> class interpolating_function_p  
  : public c2_function<float_type> {
public:
  interpolating_function_p() : c2_function<float_type>(), 
    fTransform(*new  c2_lin_lin_function_transformation<float_type>) { } 

  interpolating_function_p(const c2_function_transformation<float_type> &
                           transform) 
    : c2_function<float_type>(), 
      fTransform(transform) { } 

  interpolating_function_p<float_type> & load(
      const std::vector<float_type> &x, 
      const std::vector<float_type> &f, 
      bool lowerSlopeNatural, float_type lowerSlope, 
      bool upperSlopeNatural, float_type upperSlope, bool splined=true
        ) throw(c2_exception);

  interpolating_function_p<float_type> & load_pairs(
     std::vector<std::pair<float_type, float_type> > &data,  
     bool lowerSlopeNatural, float_type lowerSlope, 
     bool upperSlopeNatural, float_type upperSlope, bool splined=true
                                                    ) throw(c2_exception);

  interpolating_function_p<float_type> &  
  sample_function(const c2_function<float_type> &func, 
                  float_type amin, float_type amax, 
                  float_type abs_tol, float_type rel_tol,
                  bool lowerSlopeNatural, float_type lowerSlope, 
                  bool upperSlopeNatural, float_type upperSlope
                  ) throw(c2_exception);
        
  interpolating_function_p<float_type> & load_random_generator_function(
     const std::vector<float_type> &bincenters, 
     const c2_function<float_type> &binheights)
    throw(c2_exception);

  interpolating_function_p<float_type> & load_random_generator_bins(
    const std::vector<float_type> &bins, 
    const std::vector<float_type> &binheights, 
    bool splined=true)
    throw(c2_exception);
        
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception);
        
  virtual ~interpolating_function_p() { delete &fTransform; } 
        
  virtual interpolating_function_p<float_type> &clone() 
    const throw(c2_exception)
  {        return *new interpolating_function_p<float_type>(); }
        
  void get_data(std::vector<float_type> &xvals, 
                std::vector<float_type> &yvals) const throw() ;
        
  void get_internal_data(
       std::vector<float_type> &xvals, 
       std::vector<float_type> &yvals, 
       std::vector<float_type> &y2vals) const 
  { xvals=X; yvals=F; y2vals=y2; }

  void set_lower_extrapolation(float_type bound);
  void set_upper_extrapolation(float_type bound);
        
  interpolating_function_p <float_type> & 
  unary_operator(const c2_function<float_type> &source) const;

  interpolating_function_p <float_type> & 
  binary_operator(const c2_function<float_type> &rhs,
                  const c2_binary_function<float_type> *combining_stub
                  ) const;
  interpolating_function_p <float_type> & 
  add_pointwise (const c2_function<float_type> &rhs) const { 
    return binary_operator(rhs, new c2_sum_p<float_type>()); }
  interpolating_function_p <float_type> & 
  subtract_pointwise (const c2_function<float_type> &rhs) const {
                return binary_operator(rhs, new c2_diff_p<float_type>()); }
  interpolating_function_p <float_type> & 
  multiply_pointwise (const c2_function<float_type> &rhs) const { 
                return binary_operator(rhs, new c2_product_p<float_type>()); }
  interpolating_function_p <float_type> & 
  divide_pointwise (const c2_function<float_type> &rhs) const { 
                return binary_operator(rhs, new c2_ratio_p<float_type>()); }
  void clone_data(const interpolating_function_p <float_type> &rhs) {
    Xraw=rhs.Xraw; X=rhs.X; F=rhs.F; y2=rhs.y2;
    set_sampling_grid_pointer(Xraw);
  }

  const c2_function_transformation<float_type> &fTransform;
        
protected:    
  void spline(
              bool lowerSlopeNatural, float_type lowerSlope, 
              bool upperSlopeNatural, float_type upperSlope
              ) throw(c2_exception);

  static bool comp_pair(std::pair<float_type,float_type> const &i, 
                        std::pair<float_type,float_type> const &j) 
  {return i.first<j.first;}

  std::vector<float_type> Xraw, X, F, y2;
  c2_const_ptr<float_type> sampler_function;
  bool xInverted;
  mutable size_t lastKLow;
};

template <typename float_type=double> class log_lin_interpolating_function_p 
  : public interpolating_function_p <float_type> {
public:
  log_lin_interpolating_function_p() : 
    interpolating_function_p<float_type>(
       *new c2_log_lin_function_transformation<float_type>)
  { }
  virtual interpolating_function_p<float_type> &clone() 
    const throw(c2_exception)
  { return *new log_lin_interpolating_function_p<float_type>(); }
};


template <typename float_type=double> class lin_log_interpolating_function_p 
  : public interpolating_function_p <float_type> {
public:
  lin_log_interpolating_function_p() 
  : interpolating_function_p<float_type>(
      *new c2_lin_log_function_transformation<float_type>)
  { } 
  virtual interpolating_function_p<float_type> &clone() 
    const throw(c2_exception)
  { return *new lin_log_interpolating_function_p<float_type>(); }
};


template <typename float_type=double> class log_log_interpolating_function_p 
  : public interpolating_function_p <float_type> {
public:        
  log_log_interpolating_function_p() : 
    interpolating_function_p<float_type>(
      *new c2_log_log_function_transformation<float_type>)
  { } 
  virtual interpolating_function_p<float_type> &clone() 
    const throw(c2_exception)
  { return *new log_log_interpolating_function_p<float_type>(); }
};


template <typename float_type=double> class arrhenius_interpolating_function_p 
  : public interpolating_function_p <float_type> {
public:
  arrhenius_interpolating_function_p() 
  : interpolating_function_p<float_type>(
      *new c2_arrhenius_function_transformation<float_type>)
  { } 
  virtual interpolating_function_p<float_type> &clone() 
    const throw(c2_exception)
  { return *new arrhenius_interpolating_function_p<float_type>(); }        
};

template <typename float_type=double> class c2_sin_p 
  : public c2_function<float_type> {
public:
  c2_sin_p() : c2_function<float_type>() {}
        
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception)
  { float_type q=std::sin(x); 
    if(yprime) *yprime=std::cos(x); 
    if(yprime2) *yprime2=-q; 
    return q; }
    
  virtual void get_sampling_grid(float_type amin, float_type amax,  
                                 std::vector<float_type> &grid) const; 
};

template <typename float_type=double> class c2_cos_p 
  : public c2_sin_p<float_type> {
public:
  c2_cos_p() : c2_sin_p<float_type>() {}
        
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception)
  { float_type q=std::cos(x); 
    if(yprime) *yprime=-std::sin(x); 
    if(yprime2) *yprime2=-q; 
    return q; }
};

template <typename float_type=double> class c2_tan_p 
  : public c2_function<float_type> {          
public:
  c2_tan_p() : c2_function<float_type>() {}
        
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception)
  {
    float_type c=std::cos(x), ss=std::sin(x);
    float_type t=ss/c;
    float_type yp=1/(c*c);
    if(yprime) { *yprime=yp; } 
    if(yprime2){*yprime2=2*t*yp; } 
    return t; 
  }        
};

template <typename float_type=double> class c2_log_p 
  : public c2_function<float_type> {
public:
  c2_log_p() : c2_function<float_type>() {}

  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception)
  { if(yprime) *yprime=1.0/x; 
    if(yprime2) *yprime2=-1.0/(x*x); 
    return std::log(x); }        
};

template <typename float_type=double>  class c2_exp_p 
  : public c2_function<float_type> {
public:
  c2_exp_p() : c2_function<float_type>() {}

  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception)
  { float_type q=std::exp(x); 
    if(yprime) *yprime=q; 
    if(yprime2) *yprime2=q; 
    return q; }
};

template <typename float_type=double> class c2_sqrt_p 
  : public c2_function<float_type> {
public:
  c2_sqrt_p() : c2_function<float_type>() {}
        
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception)
  { float_type q=std::sqrt(x); 
    if(yprime) *yprime=0.5/q; 
    if(yprime2) *yprime2=-0.25/(x*q); 
    return q; }
};

template <typename float_type=double> class c2_recip_p 
  : public c2_function<float_type> {
public:
  c2_recip_p(float_type scale) : c2_function<float_type>(), rscale(scale) {}
        
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception)
  { 
    float_type q=1.0/x; 
    float_type y=rscale*q;
    if(yprime) *yprime=-y*q; 
    if(yprime2) *yprime2=2*y*q*q; 
    return y; 
  }
  void reset(float_type scale) { rscale=scale; } 
private:
  float_type rscale;
};

template <typename float_type=double> class c2_identity_p 
  : public c2_function<float_type> {
public:
  c2_identity_p() : c2_function<float_type>() {}
        
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception)
  { if(yprime) *yprime=1.0; if(yprime2) *yprime2=0; return x; }
};

template <typename float_type=double> class c2_linear_p 
  : public c2_function<float_type> {
public:
  c2_linear_p(float_type x0, float_type y0, float_type slope) : 
    c2_function<float_type>(), xint(x0), intercept(y0), m(slope) {}
  void reset(float_type x0, float_type y0, float_type slope) 
  { xint=x0; intercept=y0; m=slope; } 
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception)
  { if(yprime) *yprime=m; 
    if(yprime2) *yprime2=0; 
    return m*(x-xint)+intercept; }
        
private:
  float_type xint, intercept, m;
protected:
  c2_linear_p() {} 
};

template <typename float_type=double> class c2_quadratic_p 
  : public c2_function<float_type> {
public:
  c2_quadratic_p(float_type x0, float_type y0, 
                 float_type xcoef, float_type x2coef) : 
    c2_function<float_type>(), intercept(y0), center(x0), a(x2coef), b(xcoef) {}
  void reset(float_type x0, float_type y0, float_type xcoef,
             float_type x2coef) { intercept=y0; center=x0; a=x2coef; b=xcoef; } 
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception)
  { float_type dx=x-center; 
    if(yprime) *yprime=2*a*dx+b; 
    if(yprime2) *yprime2=2*a; 
    return a*dx*dx+b*dx+intercept; }
        
private:
  float_type intercept, center, a, b;
protected:
  c2_quadratic_p() {} 
};

template <typename float_type=double> class c2_power_law_p 
  : public c2_function<float_type> {
public:
  c2_power_law_p(float_type scale, float_type power) :
    c2_function<float_type>(), a(scale), b(power) {}
  void reset(float_type scale, float_type power) { a=scale; b=power; }
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception)
  { float_type q=a*std::pow(x,b-2); 
    if(yprime) *yprime=b*q*x; 
    if(yprime2) *yprime2=b*(b-1)*q; 
    return q*x*x; }
        
private:
  float_type a, b;
protected:
  c2_power_law_p() {} 
};

template <typename float_type=double> class c2_inverse_function_p 
  : public c2_function<float_type> {
public:
  c2_inverse_function_p(const c2_function<float_type> &source);  
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw(c2_exception);
    
    void set_start_hint(float_type hint) const { start_hint=hint; }
    
    virtual float_type get_start_hint(float_type x) const 
  { return hinting_function.valid()? hinting_function(x) : start_hint; } 

  void set_hinting_function(const c2_function<float_type> *hint_func) 
  { hinting_function.set_function(hint_func); }
  void set_hinting_function(const c2_const_ptr<float_type> hint_func) 
  { hinting_function=hint_func; }
        
protected:
  c2_inverse_function_p() {} 
  mutable float_type start_hint;
  const c2_const_ptr<float_type> func;
  c2_const_ptr<float_type> hinting_function;
};

template <typename float_type=double>  class accumulated_histogram 
  : public interpolating_function_p <float_type> {
public:
  accumulated_histogram(const std::vector<float_type>binedges, 
                        const std::vector<float_type> binheights,
                        bool normalize=false, 
                        bool inverse_function=false, bool drop_zeros=true);
        
};

template <typename float_type=double> class c2_connector_function_p 
  : public c2_function<float_type> {
public:        
  c2_connector_function_p(float_type x0, const c2_function<float_type> &f0, 
                          float_type x2, const c2_function<float_type> &f2, 
                        bool auto_center, float_type y1);
  c2_connector_function_p(
                float_type x0, float_type y0, float_type yp0, float_type ypp0, 
                float_type x2, float_type y2, float_type yp2, float_type ypp2, 
                bool auto_center, float_type y1);
  c2_connector_function_p(
                const c2_fblock<float_type> &fb0, 
                const c2_fblock<float_type> &fb2, 
                bool auto_center, float_type y1);

  virtual ~c2_connector_function_p();
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw (c2_exception);
protected:
        void init(
                  const c2_fblock<float_type> &fb0, 
                  const c2_fblock<float_type> &fb2, 
                  bool auto_center, float_type y1);

  float_type fhinv, fy1, fa, fb, fc, fd, fe, ff;
};

template <typename float_type=double> class c2_piecewise_function_p 
  : public c2_function<float_type> {
public:        
  c2_piecewise_function_p();
  virtual ~c2_piecewise_function_p();
  virtual float_type value_with_derivatives(
          float_type x, float_type *yprime, 
          float_type *yprime2) const throw (c2_exception);
  void append_function(const c2_function<float_type> &func) 
    throw (c2_exception);
protected:
  std::vector<c2_const_ptr<float_type> > functions;
  mutable int lastKLow;
};

#include "c2_function.icc"

#endif