Re: [Boost-users] Extending boost::function for mathematics?

Subject: Re: [Boost-users] Extending boost::function for mathematics?
From: Steven Watanabe (watanabesj_at_[hidden])
Date: 2009-01-24 20:08:51

AMDG

Jesse Perla wrote:
> Steven: beautifully simple code on your previous post.  I just fear 
> that MPL is beyond the reach of mortals such as myself because I could 
> never come up with that solution.
>
> If people are in the mood for this kind of stuff, let me post on a set 
> of questions coming up in my libraries:
> 1) QUESTION 1: MANAGING DERIVATIVES FOR FUNCTIONS
> * I am doing a lot of work with functions that may have an associated 
> derivative analytically or numerically determined.
> * This is the tip of the iceberg in associating a bunch of information 
> with a mathematical function that all needs to be bundled.
> * I have been creating functors like the following:
> template <class ParamType>
> struct square{
>   typedef double result_type; //For adaptors
>   double operator()(double x, const ParamType& params ) { return x*x;}
>   double gradient(double x, const ParamType& params ){return 2 * x;}
>   boost::function<double (double,const ParamType& 
> params)>  gradient(){return boost::bind(&square::gradient, *this, _1, 
> _2)}; //I forget the exact notation for bind with member functions, 
> but the idea is that it would return a boost::function.
> }; //Analytical derivative
>
> But I would also love the ability to use finitedifference or 
> autodifferentiation algorithms for more complicated examples.  Perhaps 
> something like:
>
> struct square : finite_differences<int Order>{
>     double operator()(double x, const ParamType& params ) { return x*x;}
> }; //Analytical derivative
>
> * Ideally here, finite_differences<int Order> would be able to 
> generate the derivative using order taylor series approx and the 
> operator() to evaluate the function.  Even here I am not sure the best 
> way to organize it.
>
> 2) QUESTION 2: ADAPTATION AND USE LIKE BOOST::FUNCTION
> Now, lets say that there was something similar to boost::function but 
> which had the "derivative" concept to call.gradient(...) to evaluate 
> grad, or return a functor which is the gradient as another function.. 
>  Lets call it math_function for now.  The behavior is things like the 
> following:
> using lambda::_1;
> using lambda::_2;
> math_function<double (double)> mysquare(_1 * _1); //This will generate 
> a math_function with the finite differences as its derivative.
> assert(mysquare(1) == 1);
> math_function<double (double)> mysquaregrad = mysquare.gradient(); 
> //Gets the derivative, which was generated through finite differences. 
>  This could then be used to get the .gradient again for example.
> assert(mysquare.gradient(.1) == mysquaregrad(.1));
>
> math_function<double (double, int)> myfunc(_1 + _2);
> math_function<double (double)> myboundfunc = bind(my_func, _1, 3);
> assert(myboundfunc.gradient(1.0) == 4.0);
>
> //And then to create one with an analytical derivative
> math_function<double (double)> mysquare(_1 * _1, 2 * _1);
> assert(mysquare.gradient(1.0) == 2.0); //evaluated through the 
> analytical derivative.
> boost::function<double (double)> myboostsquare = 
> mysquare.as_function(); //would be nice for interoperability... or 
> something similar.
>
> //Then we can create a class such as polynomial which subclasses 
> math_function

It's probably not a good idea to use inheritance.  Here's an adaptation of
the standard counter example:

polynomial_function<3> mypoly(2, 3, 5);
math_function<double(double)>& mypoly2(mypoly); // note that this is a 
reference.
mypoly2 = some function that is not a polynomial;  // Now what?  mypoly2 
and mypoly are the same object.

> polynomial_function<3> mypoly({{2, 3, 5}}); //Creates a third order 
> polynomial.  The operator() uses fancy algorithms for polynomial 
> evaluation.
> math_function<double (double)> mypoly2 = mypoly; //Can adapt to a 
> math_function for use in generic algorithms.
> math_function<double (double)> mypolygrad1 = mypoly.gradient();
> polynomial_function<2> mypolygrad2 = mypoly.gradient(); //Maybe this 
> is getting too fancy...

