AnimaL

Tutorial

Documentation

---

Documentation

• Overview
• Modules
• Namespace List
• Class Hierarchy
• Alphabetical List
• Compound List
• File List
• Namespace Members
• Compound Members
• File Members
• Examples

---

powell.h

Go to the documentation of this file.

#ifndef NR_POWELL_H
#define NR_POWELL_H

#include <math.h>
#include "utils.h"
#include "linmin.h"

namespace nr {


template < class TraitsT = Old_Traits<> >
struct powell {

    typedef TraitsT Traits;
    typedef typename Traits::Real Real;
    typedef typename Traits::Vector Vector;
    typedef typename Traits::Matrix Matrix;

    
    template <class Func>
    powell (
        Vector& p, 
        Matrix& xi, 
        unsigned int n, 
        Real ftol, 
        unsigned int *iter, 
        Real *fret, 
        Func (*func), 
        unsigned int ITMAX = 200 
    )
    {
        unsigned int i,ibig,j;
        Real del,fp,fptt,t;
        Vector pt, ptt, xit;

        Traits::alloc_vector( pt, n );
        Traits::alloc_vector( ptt, n );
        Traits::alloc_vector( xit, n );

        *fret=(*func)(p);
        for (j=0;j<n;j++) pt[j]=p[j];
        for (*iter=1;;++(*iter)) {
            fp=(*fret);
            ibig=0;
            del=0.0;
            for (i=0;i<n;i++) {
                for (j=0;j<n;j++) xit[j]=xi[j][i];
                fptt=(*fret);
                linmin<Traits>(p,xit,n,fret,func,ftol);
                if (fabs(fptt-(*fret)) > del) {
                    del=fabs(fptt-(*fret));
                    ibig=i;
                }
            }
            if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
                Traits::free_vector( pt );
                Traits::free_vector( ptt );
                Traits::free_vector( xit );
                return;
            }
            if (*iter == ITMAX) {
                throw TooManyIterations( *iter, *fret );
            }
            for (j=0;j<n;j++) {
                ptt[j]=2.0*p[j]-pt[j];
                xit[j]=p[j]-pt[j];
                pt[j]=p[j];
            }
            fptt=(*func)(ptt);
            if (fptt < fp) {
                t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt);
                if (t < 0.0) {
                    linmin<Traits>(p,xit,n,fret,func);
                    for (j=0;j<n;j++) {
                        xi[j][ibig]=xi[j][n-1];
                        xi[j][n-1]=xit[j];
                    }
                }
            }
        }
    }

    template <class Func>
    powell (
        Vector& p, 
        unsigned int n, 
        Real ftol, 
        unsigned int *iter, 
        Real *fret, 
        Func (*func),
        unsigned int ITMAX = 200 
    )
    {
        Matrix unit_xi;
        Traits::alloc_matrix( unit_xi, n, n );
        for( int i=0; i<n; ++i ){
            for( int j=0; j<n; ++j ){
                unit_xi[i][j] = i==j ? 1 : 0;
            }
        }

        powell( p, unit_xi, n, ftol, iter, fret, func, ITMAX );

        Traits::free_matrix(unit_xi,n,n);

    }
    

    
    struct TooManyIterations 
    {
        unsigned int nbIterations;   
        Real functionValue;          

        TooManyIterations ( unsigned int n, Real f )
            : nbIterations( n )
            , functionValue( f )
        {}
        


        void print ( std::ostream& str ) const
        {
            str << "powell gives up after " << nbIterations
                << " iterations, function value = " << functionValue
                << endl;
        }
        

    };

};

} // end namespace nr
#endif

---