AnimaL

Tutorial

Documentation

---

Documentation

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

---

frprMinimizer.h

Go to the documentation of this file.

#ifndef animal_minimizer_hpp
#define animal_minimizer_hpp

#include <iostream>
#include <fstream>

namespace animal {

using std::ofstream;
using std::cout;
using std::endl;
using std::flush;

template <class ValueFunction,      // a function which computes the value of the function for given parameters
          class GradientFunction,   // a function which computes the gradient of the function for given parameters
          class Parameter,          // the parameter given to the function (possibly a container of parameters)
          class Real                // floating-point values
          >
class FrprMinimizer 
{
public:
    ValueFunction& valueFunction;       
    GradientFunction& gradientFunction; 
    
    // main loop control
    Real       mainTol;                 
    unsigned int mainMaxIter;             
    int restartFreq;                     
private:
    int nitermod;                     
    bool followGivenDirection;          
public: 
    // line minimization control
    Real linminTol;   
    int linminMaxIter;   
    Real linminEps;   
    Real previousStep; 
    
    // stop criteria
    bool nullGradient;     
    bool precisionReached; 
    bool maxStepReached;  
    Real cumulatedStep;  
    bool tooManyIterationsDbrent;  
    bool negativeStep;  

    static Real nullGradientThreshold(){ return 1.0e-9; }

    // result
    Real finalValue;      
    unsigned int nbIterations;       
    bool profiling;   
    Real dumpThreshold; 
    int nbValueEvaluations;  
    int nbGradientEvaluations;  
    
public:
    
    FrprMinimizer( 
        ValueFunction& valueFunction,
        GradientFunction& gradientFunction,
        Real mainTol=0.01, int mainMaxIter=1,
        Real linminTol=0.01, int linminMaxIter=10
        , bool profiling=false
        , int restartFreq=1000000 
        , Real linminEps=1.0e-10
        , Real dumpThreshold = 1.0e-9
        ):
        valueFunction(valueFunction),
        gradientFunction(gradientFunction),
        mainTol(mainTol), mainMaxIter(mainMaxIter)
        , restartFreq(restartFreq)
        , followGivenDirection(false)
        , linminTol(linminTol), linminMaxIter(linminMaxIter)
        , linminEps(linminEps)
        , previousStep(0)
        , tooManyIterationsDbrent(false)
        , profiling(profiling)
        , dumpThreshold(dumpThreshold)
    {}
    
    void profile( bool b ){ profiling=b; }
    
    
    void operator () ( Parameter& param, Real maxStep )
    {
        
        init(animal::size(param));

        while (
            nbIterations<mainMaxIter && 
            !nullGradient && 
            !precisionReached &&
            !maxStepReached &&
            !negativeStep ) 
        {
            // current value of the function
            Real fp=(valueFunction)(param); nbValueEvaluations++;
            
            // possibly use a given initial direction
            if( followGivenDirection )
            {
                animal::v_eq(xi,initialDirection);
                followGivenDirection = false;
                //if( profiling ){
                //  cerr <<"dot product between search directions: " << animal::v_dot(xi,initialDirection) << endl;
                //}
            }
            else {
                gradientFunction(param,xi); nbGradientEvaluations++;            
            }
            
            // new search direction
            if( nitermod == 0 )  // use gradient as search direction
            {
                //cerr << "niter = " << niter << ", restart frprmn " << endl;
                for (std::size_t j=0;j<size;j++) {
                    g[j] = -xi[j];
                    xi[j]=h[j]=g[j];
                }
            }
            else    // combine gradient with previous search direction
            {
                Real dgg=0.0, gg=0.0;
                //cerr <<"continue frprmn after next computation" << endl;
                for (std::size_t j=0;j<size;j++) {
                    gg += g[j]*g[j];
                    dgg += (xi[j]+g[j])*xi[j];
                }
                //cerr << "frprmn:: linear combination done " << endl;
                if ( gg <= nullGradientThreshold() )
                    nullGradient = true;
                else {
                    Real gam=dgg/gg;
                    for (std::size_t j=0;j<size;j++) {
                        g[j] = -xi[j];
                        xi[j]=h[j]=g[j]+gam*h[j];
                    }
                }
            }

            // compute a step in the search direction
            Real step = linemin(param,xi,&finalValue,maxStep,fp,linminTol);
            previousStep = step; // (non pertinent si plusieurs iterations)
            if (step<0) {
                step=0;
                negativeStep = true;
            }
            else if( (cumulatedStep + step) > maxStep ){ // limit reached
                step = maxStep - cumulatedStep;
                maxStepReached = true;
            }
            
            // apply the step
            for (unsigned int j=0;j<size;j++) {
                xi[j] *= step;
                param[j] += xi[j];
            }
            cumulatedStep += step;

            // test if desired precision is reached
            if (2.0*Rfabs(finalValue-fp) <= mainTol*(Rfabs(finalValue)+Rfabs(fp)+linminEps))
                precisionReached = true;
            //cerr <<"continue frprmn " << endl;
            
            nbIterations++;
            nitermod++;
            if( restartFreq!=0 ) 
                nitermod = nitermod % restartFreq;
        }
        
        if( profiling ) printLog(maxStep);
    }
    