Provided that you define polynomial_function appropriately this is no 
problem.
I would probably rely on static polymorphism and type erasure to make
polynomial_function optimally efficient and at the same time allow it
to be stored in a math_function<double(double)>

> 3) QUESTION 3: WHY DO I WANT ALL THIS STUFF?
> For the same reason that boost::function is so useful for functional 
> programming. I am going to have a million permutations on this kind of 
> pattern, am writing very generic algorithms, and would prefer not to 
> have to write a superclass for every permutation of parameters.  I 
> want to be able to freely bind to it, easily use lambdas, and 
> eventually want to add a bunch of other traits such as # of 
> derivatives, range/domain, continuity, hessians, etc. that can be 
> uniformly managed.

You would need to capture all this information somehow, so would you
want a single math_function<double(double)> type to support all of it, or
would you need to be able to say something like:

math_function<double(double),
    // supported features
    boost::mpl::vector<
        gradient,
        n_derivative,
        hessian
    >
 > f;

Obviously the former is much easier to implement.

> 4) QUESTION 4: AM I INSANE?
> Is this at all possible or has something similar been done?  It seems 
> to me that it might be, but I am likely not be capable of doing it. 
>  Is piggybacking boost::function in some smart way here the best 
> approach?  Even if we can't subclass directly, can we copy the source 
> code and intelligently add in a few extra functions, etc.?

So, for an arbitrary function, you want the derivative to be calculated
using finite_difference, right?  and for polynomial_function, you want 
to use
the straightforward derivative instead of the generic version.

So, suppose that we give the gradient a default implementation
like this:

template<class T>
struct gradient {
    typedef finite_differences<T> result_type;
    finite_differences<T> operator()(const T&) const;
};

We can then provide an specialization for polynomial_function:

template<int N>
struct gradient<polynomial_function<N> > {
    typedef polynomial_function<N - 1> result_type;
    polynomial_function<N - 1> gradient(const polynomial_function<N>&) 
const;
};

Now, comes the tricky part.  The obvious method of storing both the function
and its derivative doesn't quite work, because that doesn't allow 
multiple derivatives.

Instead, we can store the function and another function that computes 
the derivative.

template<class F>
struct default_gradient_impl {
    typedef typename gradient<F>::result_type result_type;
    template<class Signature>
    result_type operator()(const boost::function<Signature>& f) {
        return(gradient<F>()(*f.target<F>()));
    }
};

template<class Signature>
class math_function {
public:
    template<class F>
    math_function(const F& f) : f_(f), gradient(default_gradient<F>()) {}
private:
    boost::function<Signature> f_;
    
boost::function<math_function<Signature>(boost::function<Signature>&)> 
gradient_;
};

template<class Signature>
struct gradient<math_function<Signature> > {
    typedef math_function<Signature> result_type;
    result_type operator()(const math_function<Signature>& f) { 
return(f.gradient()); }
};

One more subtlety: as it stands, the code for derivatives of all orders 
has to be
compiled.  This is going to be a problem if the derivative of F is 
computed by
finite_difference<F>, because we'll end up instantiating 
finite_difference<F>,
finite_difference<finite_difference<F> >, 
finite_difference<finite_difference<finite_difference<F> > >,
and so on ad infintum.  You'll have to provide some way to terminate the
recursion.

Finally, this implementation allows any number of derivatives to be 
specified explicitly
by adding appropriate constructors.

    template<class F, class D1>
    math_function(const F& f, const D1& d1)
        : f_(f), gradient(boost::lambda::constant(d1)) {}
    template<class F, class D1, class D2>
    math_function(const F& f, const D1& d1, const D2& d2)
        : f_(f), gradient(boost::lambda::constant(math_function(d1, d2))) {}

In Christ,
Steven Watanabe