Changeset 6574

Timestamp:

11/15/10 09:57:00 (14 years ago)

Author:

Mathieu Morlighem

Message:

Extended cm_jump so that whenever f(x)/f(0)<cm_jump, return

File:

• issm/trunk/src/c/shared/Numerics/BrentSearch.cpp (modified) (11 diffs)

issm/trunk/src/c/shared/Numerics/BrentSearch.cpp

@@ -20 +20 @@
          * (Golden or parabolic procedure)*/
 
-        /*optimization variable: */
-        double si;
-        double gold;
-        double intervalgold;
-        double oldintervalgold;
+        /*Intermediary*/
+        double si,gold,intervalgold,oldintervalgold;
         double parab_num,parab_den;
-        double distance;
-
-        /*function values: */
+        double distance,cm_jump;
         double fxmax,fxmin,fxbest;
         double fx,fx1,fx2;
-
-        /*x : */
         double xmax,xmin,xbest;
         double x,x1,x2,xm;
-
-        /*tolerances: */
         double tol1,tol2,seps,tolerance;
-        int    maxiter;
-        double cm_jump;
-
-        /*counters: */
-        int iter,goldenflag,loop;
+        int    maxiter,iter;
+        bool   loop,goldenflag;
 
         /*Recover parameters:*/
@@ -51 +39 @@
         cm_jump=optpars->cm_jump;
         
-        //initialize counter and boundaries
+        /*initialize counter and get response at the boundaries*/
         iter=0;
-
-        //get the value of the function at the first boundary
         fxmin = (*f)(xmin,optargs);
         if isnan(fxmin) _error_("Function evaluation returned NaN");
-
-        //display result
         _printf_(VerboseControl(),"\n        Iteration         x           f(x)       Tolerance         Procedure\n\n");
         _printf_(VerboseControl(),"        %s    %12.6g  %12.6g  %s","   N/A",xmin,fxmin,"         N/A         boundary\n");
-
-        //get the value of the function at the first boundary xmax and display result
         fxmax = (*f)(xmax,optargs);
         if isnan(fxmax) _error_("Function evaluation returned NaN");
         _printf_(VerboseControl(),"        %s    %12.6g  %12.6g  %s","   N/A",xmax,fxmax,"         N/A         boundary\n");
 
-        //test if jump option activated and xmin==0
-        if (!isnan(cm_jump) && (xmin==0)){
-
-                //test ration between fxmin and fxmax. if < criterion, return
-                if((fxmax/fxmin)<cm_jump){
-
-                        /*Assign output pointers: */
-                        *psearch_scalar=xmax;
-                        *pJ=fxmax;
-                        return;
-                }
+        /*test if jump option activated and xmin==0*/
+        if (!isnan(cm_jump) && (xmin==0) && (fxmax/fxmin)<cm_jump){
+                *psearch_scalar=xmax;
+                *pJ=fxmax;
+                return;
         }
 
-        //initialize the other variables
+        /*initialize optimization variables*/
         seps=sqrt(DBL_EPSILON);    //precision of a double
         distance=0.0;              //new_x=old_x + distance
@@ -86 +62 @@
         intervalgold=0.0;          //distance used by Golden procedure
 
-        //Compute initial point
-        
-        //1: initialize the value of the 4 x needed (x1,x2,x,xbest)
+        /*1: initialize the values of the 4 x needed (x1,x2,x,xbest)*/
         x1=xmin+gold*(xmax-xmin);
         x2=x1;
@@ -94 +68 @@
         x=xbest;
 
-        //2: call the function to be evaluated
+        /*2: call the function to be evaluated*/
         fxbest = (*f)(x,optargs);
-        if isnan(fxbest) _error_("Function evaluation returned NaN");
+        if(isnan(fxbest)) _error_("Function evaluation returned NaN");
         iter=iter+1;
 
-        //3: update the other variables
+        /*3: update the other variables*/
         fx1=fxbest;
         fx2=fxbest;
-        //xm is always in the middle of a and b
+        /*xm is always in the middle of a and b*/
         xm=0.5*(xmin+xmax);                           
-        //update tolerances
+        /*update tolerances*/
         tol1=seps*sqrt(pow(xbest,2))+tolerance/3.0;
         tol2=2.0*tol1;
 
-        //4: print result
+        /*4: print result*/
         _printf_(VerboseControl(),"         %5i    %12.6g  %12.6g  %12.6g  %s\n",iter,xbest,fxbest,pow(pow(xbest-xm,2),0.5),"       initial");
 
-        //Main Loop
-        loop=1;
+        loop=true;
         while(loop){
 
-                goldenflag=1;
+                goldenflag=true;
 
                 // Is a parabolic fit possible ?
@@ -121 +94 @@
 
                         // Yes, so fit parabola
-                        goldenflag=0;
+                        goldenflag=false;
                         parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);;
                         parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1);