    void initSearchDirection( const Parameter& point )
    {
        followGivenDirection = true;
        gradientFunction(point); nbValueEvaluations++;
        initialDirection.resize(animal::size(point));
        xi.resize(animal::size(point));
        gradientFunction(point,initialDirection); nbGradientEvaluations++;
    }
    
    void printFunction( Parameter& p )
    {
        init( animal::size(p) );
        xi.resize( animal::size(p) );
        cout << "function value: " << valueFunction(p) << flush << endl; 
        nbValueEvaluations++;
        gradientFunction(p,xi); nbGradientEvaluations++;
        cout << "gradient: \n" << v_output(xi) << flush << endl;
    }


    
private:

    void printLog( Real ) 
    {

        std::cerr<< nbValueEvaluations << " function evaluations " << endl;
        std::cerr<< nbGradientEvaluations << " gradient evaluations " << endl;
        if( maxStepReached ) std::cerr<< "max step reached " << endl;
        if( precisionReached ) std::cerr<< "precision reached " << endl;
        if( nullGradient ) std::cerr<< "nullGradient reached " << endl;
        if( nbIterations >= mainMaxIter ) std::cerr<< "maxiter reached " << endl;
        if( negativeStep ) std::cerr<< "negative step encountered " << endl;
        std::cerr << "minimization terminated in " << nbIterations << " iterations, endValue = " << finalValue << endl;
        if( cumulatedStep<0 ) cerr<<"*********** negative step ************** " << endl;
        std::cerr << "cumulated step = " << sqrt(cumulatedStep) << endl;
        valueFunction.display();
        //dumpPotential(-maxStep*1000, maxStep*1000, 1000, "potential.dmp" );
            //cerr << "potential dumped to file " << dumpfile << endl;
        std::cerr << endl;
    }
    
    void init(int s) 
    {
        size=s;
        g.resize(size);
        h.resize(size);
        if( !followGivenDirection ) xi.resize(size);
        pcom.resize(size);
        xicom.resize(size);
        xt.resize(size);
        df.resize(size);
        initialDirection.resize(size);

        maxStepReached = false;
        precisionReached = false;
        nullGradient = false;
        tooManyIterationsDbrent = false;
        negativeStep = false;
        nbIterations = 0;
        nitermod = 0;
        cumulatedStep = 0.0;
        
        nbValueEvaluations = nbGradientEvaluations = 0;
    }
    

