Numerical Recipes Forum  

Go Back   Numerical Recipes Forum > Numerical Recipes Third Edition Forum > Methods: All Chapters in NR3

Reply
 
Thread Tools Rate Thread Display Modes
  #1  
Old 06-23-2008, 06:31 PM
dgholstein dgholstein is offline
Registered User
 
Join Date: Jun 2008
Location: Sparta, NJ
Posts: 2
Powell minimization

I've used NR3 interpolation classes to construct the following error function:
Code:
//    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
Reply With Quote
  #2  
Old 06-23-2008, 10:23 PM
davekw7x davekw7x is offline
Registered User
 
Join Date: Jan 2008
Posts: 453
Quote:
Originally Posted by dgholstein
Code:
Doub error_sq_func(VecDoub x)
{
.
.
.
}
...attempts at building a constructor come up empty
With that function, the constructor could be like this:
Code:
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):

Code:
Doub error_sq_func(VecDoub_I & x)
{
.
.
.
}
In which case the constructor could be like this:
Code:
Powell <Doub(VecDoub_I &)> powell(error_sq_func);
The function implementation code would be exactly the same, whichever style you use.


Quote:
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:

Code:
#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
Code:
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

Last edited by davekw7x; 06-24-2008 at 11:07 AM.
Reply With Quote
  #3  
Old 06-25-2008, 09:40 AM
dgholstein dgholstein is offline
Registered User
 
Join Date: Jun 2008
Location: Sparta, NJ
Posts: 2
Thumbs up Thanks

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
Reply With Quote
Reply

Thread Tools
Display Modes Rate This Thread
Rate This Thread:

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off

Forum Jump


All times are GMT -5. The time now is 09:40 AM.


Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2014, Jelsoft Enterprises Ltd.