Hi,

The possibility of using functors with many of the NR3 routines is one of the features I love of this new edition. However, I ran into trouble in a particular case, when using the Stiel class to obtain abscissas and weights for Gaussian quadrature. The compiler cannot find an appropriate constructor that takes a functor instead of a function. Is there a way of doing this without changing the Stiel class?

Below you'll find the full error message followed by the full test code. If I replace the functor with a normal function it compiles and works fine.

Thanks in advance,
Jose Miguel Pasini

Error message:

stiel_test.cpp: In function ‘void gauss_hermite(Int, Doub, Int)’:
stiel_test.cpp:32: error: no matching function for call to ‘Stiel::Stiel(const Int&, Doub&, Doub&, ShiftedGaussian&, Doub (&)(Doub), Doub (&)(Doub))’
../NR_C301/code/stiel.h:69: note: candidates are: Stiel::Stiel(Int, Doub, Doub, Doub (*)(Doub), Doub (*)(Doub), Doub (*)(Doub))
../NR_C301/code/stiel.h:61: note: Stiel::Stiel(Int, Doub, Doub, Doub, Doub (*)(Doub, Doub))
../NR_C301/code/stiel.h:1: note: Stiel::Stiel(const Stiel&)
make: *** [stiel_test] Error 1


Full test code:

#include "nr3.h"
#include "sort.h"
#include "quadrature.h"
#include "derule.h"
#include "eigen_sym.h"
#include "gauss_wgts2.h"
#include "stiel.h"

Doub map_semiinf_x(const Doub t)
{
return exp(t-exp(-t));
}

Doub map_semiinf_dxdt(const Doub t)
{
return exp(t-exp(-t))*(1.0 + exp(-t));
}

struct ShiftedGaussian {
Doub x0; // position of peak
ShiftedGaussian(const Doub xx0) : x0(xx0) { }
Doub operator()(const Doub x)
{
return exp(-(x-x0)*(x-x0));
}
};

void gauss_hermite(const Int npoints, const Doub x0, const Int prec)
{
Doub a = -6.0, b = 6.0;
ShiftedGaussian sg(x0);
Stiel gauss(npoints,a,b,sg,map_semiinf_x,map_semiinf_dxd t);
VecDoub x(npoints), w(npoints);
gauss.get_weights(x,w);
}

int main(int argc, char *argv[])
{
gauss_hermite(4,2.0,5);
return 0;
}

The compiler cannot find an appropriate constructor that takes a functor instead of a function. Is there a way of doing this without changing the Stiel class?Haven't you answered your own question? (Actually the compiler answered it for us.) The Stiel class is not templated, so there is no way to use a functor-type argument.

If you don't want to change the class (a Good Decision, I'm thinking), then you must give its constructor an argument that is a pointer to a function. For your object, the function's return type must be a Doub and it must have one argument that is a Doub.

Instead of the functor class that you showed us, maybe you could try something like

//
// Keeps value of x0 for ShiftedGaussian
//
Doub xMinusX0(Doub x, bool init=false)
{
static Doub x0 = 0; // not shifted unless re-initialized with init == true
if (init) {
x0 = x;
}
return (x - x0);
}

//
//Calls xMinusX0 to get shifted argument
//
Doub ShiftedGaussian(Doub x)
{
Doub xmx0 = xMinusX0(x);
return exp(-(xmx0*xmx0));
}
.
.
.
void gauss_hermite(const Int npoints, const Doub x0, const Int prec)
{
Doub a = -6.0, b = 6.0;

xMinusX0(x0, true); // Initialise x0 for ShiftedGaussian
Stiel gauss(npoints, a, b, ShiftedGaussian, map_semiinf_x, map_semiinf_dxdt);
.
.
.
}


Regards,

Dave

Instead of the functor class that you showed us, maybe you could try something like
...
Regards,

Dave

Dave,

Thank you very much for the suggestion. That's a good working solution: I had considered having a global variable x0 (a solution I disliked very much), but I hadn't considered using the function xminusx0 with a static variable to achieve the same goal. However, there's one thing that I don't like: in effect, xminusx0 now behaves like an interface to a global variable x0, which is what I wanted to avoid in the first place: encapsulation has now been broken and we cannot be using xminusx0 with different values of x0 at the same time.

I know it sounds like I'm being too picky, but I'd like to learn how to do this correctly.

Thanks again,
Jose Miguel

in effect, xminusx0 now behaves like an interface to a global variable x0
From a debugging point of view, I find use of a static variable in a public function less onerous than use of global variables.

There is only one way to change x0 and that is to call the function, whereas with a global variable, a typographical error inside a block of code can inadvertently change the value. I think that it's easier to trap a spurious function call than put a watch point on some memory location in a debugger---but that depends on your development habits.

However...

I agree that a static variable in a public function is not ideal from a picky (thread-safe or other hoity-toity point view). In my experience, thread-safety is generally not an issue with my programs in which I would be using this class, but that might not always be the case.

So...

If you want to make the Stiel class constructor accept something like a functor that can keep the value of an internal (private) variable like x0, change it to a templated class or create a new templated class with the whatever functionality you need. Then the argument that you need can be something other than Doub(*)(Doub).

If it's worth it to you, then change the class. If not, then use a global variable or a static variable in a function.

Regards,

Dave

...

I think that it's easier to trap a spurious function call than put a watch point on some memory location in a debugger---but that depends on your development habits.

...

In my experience, thread-safety is generally not an issue with my programs in which I would be using this class, but that might not always be the case.

...

Regards,

Dave

Dave,

Thanks again for the advice. I think you're right in that it's less error-prone and that in any case thread-safety is not an issue for this particular class.

I'm new to NR3, so I feared that templating Stiel would require further templating up the chain, but I checked now and it appears not to be necessary: just templating Stiel should do the trick.

Thanks again,
Jose Miguel