Accept std::function with arbitrary inputs as input w/o Templates

For Learning Purposes:

I am creating a small numerical methods library and I am trying to implement the gradient currently I have done 2D gradient and 3D gradient . But I want to generalize this to higher dimensions.

Currently I have :

matrix<double> analysis::gradient_2D(std::function<double(double, double)> fn, double x, double y)
{
    matrix<double> R(2, 1);

    std::function<double(double)> fnX = [fn, y](double xVar){ return fn(xVar, y); };
    std::function<double(double)> fnY = [fn, x](double yVar){ return fn(x, yVar); };

    R(1, 1) = differentiateBest(fnX, x);
    R(1, 2) = differentiateBest(fnY, y);

    return R;
}

matrix<double> analysis::gradient_3D(std::function<double(double, double, double)> fn, double x, double y, double z)
{
    matrix<double> R(3, 1);

    std::function<double(double)> fnX = [fn, y, z](double xVar){ return fn(xVar, y,z); };
    std::function<double(double)> fnY = [fn, x, z](double yVar){ return fn(x ,yVar, z); };
    std::function<double(double)> fnZ = [fn, x, y](double zVar){ return fn(x, y, zVar); };

    R(1, 1) = differentiateBest(fnX, x);
    R(1, 2) = differentiateBest(fnY, y);
    R(1, 3) = differentiateBest(fnZ, z);

    return R;
}

// Where 

double analysis::differentiateBest(std::function<double(double)> fn, double x)
{
    return derivative_1_Richardson_6O(fn,  x);
}
// For brevity , derivative_1_Richardson_6O also has the same input as differentiateBest

I know it is verbose , but I like it

Question What I would like to do is to create a

// What do I do at the ...  ? 
matrix<double> analysis::gradient_ND(std::function<double(...)> fn, matrix<double>)

So that I can pass a std::function with arbitrary input say N and I will pass a vector which has N values.

How will I go about doing this ? If the answer is too long , links will be appreciated too . Thank you.

PS: I saw a method using Templates , but if I use templates in the implementation , the I will have to change the .cpp file to something else right ? I would like to avoid. If using templates is the only way , then I will have to compromise. Thanks.

Solution

template<class T>
struct array_view {
  T* b = nullptr; T* e = nullptr;
  T* begin() const { return b; }
  T* end() const { return e; }

  size_t size() const { return end()-begin(); }
  bool empty() const { return begin()==end(); }
  T& front()const{return *begin(); }
  T& back()const{return *std::prev(end()); }

  T& operator[](size_t i)const{return begin()[i]; }

  array_view( T* s, T* f ):b(s),e(f) {};
  array_view() = default;

  array_view( T* s, size_t l ):array_view(s, s+l) {}

  using non_const_T = std::remove_const_t<T>;
  array_view( std::initializer_list<non_const_T> il ):
    array_view(il.begin(), il.end()) {}
  template<size_t N>
  array_view( T(&arr)[N] ):array_view(arr, N){}
  template<size_t N>
  array_view( std::array<T,N>&arr ):array_view(arr.data(), N){}
  template<size_t N>
  array_view( std::array<non_const_T,N> const&arr ):
    array_view(arr.data(), N){}
  template<class A>
  array_view( std::vector<T,A>& v):
    array_view(v.data(), v.size()){}
  template<class A>
  array_view( std::vector<non_const_T,A> const& v):
    array_view(v.data(), v.size()){}
};

an array_view is a non-owning view into a contiguous array of T. It has converting constructors from a myriad of contiguous containers (if matrix is contiguous, a converter to array_view should be written).

Then:

matrix<double> analysis::gradient_ND(std::function<double(array_view<const double>)>, array_view<const double> pt)

is reasonable.

Using array_view causes no problem with .h and .cpp code splitting. The only templates involved are either fixed (the array_view itself), or are resolved when constructing the array_view.

matrix<double> R(pt.size(), 1);

auto make_tmp = [pt]{
  std::vector<double> tmp;
  tmp.reserve(pt.size());
  for (double x:pt)
    tmp.push_back(x);
  return tmp;
};

std::vector<std::function<double(double)>> partials;
partials.reserve(pt.size());
for (size_t i = 0; i < pt.size(); ++i) {
  partials.push_back(
    [&,i](double x){ auto tmp=make_tmp(); tmp[i]=x; return fn(tmp); };
  );
}
for (size_t i = 0; i < pt.size(); ++i) {
  R(1, i) = differentiateBest(partials[i], pt[i]);
}

return R;

note that creating array of partials is not actually needed. You could just directly differentiateBest.

There is an inefficiency where partials reallocate each call. If you are ok with making reentrancy not work (which will often be ok), creating a tmp and capturing it by-value, and modifying and restoring it after each call to fn, would boost performance.

  [&,i,tmp=make_tmp()](double x){ std::swap(tmp[i],x); double r=fn(tmp); std::swap(tmp[i],x); return r; };

is C++14 version. C++11 version would create a tmp variable and capture it by-value.