00001 /***************************************************************************
00002                           linesearch.h  -  description
00003                              -------------------
00004     begin                : Tue May 2 2000
00005     copyright            : (C) 2000 by Michael Peeters
00006     email                : Michael.Peeters@vub.ac.be
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017 
00018 #ifndef LINESEARCH_H
00019 #define LINESEARCH_H
00020 
00021 #include "numerictypes.h"
00022 #include "normfunction.h"
00023 #include "utility.h"
00024 #include <stdexcept>
00025 // #include "debugmacro.h"
00026 
00027 namespace MODEL{
00028 
00032   template <integer dims, class NT = NumericTraits<number,dims> >
00033   class LineSearch // Could we derive it off vectorfunction - that is what it is... ???
00034   {
00035   public:
00036         typedef typename NT::number                     numT;
00037         typedef typename NT::vect                       vect;   
00038         typedef NormFunction<dims,NT>           nf;
00039         typedef typename NT::vf                         vf;
00040                 
00041         LineSearch(vf& fu, bool own=false)
00042           :  func(own?fu.clone():&fu), fmin(nf(*func,false,0.5)), tooclose(false), owned(own) {};
00043         ~LineSearch() {if(owned)delete func;}
00044                 
00056         void operator()                 (       const vect& uold, const numT fold,
00057                                                                 vect& grad, vect& p,
00058                                                                 vect& u, numT& f, numT maxstep);
00059 
00060         bool converged(void) {return tooclose;}
00061                         
00064         nf&             norm(void){return fmin;}
00065         
00066         LineSearch(const LineSearch& ls) :      func(ls.func), fmin(ls.fmin), tooclose(ls.tooclose),
00067                                                                                 owned(false) {}
00068         // to be safe
00069         PRIVATE_ASSIGN(LineSearch);
00070         
00071 
00072   private:
00073         vf*             func;
00074         nf              fmin;
00075         bool    tooclose;
00076         bool    owned;
00077         
00078   public:
00080         static const number     alfa=1e-4;      
00082         static const number     tolerance=1e-7;
00083   };
00084 
00085   template <integer dims, class NT>
00086   void LineSearch<dims,NT>::operator()                  
00087         (       const vect& uold, const numT fold,
00088                 vect& grad, vect& p,
00089                 vect& u, numT& f, numT maxstep)
00090   {
00091         tooclose=false;                         // in NRC called: check
00092         
00093         // rescale if too large
00094         number  pl(length(p)); // in NRC called: sum
00095         if(pl>maxstep) p*=maxstep/pl;   
00096 
00097         // scalar poduct of the gradient and the direction      
00098         number  slope(grad*p);          
00099         
00100         // calculate minimal lambda using heuristic
00101         number  lambdascale(0.0);       // in NRC called: test
00102         for (integer i=0;i<dims;i++)
00103           {
00104                 number temp;
00105                 temp=abs(p[i])/max(abs(uold[i]),number(1.0));
00106                 if (temp > lambdascale) lambdascale=temp;
00107           }
00108         
00109         number  lambdamin(tolerance/lambdascale),lambda(1.);            // in NRC called alamin,alam
00110         
00111         number lambda2(0.),templambda(0.),f2(0.);
00112         while (true) // Search forever
00113           {
00114                         
00115                 u=uold; u+=lambda*p;
00116                 f=fmin(u);
00117                         
00118                 //              DEBUG_Microscopic << "linesearch:" << u << " | " << f << " | "
00119                 //                        << lambda*p <<endl;
00120                                                 
00121                 if (lambda < lambdamin) { tooclose=true; u=uold ;return;}       // wooha - we are on top of it
00122                 else if (f <= fold+alfa*lambda*slope) return;                   // step is OK (hopefully this happens)
00123                 else
00124                   {
00125                         // First step ? We use a quadratic model for lambda
00126                         if (lambda == 1.0) templambda = -slope/(2.0*(f-fold-slope));
00127                         // All other steps: use a cubic (a few more calculations)
00128                         else
00129                           {     
00130                                 number  rhs1 = f-fold-lambda*slope;
00131                                 number  rhs2 = f2-fold-lambda2*slope;
00132                                         
00133                                 number  a=(rhs1/(lambda*lambda)-rhs2/(lambda2*lambda2))/(lambda-lambda2);
00134                                 number  b=(-lambda2*rhs1/(lambda*lambda)
00135                                                    +lambda*rhs2/(lambda2*lambda2))/(lambda-lambda2);
00136                                         
00137                                 if (a == 0.0) templambda = -slope/(2.0*b); // easy roots
00138                                 else
00139                                   {     // use discriminant
00140                                         number  disc=b*b-3.0*a*slope;
00141                                         if (disc<0.0) templambda=0.5*lambda;
00142                                         else if (b <= 0.0) templambda=(-b+sqrt(disc))/(3.0*a);
00143                                         else templambda=-slope/(b+sqrt(disc));                                                                                                  
00144                                   }
00145                                 if (templambda>0.5*lambda) templambda=0.5*lambda; // don't overdo it.
00146                           }     
00147                   }
00148                 lambda2=lambda;
00149                 f2 = f;
00150                 lambda=max(templambda,0.1*lambda);      // don't overdo it.
00151           }
00152   }
00153 
00154 } // end MODEL
00155 #endif
00156 
00157                 
00158 
00159 /*********************************************************************
00160 $Id: linesearch.h,v 1.1 2001/05/22 10:54:55 mpeeters Exp $
00161 **********************************************************************
00162 
00163 $Log: linesearch.h,v $
00164 Revision 1.1  2001/05/22 10:54:55  mpeeters
00165 Moved sources and headers for libModel to model/ subdirectory, in an attempt to rationalize the source tree. This should make things "netter".
00166 
00167 Revision 1.2  2000/09/15 10:26:31  mpeeters
00168 Cleaned out header and added CVS tails to files, removed superfuous
00169 @author comments, inserted dates
00170 
00171 *********************************************************************/