    // auxiliary values
    unsigned int size;  
    Parameter g;  
    Parameter h;  
    Parameter xi;  
    Parameter pcom;  
    Parameter xicom;  
    Parameter xt;  
    Parameter df;  
    Parameter initialDirection;  

Real f1dim (Real x)
{
    if( profiling ) cerr << "compute f1dim(" << x<< ")" << endl;
    for (unsigned int j=0;j<size;j++) 
        xt[j]=pcom[j]+x*xicom[j];
     nbValueEvaluations++;
     return (valueFunction)(xt);
}

Real df1dim (Real x)
{
    unsigned int j;
    Real df1=0;

    for (j=0;j<size;j++) 
        xt[j]=pcom[j]+x*xicom[j];
    (gradientFunction)(xt,df); nbGradientEvaluations++;
    for (j=0;j<size;j++) 
        df1 += df[j]*xicom[j];
    return df1;
}

Real linemin (
    Parameter& p, 
    Parameter& xi, 
    Real *fret,
    Real maxStep,
    Real finit,
    Real TOL = 2.0e-4
)
{
    unsigned int j;
    Real xx,xmin,fx,fb,fa,bx,ax;
        
    for (j=0;j<size;j++) {   // pour df1dim
        pcom[j]=p[j];
        xicom[j]=xi[j];
    }
    ax=0.0;     fa=finit;
    xx=maxStep; fx=f1dim(xx);
    
    if(profiling){
        //cerr<< "linemin from " << v_output(p)  << endl;
        //cerr<< "linemin along " << v_output(xi)  << endl;
        //cerr << "linmin, f(ax) = "<< f1dim(ax) <<", f1dim(maxStep) = " << f1dim(maxStep) << endl;
    }

    //if( fx<fa ){      // euler step is OK
    //  *fret = fx;
    //  return maxStep;
    //}
    
    // bracket the minimum
    Real fcoherent = f1dim(previousStep); // try temporal coherency
    if( fcoherent<fa && fcoherent< fx ) 
    {                                       // minimum is bracketed 
        bx=xx; fb=fx;
        xx=previousStep; fx=fcoherent;
        if(profiling) cerr << "temporal coherency succeeded" << endl;
    }
    else 
    mnbrak(&ax,&xx,&bx,&fa,&fx,&fb);

    if( profiling ) {
        cout<<"bracketed: ("<<ax<<", "<<fa
        <<"), ("<<xx<<", "<<fx
        <<"), ("<<bx<<", "<<fb << ")" << endl;
    }

    // find the minimum
    *fret= dbrent(ax,xx,bx,TOL,&xmin); // attention, appliquer df1dim suppose que le dernier calcul de valeur dans mnbrak s'est effectué en xx. A verifier.
    //animal::v_eq(initialDirection,df);  // to reuse the gradient at the next step. This value will be used as initial search direction
    //cerr << "positions used for temporal cherency, method 2 " << endl << v_output(xt) << endl;
    if( profiling )
        cout<<"dbrent finds the minimum at (" << xmin <<","<<*fret<<")" <<endl;

    return xmin;
}
    

void mnbrak
(
    Real *ax, Real *bx, Real *cx,
    Real *fa, Real *fb, Real *fc
)
{
    const Real GOLD = 1.618034;
    const Real GLIMIT = 100.0;
    const Real TINY = 1.0e-20;
    const Real TWO = 2.0;
    #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
    #define SIGN(a,b) ((b) >= 0.0 ? Rfabs(a) : -Rfabs(a))


    int nbfunc=nbValueEvaluations;     // for profiling
    int nbgrad=nbGradientEvaluations;

    Real ulim,u,r,q,fu,dum;

    //*fa=f1dim(*ax);   calcules dans linemin
    //*fb=f1dim(*bx);
    if (*fb > *fa) {
        SHFT(dum,*ax,*bx,dum)
        SHFT(dum,*fb,*fa,dum)
    }
    *cx=(*bx)+GOLD*(*bx-*ax);
    *fc=f1dim(*cx);
    while (*fb > *fc) {
        r=(*bx-*ax)*(*fb-*fc);
        q=(*bx-*cx)*(*fb-*fa);
        u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/(TWO*SIGN(MAX(Rfabs(q-r),TINY),q-r));
        ulim=(*bx)+GLIMIT*(*cx-*bx);
        if ((*bx-u)*(u-*cx) > 0.0) {
            fu=f1dim(u);
            if (fu < *fc) {
                *ax=(*bx);
                *bx=u;
                *fa=(*fb);
                *fb=fu;
                if( profiling ) cerr <<"mnbrak: " << nbValueEvaluations-nbfunc << " value evaluations, " << nbGradientEvaluations-nbgrad << " gradient evaluations" << endl;
                return;
            } else if (fu > *fb) {
                *cx=u;
                *fc=fu;
                if( profiling ) cerr <<"mnbrak: " << nbValueEvaluations-nbfunc << " value evaluations, " << nbGradientEvaluations-nbgrad << " gradient evaluations" << endl;
                return;
            }
            u=(*cx)+GOLD*(*cx-*bx);
            fu=f1dim(u);
        } else if ((*cx-u)*(u-ulim) > 0.0) {
            fu=f1dim(u);
            if (fu < *fc) {
                SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
                SHFT(*fb,*fc,fu,f1dim(u))
            }
        } else if ((u-ulim)*(ulim-*cx) >= 0.0) {
            u=ulim;
            fu=f1dim(u);
        } else {
            u=(*cx)+GOLD*(*cx-*bx);
            fu=f1dim(u);
        }
        SHFT(*ax,*bx,*cx,u)
        SHFT(*fa,*fb,*fc,fu)
    }
    #undef SHFT 
    #undef SIGN 
    if( profiling ) cerr <<"mnbrak: " << nbValueEvaluations-nbfunc << " value evaluations, " << nbGradientEvaluations-nbgrad << " gradient evaluations" << endl;
    
}

template <class T>
inline T 
MAX( const T& a, const T& b ) { return a > b ? a : b; }

template <class T>
inline T 
SIGN(const T& a, const T& b) { return ((b) >= 0.0 ? Rfabs(a) : -Rfabs(a)); }


Real dbrent (
    Real ax, Real bx, Real cx, 
    //Real fb, Real derivb,        // plus utilises
    Real tol, 
    Real *xmin,
    Real ZEPS = 1.0e-10
)
{
    #define MOV3(a,b,c, d,e,f) (a)=(d);(b)=(e);(c)=(f);
    int ok1,ok2;
    Real a,b,d,d1,d2,du,dv,dw,dx,e=0.0;
    Real fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,v,w,x,xm;

    int nbfunc=nbValueEvaluations;     // for profiling
    int nbgrad=nbGradientEvaluations;

    *xmin = bx; // in case no minimum is found
    a=(ax < cx ? ax : cx);
    b=(ax > cx ? ax : cx);
    x=w=v=bx;
    fw=fv=fx=f1dim(x);
    dw=dv=dx=df1dim(x);
    //fw=fv=fx=fb;
    //dw=dv=dx=derivb;
    Real defaultValue=fx;
    for (int iter=1;iter<=linminMaxIter;iter++) 
    {
        xm=0.5*(a+b);
        tol1=tol*Rfabs(x)+ZEPS;
        tol2=2.0*tol1;
        if (Rfabs(x-xm) <= (tol2-0.5*(b-a))) {
            *xmin=x;    
            if( profiling ) cerr <<"dbrent: " << nbValueEvaluations-nbfunc << " value evaluations, " << nbGradientEvaluations-nbgrad << " gradient evaluations" << endl;
            return fx;
        }
        
        if (Rfabs(e) > tol1) {
            d1=2.0*(b-a);
            d2=d1;
            if (dw != dx) d1=(w-x)*dx/(dx-dw);
            if (dv != dx) d2=(v-x)*dx/(dx-dv);
            u1=x+d1;
            u2=x+d2;
            ok1 = (a-u1)*(u1-b) > 0.0 && dx*d1 <= 0.0;
            ok2 = (a-u2)*(u2-b) > 0.0 && dx*d2 <= 0.0;
            olde=e;
            e=d;
            if (ok1 || ok2) {
                if (ok1 && ok2)
                    d=(Rfabs(d1) < Rfabs(d2) ? d1 : d2);
                else if (ok1)
                    d=d1;
                else
                    d=d2;
                if (Rfabs(d) <= Rfabs(0.5*olde)) {
                    u=x+d;
                    if (u-a < tol2 || b-u < tol2)
                        d=SIGN(tol1,xm-x);
                } else {
                    d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
                }
            } else {
                d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
            }
        } else {
            d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
        }
        
        if (Rfabs(d) >= tol1) {
            u=x+d;
            fu=f1dim(u);
        } else {
            u=x+SIGN(tol1,d);
            fu=f1dim(u);
            if (fu > fx) {
                *xmin=x;
                if( profiling ) cerr <<"dbrent: " << nbValueEvaluations-nbfunc << " value evaluations, " << nbGradientEvaluations-nbgrad << " gradient evaluations" << endl;
                return fx;
            }
        }
        
        du=df1dim(u);
        if (fu <= fx) {
            if (u >= x) a=x; else b=x;
            MOV3(v,fv,dv, w,fw,dw)
            MOV3(w,fw,dw, x,fx,dx)
            MOV3(x,fx,dx, u,fu,du)
        } else {
            if (u < x) a=u; else b=u;
            if (fu <= fw || w == x) {
                MOV3(v,fv,dv, w,fw,dw)
                MOV3(w,fw,dw, u,fu,du)
            } else if (fu < fv || v == x || v == w) {
                MOV3(v,fv,dv, u,fu,du)
            }
        }
    }
    if( profiling ) cerr <<"Too many iterations in routine dbrent: " << nbValueEvaluations-nbfunc << " value evaluations, " << nbGradientEvaluations-nbgrad << " gradient evaluations" << endl;
    tooManyIterationsDbrent = true;
    return defaultValue;
    #undef MOV3
}

void dumpPotential ( Real stepMin, Real stepMax, int resolution, const char* filename = "potential.plo" )
{
    ofstream out(filename);
    if(!out) {
        cerr << "could not open " << filename << "for writing" << endl;
        return;
    }
    
    Real delta = (stepMax-stepMin)/resolution;
    for( Real x=stepMin; x<stepMax; x+=delta ){
        out << x << "\t" << f1dim(x) << endl;
    }
    valueFunction.display(); // to reset the counter
}

Real Rfabs( Real a ) { return a<0 ? -a : a; }


