[Toon-members] TooN/optimization conjugate_gradient.h

[Edward Rosten]

CVSROOT:        /cvsroot/toon
Module name:    TooN
Changes by:     Edward Rosten <edrosten>        10/01/21 22:45:40

Modified files:
        optimization   : conjugate_gradient.h 

Log message:
        Add some NaN-proofing to conjugate gradient optimizer. 
        
        Still no max iterations for the bracketing step.

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/optimization/conjugate_gradient.h?cvsroot=toon&r1=1.9&r2=1.10

Patches:
Index: conjugate_gradient.h
===================================================================
RCS file: /cvsroot/toon/TooN/optimization/conjugate_gradient.h,v
retrieving revision 1.9
retrieving revision 1.10
diff -u -b -r1.9 -r1.10
--- conjugate_gradient.h        19 Nov 2009 18:19:02 -0000      1.9
+++ conjugate_gradient.h        21 Jan 2010 22:45:38 -0000      1.10
@@ -60,14 +60,31 @@
                b = lambda;
                b_val = func(b);
 
+               while(std::isnan(b_val))
+               {
+                       //We've probably gone in to an invalid region. This can 
happen even 
+                       //if following the gradient would never get us there.
+                       //try backing off lambda
+                       lambda*=.5;
+                       b = lambda;
+                       b_val = func(b);
+
+               }
+
+
                if(b_val < a_val) //We've gone downhill, so keep searching 
until we go back up
                {
+                       double last_good_lambda = lambda;
+                       
                        for(;;)
                        {
                                lambda *= 2;
                                c = lambda;
                                c_val = func(c);
 
+                               if(std::isnan(c_val))
+                                       break;
+                               last_good_lambda = lambda;
                                if(c_val >      b_val) // we have a bracket
                                        break;
                                else
@@ -79,6 +96,27 @@
 
                                }
                        }
+
+                       //We took a step too far.
+                       //Back up: this will not attempt to ensure a bracket
+                       if(std::isnan(c_val))
+                       {
+                               double bad_lambda=lambda;
+                               double l=1;
+
+                               for(;;)
+                               {
+                                       l*=.5;
+                                       c = last_good_lambda + (bad_lambda - 
last_good_lambda)*l;
+                                       c_val = func(c);
+
+                                       if(!std::isnan(c_val))
+                                               break;
+                               }
+
+
+                       }
+
                }
                else //We've overshot the minimum, so back up
                {
@@ -280,6 +318,17 @@
                        return;
 
                //We should have a bracket here
+
+               if(c < b)
+               {
+                       //Failed to bracket due to NaN, so c is the best known 
point.
+                       //Simply go there.
+                       x-=h * c;
+                       y=c_val;
+
+               }
+               else
+               {
                assert(a < b && b < c);
                assert(a_val > b_val && b_val < c_val);
 
@@ -293,6 +342,7 @@
                x -= m[0] * h;
                y = m[1];
        }
+       }
 
        ///Check to see it iteration should stop. You probably do not want to 
use
        ///this function. See iterate() instead. This function updates nothing.