I developed a rudimentary class Optim (see below) for use with the NR minimizing functions but cannot seem to get the C++ compiler to recognize the function pointers (usually called "func" and "dfunc") used as arguments in the NR minimization routines. For example, member functions of class Optim specify the objective function, the gradient, and the quasi-Newton algorithm illustrated in the Numerical Recipes Example file "xdfpmin.cpp". The compiler cannot match the func and dfunc arguments in Optim::dfpmin with the parameters objFunc and objGrad that are passed in the calling function Optim::fit.

I'm not sure why this doesn't work. If objFunc and objGrad are declared as static members, then the source code compiles and runs properly! Of course, this solution is not really preferred.

This seems like a fairly general problem:

How can class member functions be used with NR functions whose arguments are function pointers?

**************************************************

****************************
Bob Dorazio
United States Geological Survey
7920 NW 71 Street
Gainesville, Florida 32653
Phone: (352) 378-8181 x373
Fax: (352) 378-4956
Email: bdorazio@usgs.gov
Internet: www.fcsc.usgs.gov
**************************************************

****************************
#include <cmath>
#include <limits>
#include "nr.h"
using NR::lnsrch;
using NR::nrerror;
using namespace std;

#include "mtl/mtl.h"
#include "mtl/dense1D.h"
#include "mtl/utils.h"
using namespace mtl;

class Optim {

public:
Optim(const dense1D<double>& p, const DP GTOL);

DP objFunc(Vec_I_DP &x);
void objGrad(Vec_I_DP &x, Vec_O_DP &df);
void dfpmin(Vec_IO_DP &p, DP gtol, int &iter, DP &fret,
DP func(Vec_I_DP &), void dfunc(Vec_I_DP &, Vec_O_DP &));

void fit();
void PrintParamValue();
void PrintFuncValue();

private:
int niter;
DP funcValue;
DP tolerance;
dense1D<double> pGuess;
};


Optim::Optim(const dense1D<double>& p, const DP GTOL) {
pGuess.resize(p.size());
copy(p,pGuess);
tolerance = GTOL;
niter = 0;
}

void Optim::PrintParamValue() {
print_vector(pGuess);
}

void Optim::PrintFuncValue() {
cout << objFunc(pGuess) << endl;
}


void Optim::fit() {
dfpmin(pGuess,tolerance,niter,funcValue,objFunc,ob

jGrad);
}


DP Optim:bjFunc(Vec_I_DP &x)
{
DP x1p2sqr=SQR(2.0+x[0]);

return 10.0*SQR(SQR(x[1])*(3.0-x[0])-SQR(x[0])*(3.0+x[0]))+
x1p2sqr/(1.0+x1p2sqr);
}

void Optim:bjGrad(Vec_I_DP &x, Vec_O_DP &df)
{
DP x1sqr=SQR(x[0]),x2sqr=SQR(x[1]),x1p2=x[0]+2.0;
DP x1p2sqr=SQR(x1p2);

df[0]=20.0*(x2sqr*(3.0-x[0])-x1sqr*(3.0+x[0]))*
(-x2sqr-6.0*x[0]-3.0*x1sqr)+2.0*x1p2/(1.0+x1p2sqr)
-2.0*x1p2*x1p2sqr/SQR((1.0+x1p2sqr));
df[1]=40.0*(x2sqr*(3.0-x[0])-x1sqr*(3.0+x[0]))*x[1]*(3.0-x[0]);
}


