Powell minimization

[dgholstein]

I've used NR3 interpolation classes to construct the following error function:

//    function calculates error function as a function of sigma and dielectric
//    result is square of differences
Doub error_sq_func(VecDoub x)
{
int i, j;
double sigma, dielectric, result, sum = 0.0;

sigma = x[0]; dielectric = x[1];
for (i=0; i<func_parameters.n_splines; i++) {
	if (eval_spline((unsigned long *) func_parameters.spline+i, sigma, dielectric, &result, func_parameters.error) < 0)
		return NaN;
	sum += sqr(func_parameters.vals[i] - result);
}

Of course, most of the code isn't shown here, but it has been tested and works well.

I just can't figure out how to call the Powell minimizer, all my attempts at building a constructor come up empty -- I've been trying to construct the second instance mentioned on page 513 of the book (with simple C++ function instead of functor).

Does anyone have an example of Powell minimization? Is there another place for NR3 examples?

...Dan

[davekw7x]

Originally Posted by dgholstein

Doub error_sq_func(VecDoub x)
{
.
.
.
}

...attempts at building a constructor come up empty

With that function, the constructor could be like this:

Powell <Doub(VecDoub)> powell(error_sq_func); // Make sure there is no '&' in the template function argument!

Although, I think many (most?) C++ programmers would prefer the form shown on page 513 (Pass a const reference to an object rather than passing a copy of the object):

Doub error_sq_func(VecDoub_I & x)
{
.
.
.
}

In which case the constructor could be like this:

Powell <Doub(VecDoub_I &)> powell(error_sq_func);

The function implementation code would be exactly the same, whichever style you use.

Originally Posted by dgholstein

Does anyone have an example of Powell minimization?

A program based on the one in the Numerical Recipes 2 C++ Examples:

#include "../code/nr3.h"
#include "../code/bessel.h"
#include "../code/mins.h"
#include "../code/mins_ndim.h"

using namespace std;

// Driver for Powell minimization
//
// Find minimum value of 0.5 - J0[(x-1)^2 +(y-2)^2 +(z-3)^2]
//

Doub func(VecDoub_I & x)
{
    Bessjy bessjy;
    return 0.5 - bessjy.j0(SQR(x[0]-1.0) + SQR(x[1]-2.0) + SQR(x[2]-3.0));
}

int main()
{

    int NDIM = 3;
    Doub initx[] = { // Initial guess: has NDIM elements
        1.5, 1.5, 2.5 
    };
    VecDoub p(NDIM, initx); // Initialize the vector with the guess
    VecDoub pp;

    // Feed the function to the constructor
    Powell <Doub(VecDoub_I &)> powell(func);

    // Perform the minimization
    pp = powell.minimize(p);

    cout << "Number of iterations =  " << powell.iter << endl << endl;

    cout << scientific;

    cout << "Minimum function value = ";
    cout << powell.fret
         << ", found at: " << endl;
    for (int i = 0; i < NDIM; i++) {
        cout << setw(15) << pp[i];
    }
    cout << endl << endl;


    cout << "Actual minimum value = "
         << -0.5
         << " at:" << endl;
    cout << setw(15) << 1.0 
         << setw(15) << 2.0 
         << setw(15) << 3.0 
         << endl;
    return 0;
}

Output

Number of iterations =  1

Minimum function value = -5.000000e-01, found at: 
   1.000000e+00   2.000000e+00   2.999954e+00

Actual minimum value = -5.000000e-01 at:
   1.000000e+00   2.000000e+00   3.000000e+00

Regards,

Dave

[dgholstein]

davekw7x;

I modeled my code based on your recommendations and it compiles and appears to work (I haven't yet rigorously tested it and still have to obtain experimental data for verification). Anyway, a million thanks!

C++ is difficult enough, the added layers of templates, typedefs, and polymorphism of the library make this old FORTRAN programmer (30 years, 20 years for C) strain to make sense of it.

Again, thanks.

...Dan