-2

I am creating a class which takes as input arguments either 1 or 2 functions.

My goal is that if only one function func is given then the member function dfunc is calculated using num_dfunc (which is the numerical derivative of func and is hardcoded inside the class). If two functions are given func and analytical_dfunc then the derivative is calculated using the analytical version. What is the best way to achieve this?

This is a portion of my code

class MyClass
{
public:
    int dim = 2;
    vector<double> num_dfunc(vector<double> l0)
    {
        // Numerical gradient of the potential up to second order
        // #TODO This should be rewritten!
        vector<double> result(dim);
        double eps = 0.001;

        for (int i = 0; i < dim; i++)
        {
            vector<double> lp2 = l0;
            lp2[i] += 2 * eps;
            vector<double> lp1 = l0;
            lp1[i] += eps;
            vector<double> lm1 = l0;
            lm1[i] -= eps;
            vector<double> lm2 = l0;
            lm2[i] -= 2 * eps;
            result[i] = (-func(lp2) + 8 * func(lp1) - 8 * func(lm1) + func(lm2)) / (12 * eps);
        }
        return result;
    }
    double (*func)(vector<double>);   // Potencial pointer
    vector<double> (*dfunc)(vector<double>);  // Gradient pointer
    MyClass(double (*func)(vector<double>))
    {
        this->func = func;
        // THIS IS WRONG
        this->dfunc = num_dfunc;
    }
    MyClass(double (*func)(vector<double>),double (*analytical_dfunc)(vector<double>))
    {
        this->func = func;
        // THIS IS WRONG
        this->dfunc = analytical_dfunc;
    }

This is a, somewhat, pythonic way of what I want to do.

PS: This is not what I have so far, I tried many things but none of them worked.

Edit : Error, on dfunc return type. Typo on analytical_dfunc

João Viana
  • 31
  • 1
  • 5
  • 1
    `num_dfunc` and `analytical_func` have different return types. Did you mean for `analytical_func` to return a `vector`? – Nelfeal Sep 21 '22 at 15:33
  • Side note: `vector num_dfunc(vector l0)` had me wondering vector *ten*?!?! Dah hell? I probably won't be the only one, so you might want to make that identifier a bit less dependent on a good font. – user4581301 Sep 21 '22 at 15:41
  • @Nelfeal Correct, I will correct it. – João Viana Sep 21 '22 at 15:44

1 Answers1

0

I would make num_dfunc either a static member function, or probably even better a free function, that has nothing to do with MyClass. I would modify it, such that it also takes the zeroth-order function as well as the dimension as an input.

If your constructor is called with just one argument, you can create the second function from the static num_dfunc using a lambda.

[=](auto const& arg){
    return num_dfunc(func, dim, arg);
}

the notation [=] will capture all needed variables (func, dim) in the lambda body, that are not part of the argument list (auto const& arg).

Also I would make your class a class template accepting any type of function F, C++ has many different kinds of callable objects, of which raw function pointers are just one (there are also functors, lambdas, std::function, ...). If you make your class a template, it will work with all kinds of function types.

I haven't tested this, but basically your class would look like this:

#include <vector>


template <typename F>
class MyClass
{
private:

    static std::vector<double> num_dfunc(
        F const& f,
        int dim,
        std::vector<double> const& l0
    )
    {
        // Numerical gradient of the potential up to second order
        // #TODO This should be rewritten!
        std::vector<double> result(dim);
        double eps = 0.001;

        for (int i = 0; i < dim; i++)
        {
            std::vector<double> lp2 = l0;
            lp2[i] += 2 * eps;
            std::vector<double> lp1 = l0;
            lp1[i] += eps;
            std::vector<double> lm1 = l0;
            lm1[i] -= eps;
            std::vector<double> lm2 = l0;
            lm2[i] -= 2 * eps;
            result[i] = (-f(lp2) + 8 * f(lp1) - 8 * f(lm1) + f(lm2)) / (12 * eps);
        }
        return result;
    }
    
    int dim;
    F func;   // Potencial pointer
    F dfunc;  // Gradient pointer

public:
    MyClass(F const& fun)
     : dim(2)
     , func(fun)
     , dfunc(
            [=](auto const& arg){
                return num_dfunc(func, dim, arg);
            }
      )
    {
    }
    MyClass(F const& fun, F const& analytical_fun)
     : dim(2)
     , func(fun)
     , dfunc(analytical_fun)
    {
    }
};

You can play around on compiler explorer with this: https://godbolt.org/z/EMxqr7KE4

joergbrech
  • 2,056
  • 1
  • 5
  • 17