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