    //  Deprecated ! Use operator() instead.
    //  Optimize the parameters in the steepest direction starting from given value(s). Limit the displacement to maxStep. Perform only one step. Dump the one-dimensional potential function to file pot.dmp if the step is smaller than dumpThreshold.
    //  */
    //void oneStep( Parameter& param, Real maxStep, Real dumpThreshold=1.0e-10 )
    //{
    //  init(animal::size(param));
    //  
    //  //oneCubicStep( param, maxStep );
//
    //  (valueFunction)(param);
    //  (gradientFunction)(param,xi);
//
    //  // compute a step in the search direction
    //  Real step = linemin(param,xi,&finalValue,linminTol);
    //  //cout << "oneStep: optimal stepRatio = " << step/maxStep << endl;
    //  if( step < -maxStep ){ // limit reached (step is negative)
    //      step = -maxStep ;
    //  }
//
    //  // apply the step
    //  for (unsigned int j=0;j<size;j++) {
    //      xi[j] *= step;
    //      param[j] += xi[j];
    //  }
    //  
    //  //
    //  if( Rfabs(step)<dumpThreshold ){
    //      char* dumpfile = "pot.dmp";
    //      dumpPotential(-maxStep, maxStep, 100, dumpfile );
    //      cerr << "potential dumped to file " << dumpfile << endl;
    //  }
    //}
//
//
    //  Optimize the parameters in the steepest direction starting from given value(s) by finding the minimu of a Hermite spline defined by the values and derivatives at step 0 and maxStep. Take the sortest time step. Limit the displacement to maxStep.
    //  */
    //void oneCubicStep( const Parameter& param, Real maxStep )
    //{
    //  init(animal::size(param));      
//
    //  pcom = param;
//
    //  Real a0 = (valueFunction)(pcom);
    //  (gradientFunction)(pcom,xi);
    //  animal::v_teq(xi,-1.0);
    //  Real a1 = -animal::v_dot(xi,xi);
    //  
    //  animal::v_peq(pcom,maxStep,xi);      // positions at maxStep
    //  Real y1 = (valueFunction)(pcom); // potential at maxStep
    //  (gradientFunction)(pcom,xt);
    //  animal::v_teq(xt,-1.0);
    //  Real y_1 = -animal::v_dot(xi,xt);  // derivative at maxStep
    //  
    //  Real& s = maxStep;
    //  Real s2 = s*s;
    //  Real s3 = s2*s;
    //  
    //  Real a2 = (-3*a0 -2*s*a1 +3*y1 -s*y_1)/s2;
    //  Real a3 = ( 2*a0 +  s*a1 -2*y1 +s*y_1)/s3;
    //  
//
    //  
    //  
    //  Real step;
    //  Real d = a2*a2 - 3*a1*a3;  // determinant
    //  if( d<0 ){
    //      // no minimum. apply max step
    //      cerr << "no minimum, apply maxStep" << endl;
    //      step = s;
    //  }
    //  else {
    //      d = sqrt(d)/(3*a3);  // +/-
    //      step = -a2/(3*a3) -d;
    //      if( step>s)
    //          step = s;
    //      else if( step<0 ){
    //          step += 2*d;
    //          if( step>s)
    //              step = s;               
    //      }                   
    //  }
    //  //assert( stepRatio >= 0 );
//
    //  // apply the step
    //  //animal::v_peq( positions, maxStep*stepRatio, xi );
    //}
//
//
};

}// animal
#endif

---