Wrapper class for a function with derivative

epiliam 0
Viewed 29

I want to write a wrapper class for a function which has a value and a derivative. What I have gone with so far is this:

template<typename input, typename output, typename derivative_output>
class Function
{
    public:
        Function() {}
        ~Function() {}
        Function(std::function<output(input)> &val, std::function<derivative_output(input)> &deriv) : __value(val), __derivative(deriv) {}

        output value(input x)
        {
            return __value(x);
        }
        derivative_output derivative(input x)
        {
            return __derivative(x);
        }

    private:
        std::function<output(input)> __value;
        std::function<derivative_output(input)> __derivative;
};

template<typename input, typename output, typename derivative_output>
Function<input, output, derivative_output> operator +(Function<input, output, derivative_output> func1, Function<input, output, derivative_output> func2) // addition of two functions
{
    std::function<output(input)> value_sum = [&func1, &func2](input x) {return func1.value(x) + func2.value(x);};
    std::function<derivative_output(input)> derivative_sum = [&func1, &func2](input x) {return func1.derivative(x) + func2.derivative(x);};
    return Function<input, output, derivative_output>(value_sum, derivative_sum);
}

template<typename input, typename output, typename derivative_output>
Function<input, output, derivative_output> operator *(double scalar, Function<input, output, derivative_output> func) // scalar multiplication
{
    std::function<output(input)> scaled_value = [&scalar, &func](input x) {return scalar * func.value(x);};
    std::function<derivative_output(input)> scaled_derivative = [&scalar, &func](input x) {return scalar * func.derivative(x);};
    return Function<input, output, derivative_output>(scaled_value, scaled_derivative);
}

using ScalarFunction1D = Function<double, double, double>; // maps scalars to scalars
using ScalarFunction2D = Function<Eigen::Vector2D, double, Eigen::Vector2D>; // maps vectors (in 2D space) to scalars
using VectorFunction1D = Function<double, Eigen::Vector2D, Eigen::Vector2D>; // maps scalars to vectors (in 2D space)
using VectorFunction2D = Function<Eigen::Vector2D, Eigen::Vector2D, Eigen::Matrix2D>; // maps vectors (in 2D space) to vectors (in 2D space)

ScalarFunction1D ScalarMonomial1D(int degree)
{
    std::function<double(double)> val = [&degree](double x) {return std::pow(x, degree);};
    std::function<double(double)> deriv = [&degree](double x) {return (degree == 0 ? 0.0 : double(degree) * std::pow(x, degree - 1));};
    return ScalarFunction1D(val, deriv);
}

I have a quick test case as follows

int main()
{
     ScalarFunction1D poly_1 = ScalarMonomial1D(2) + 5.0 * ScalarMonomial1D(1) + 6.0 * ScalarMonomial1D(0); // x^2 + 5x + 6
     ScalarFunction1D poly_2 = ScalarMonomial1D(2); // x^2

     std::cout << poly_1.value(0.0) << "\n"; //should be 6
     std::cout << poly_1.derivative(0.0) << "\n"; //should be 5
     std::cout << "\n";
     std::cout << poly_2.value(2.0) << "\n"; //should be 4
     std::cout << poly_2.derivative(2.0) << "\n"; //should be 4
}

It compiles and runs, but the output I get is

0
0

Inf
Inf

I have two questions: Firstly, why is this failing? Secondly, is the way I have designed my wrapper class good, or are there better alternatives? I thought about having function pointers or lambdas for a while, but ended up going with std::function. I want to be able to add and scalar multiply functions and will also add thing like scalar product and dot product eventually. I also want for the class to be completely abstract (in the mathematical not programming sense) so I can have any sort of instance of it e.g. polynomials, exponentials, whatever. Both memory and cpu efficiency may be an issue at some stage, although probable cpu over memory.

Regarding the first question, my guess is it is to do with how I am capturing variables in the std::function. This is an issue in the operators * and + and the function ScalarMonomial1D(int). My guess is that when I have std::function<double(double)> val = [&degree](double x) {return std::pow(x, degree);}; degree is a local variable so when it is eventually called it no longer exists. Same goes for the + and * operators. Is this the case and if so how can I circumvent it?

c++
0 Answers
Related