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 66995 2013-01-29 14:46:45Z gcosmo $

#ifndef __has_c2_function_hh
#define __has_c2_function_hh 1

// MSVC does not automatically define numerical constants such as M_PI without this.
// 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 in c2_function.cc before

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 happy.
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 ; // { return 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) ; // solve f(x)=value

        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; } // evaluations are 'invisible' to constant
        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)
                {} // prevent accidental naked construction (impossible any since this has pure virtual methods)
        
        // this should only be called very early on, by a constructor, before anyone else
        // sets a sampling grid, or it will leak memory
        virtual void set_sampling_grid_pointer(std::vector<float_type> &grid) 
                { 
                        if (sampling_grid && !no_overwrite_grid) delete sampling_grid; // grid was ours, lose it. 
                        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() {} // hide default constructor, since its use is almost always an error.
        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() {} // hide default constructor, since its use is almost always an error.
        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) {} ; // create a stub just for the combiner to avoid statics

        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) {} ; // create a stub just for the combiner to avoid statics

        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) {} ; // create a stub just for the combiner to avoid statics

        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) {} ; // create a stub just for the combiner to avoid statics
        
        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);
        
        // these functions correctly combine the interpolating function with another interpolating function 
        // preserving the X bounds and mapping functions of the host (left hand) function.
        
        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);

        // This is for sorting the data. It must be static if it's going to be a class member.
        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() {} // do not allow naked construction... it is usually an accident.
};

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() {} // do not allow naked construction... it is usually an accident.
};

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() {} // do not allow naked construction... it is usually an accident.
};

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() {} // do not allow naked construction... it is usually an accident.
        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