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.

How can class member functions be used with the NR minimizing functions, such as NR::dfpmin?

**************************************************
****************************
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");
}

Hi, as it's been a couple of years since you posted this message, you've probably figured it out... but I ran into the same issue when I was implementing a couple of the routines here, so here's my solution. Much of this was learned from reading Stroustrup's book and poking around on my own.

Pointers to non-static member functions need to have a 'this' pointer associated with them, so that they can access the member objects (and make the right function call!) for the instantiation you are using. That's why static member functions work without a problem: they (by definition) do not have a 'this' pointer, and as a result cannot modify instantiation-specific variables. But naturally when we define objects that do optimization for us, we want to keep context within the object, and have that context available to our 'pointed-to' member functions.

There is specific C++ syntax for a pointer-to-member-function, and for dereferencing such a pointer while associating an appropriate 'this' with it. You probably tried what you know from C, as did I:

#include <vector>
class Optim {
public:
typedef std::vector<double> v_type;

double eval1(const v_type& p);
double eval2(const v_type& p);

void optimize(const v_type& parms, double func(const v_type& p)) {

double result = func(parms);

double result = (*func)(parms); // alternate syntax
...
}
};
...
Optim o;
Optim::v_type q;
...
o.optimize(q, &Optim::eval1);

But this won't compile in C++. The solution is to use the special syntax, which I show below. I recommend the use of a typedef for readability:

#include <vector>
class Optim {
public:
typedef std::vector<double> v_type;

double eval1(const v_type& p);
double eval2(const v_type& p);

// ptr_func is the type of a pointer to one of our two
// eval member functions within class Optim.

typedef double (Optim::* ptr_func)(const v_type&);

// Use the typedef in the declaration of the function

void optimize(const v_type& parms, ptr_func func) {

// In the function body, you want 'func' to be
// evaluated within the same Optim instantiation
// as the call to 'optimize', that is, you want
// 'optimize' to share its 'this' pointer with 'func'.
// Use ->* or .* to supply the 'this' pointer,
// depending on whether you need to dereference.
// For our example, we do (this->*func)().

double result = (this->*func)(parms);

...
}
};
...
Optim o;
Optim::v_type q;
Optim::ptr_func the_eval_func_of_choice = &Optim::eval1;
o.optimize(q, the_eval_func_of_choice);
o.optimize(q, &Optim::eval1); // alternate
...


Good luck!