Hi

Please help me. For a week or so I have been trying
to figure out how to use an optimizer (NR::amoeba)
to max a function say opt_func in my class called
model.

Can anybody tell me how it should be decleared, and
how the amoeba call should look like?

Please help me.

Best regards Asger Lunde

#########################

Asger Lunde
Associate Professor, Ph.D.

Business Address:
Department of Information Science
The Aarhus School of Business
Fuglesangs Allé 4
8210 Aarhus V
Denmark

Tlf work: (+45) 8948-6360
Fax: (+45) 8615-3792

E-mail: alunde@asb.dk
Web page: www.asb.dk/~alunde

#########################

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::objFunc(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::objGrad(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");
}