void Optim::dfpmin(Vec_IO_DP &p, DP gtol, int &iter, DP &fret,
DP func(Vec_I_DP &), void dfunc(Vec_I_DP &, Vec_O_DP &))
{
const int ITMAX=200;
const DP EPS=numeric_limits<DP>::epsilon();
const DP TOLX=4*EPS,STPMX=100.0;
bool check;
int i,its,j;
DP den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp ,test;

int n=p.size();
Vec_DP dg(n),g(n),hdg(n),pnew(n),xi(n);
Mat_DP hessin(n,n);
fp=func(p);
dfunc(p,g);
for (i=0;i<n;i++) {
for (j=0;j<n;j++) hessin[i][j]=0.0;
hessin[i][i]=1.0;
xi[i] = -g[i];
sum += p[i]*p[i];
}
stpmax=STPMX*MAX(sqrt(sum),DP(n));
for (its=0;its<ITMAX;its++) {
iter=its;
lnsrch(p,fp,g,xi,pnew,fret,stpmax,check,func);
fp=fret;
for (i=0;i<n;i++) {
xi[i]=pnew[i]-p[i];
p[i]=pnew[i];
}
test=0.0;
for (i=0;i<n;i++) {
temp=fabs(xi[i])/MAX(fabs(p[i]),1.0);
if (temp > test) test=temp;
}
if (test < TOLX)
return;
for (i=0;i<n;i++) dg[i]=g[i];
dfunc(p,g);
test=0.0;
den=MAX(fret,1.0);
for (i=0;i<n;i++) {
temp=fabs(g[i])*MAX(fabs(p[i]),1.0)/den;
if (temp > test) test=temp;
}
if (test < gtol)
return;
for (i=0;i<n;i++) dg[i]=g[i]-dg[i];
for (i=0;i<n;i++) {
hdg[i]=0.0;
for (j=0;j<n;j++) hdg[i] += hessin[i][j]*dg[j];
}
fac=fae=sumdg=sumxi=0.0;
for (i=0;i<n;i++) {
fac += dg[i]*xi[i];
fae += dg[i]*hdg[i];
sumdg += SQR(dg[i]);
sumxi += SQR(xi[i]);
}
if (fac > sqrt(EPS*sumdg*sumxi)) {
fac=1.0/fac;
fad=1.0/fae;
for (i=0;i<n;i++) dg[i]=fac*xi[i]-fad*hdg[i];
for (i=0;i<n;i++) {
for (j=i;j<n;j++) {
hessin[i][j] += fac*xi[i]*xi[j]
-fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
hessin[j][i]=hessin[i][j];
}
}
}
for (i=0;i<n;i++) {
xi[i]=0.0;
for (j=0;j<n;j++) xi[i] -= hessin[i][j]*g[j];
}
}
nrerror("too many iterations in dfpmin");
}

According to Barton and Nackman (1994), function pointers cannot point to a non-static member function because a non-static member can only be invoked with an instance of its class. So instead of using function pointers, one solution is to remove the pointers from the argument list then simply declare these as pure virtual functions of an abstract class, which can be overriden by a derived class's member functions. This solution was suggested in an earlier posting "passing member functions ..." by mpearl.

Bob:

You cannot use a non-static member function because such functions implicitly require the "this" pointer to identify the object.

The solution given by mpearl is one way to go, but there are a couple of other solutions which do less violence to the organization of the Numerical Recipes library (you can leave the code as standalone functions):

1. You can use a static member function of your class. This works, but if your function has auxiliary parameters, then they must be static members of your class, and so not unique to your particular object.

2. You can use a function object. You can read about these in Stroustrup's book, or in any book on the Standard Template Library. This method requires you to modify the prototype (but not the code) of the Numerical Recipes routines, but it is (IMO) a much better solution, because it avoids the use of static data. Very briefly, a function object is an instance of a class which overloads the () operator. First make an abstract base class, say

namespace NR
{
class ScalFunc
{
virtual DP operator()(const DP x)=0;
virtual ~ScalFunc(){}
};
}

You can put this class at the top of your nr.h. Now in your Numerical Recipes routine, change the argument

DP f(const DP) // pointer to function

to

ScalFunc &f // reference to function object

Remember to this both in NR.h and in the code (and do likewise for any subsidiary functions if necessary). Both the old and new f have the same call semantics, so you do not have to change anything else.

To use the modified NR routine, derive from ScalFunc your own function class which implements the () operator, and pass an instance of this class to the NR routine. You can pass any additional parameters in the constructor of your function object. Or you can pass a pointer or reference to your original object (of class Optim, say) and use this pointer to call a non-static member function as you wanted to originally.

You will see that both the function object and the "mpearl" method make use of a virtual function to do the function evaluation. Take your choice.

David Wilkinson

Is there any criterion to set the stpmax using following expression:

stpmax=STPMX*MAX(sqrt(sum),DP(n));

in routine newt()?