@@ -145 +118 @@
                                 // f must not be evaluated too close to min_x or max_x
                                 if (((x-xmin)<tol2) || ((xmax-x)<tol2)){
-
-                                        if ((xm-xbest)<0.0){
-                                                si=-1;
-                                        }
-                                        else{
-                                                si=1;
-                                        }
-
+                                        if ((xm-xbest)<0.0) si=-1;
+                                        else                si=1;
                                         //compute new distance
                                         distance=tol1*si;
@@ -158 +125 @@
                         }
                         else{
-
                                 // Not acceptable, must do a golden section step
-                                goldenflag=1;
+                                goldenflag=true;
                         }
                 }
@@ -166 +132 @@
                 //Golden procedure
                 if(goldenflag){
-
                         // compute the new distance d
                         if(xbest>=xm){
@@ -178 +143 @@
 
                 // The function must not be evaluated too close to xbest
-                if(distance<0){
-                        si=-1;
-                }
-                else{
-                        si=1;
-                }
-                if (sqrt(pow(distance,2))>tol1){
-                        x=xbest+si*sqrt(pow(distance,2));
-                }
-                else{
-                        x=xbest+si*tol1;
-                }
+                if(distance<0) si=-1;
+                else           si=1;
+                if(sqrt(pow(distance,2))>tol1) x=xbest+si*sqrt(pow(distance,2));
+                else                           x=xbest+si*tol1;
 
                 //evaluate function on x
                 fx = (*f)(x,optargs);
-                if isnan(fx) _error_("Function evaluation returned NaN");
+                if(isnan(fx)) _error_("Function evaluation returned NaN");
                 iter=iter+1;
 
                 // Update a, b, xm, x1, x2, tol1, tol2
                 if (fx<=fxbest){
-                        if (x>=xbest){
-                                xmin=xbest;
-                        }
-                        else{
-                                xmax=xbest;
-                        }
+                        if (x>=xbest) xmin=xbest;
+                        else          xmax=xbest;
                         x1=x2;    fx1=fx2;
                         x2=xbest; fx2=fxbest;
                         xbest=x;  fxbest=fx;
                 }
-
                 else{ // fx > fxbest
-                        if (x < xbest){
-                                xmin=x;
-                        }
-                        else{
-                                xmax=x;
-                        }
+                        if (x<xbest) xmin=x;
+                        else         xmax=x;
                         if ((fx<=fx2) || (x2==xbest)){
                                 x1=x2; fx1=fx2;
@@ -227 +175 @@
                 tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0;
                 tol2=2.0*tol1;
-
-                //print result
-                if (goldenflag){
-                        _printf_(VerboseControl(),"         %5i    %12.6g  %12.6g  %12.6g  %s\n",iter,x,fx,pow(pow(xbest-xm,2),0.5),"       golden");
-                }
-                else{
-                        _printf_(VerboseControl(),"         %5i    %12.6g  %12.6g  %12.6g  %s\n",iter,x,fx,pow(pow(xbest-xm,2),0.5),"       parabolic");
-                }
-
-                //Stop the optimization?
+                _printf_(VerboseControl(),"         %5i    %12.6g  %12.6g  %12.6g  %s\n",iter,x,fx,pow(pow(xbest-xm,2),0.5),goldenflag?"       golden":"       parabolic");
+
+                /*Stop the optimization?*/
                 if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(xmax-xmin))){
                         _printf_(VerboseControl(),"      %s%g\n","optimization terminated: the current x satisfies the termination criteria using 'tolx' of" ,tolerance);
-                        loop=0;
+                        loop=false;
                 }
                 else if (iter>=maxiter){
                         _printf_(VerboseControl(),"      %s\n","exiting: Maximum number of iterations has been exceeded  - increase 'maxiter'\n");
-                        loop=0;
+                        loop=false;
+                }
+                else if (!isnan(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
+                        _printf_(VerboseControl(),"      %s%g\n","optimization terminated: the current x satisfies the termination criteria using 'cm_jump' of" ,cm_jump);
+                        loop=false;
                 }
                 else{
                         //continue
-                        loop=1;
+                        loop=true;
                 }
         }//end while
@@ -253 +198 @@
         //Now, check that the value on the boundaries are not better than current fxbest
         if (fxbest>fxmin){
-                xbest=optpars->xmin;;
-                fxbest=fxmin;
+                xbest=optpars->xmin; fxbest=fxmin;
         }
         if (fxbest>fxmax){
-                xbest=optpars->xmax;;
-                fxbest=fxmax;
+                xbest=optpars->xmax; fxbest=fxmax;
         }