ConstrainedOptPack_DirectLineSearchArmQuad_Strategy.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <ostream>
00030 #include <iomanip>
00031 #include <sstream>
00032 
00033 #include "ConstrainedOptPack_DirectLineSearchArmQuad_Strategy.hpp"
00034 #include "ConstrainedOptPack_MeritFuncCalc1D.hpp"
00035 #include "check_nan_inf.h"
00036 
00037 namespace ConstrainedOptPack {
00038 inline value_type min(value_type v1, value_type v2) {
00039   return (v1 < v2) ? v1 : v2;
00040 }
00041 inline value_type max(value_type v1, value_type v2) {
00042   return (v1 > v2) ? v1 : v2;
00043 }
00044 }
00045 
00046 ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::DirectLineSearchArmQuad_Strategy(
00047   int           max_iter
00048   ,value_type   eta
00049   ,value_type   min_frac
00050   ,value_type   max_frac
00051   ,bool         max_out_iter
00052   )
00053   :max_iter_(max_iter)
00054   ,eta_(eta)
00055   ,min_frac_(min_frac)
00056   ,max_frac_(max_frac)
00057   ,max_out_iter_(max_out_iter)
00058 {}
00059 
00060 void ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::set_max_iter(int max_iter)
00061 {
00062   max_iter_ = max_iter;
00063 }
00064 
00065 int ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::max_iter() const
00066 {
00067   return max_iter_;
00068 }
00069 
00070 int ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::num_iterations() const
00071 {
00072   return num_iter_;
00073 }
00074 
00075 bool ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::do_line_search(
00076     const MeritFuncCalc1D& phi, value_type phi_k
00077   , value_type* alpha_k, value_type* phi_kp1
00078   , std::ostream* out )
00079 {
00080   using std::setw;
00081   using std::endl;
00082 
00083   if(*alpha_k < 0.0) {
00084     throw std::logic_error( "DirectLineSearchArmQuad_Strategy::do_line_search(): "
00085       "alpha_k can't start out less than 0.0" );
00086   }
00087 
00088   validate_parameters();
00089   
00090   int w = 20;
00091   int prec = 8;
00092   int prec_saved;
00093   if(out) {
00094     prec_saved = out->precision();
00095     *out  << std::setprecision(prec)
00096         << "\nStarting Armijo Quadratic interpolation linesearch ...\n";
00097   }
00098 
00099   // Loop initialization (technically the first iteration)
00100 
00101   const value_type Dphi_k = phi.deriv();
00102 
00103   // output header
00104   if(out) {
00105     if(max_out_iter())
00106       *out << "\nmax_out_iter == true, maxing out the number of iterations!"; 
00107     *out  << "\nDphi_k = "  << Dphi_k
00108         << "\nphi_k = "   << phi_k << "\n\n"
00109         << setw(5)      << "itr"
00110         << setw(w)      << "alpha_k"
00111         << setw(w)      << "phi_kp1"
00112         << setw(w)      << "phi_kp1-frac_phi\n"
00113         << setw(5)      << "----"
00114         << setw(w)      << "----------------"
00115         << setw(w)      << "----------------"
00116         << setw(w)      << "----------------\n";
00117   }
00118 
00119   // Check that this is a decent direction
00120   if(Dphi_k >= 0) throw DirectLineSearch_Strategy::NotDescentDirection(
00121     "DirectLineSearchArmQuad_Strategy::do_line_search(): "
00122     "The given d_k is not a descent direction for the given "
00123     "phi (phi.deriv() is >= 0)"     );
00124 
00125   // keep memory of the best value
00126   value_type  best_alpha = *alpha_k, best_phi = *phi_kp1;
00127 
00128   // Perform linesearch.
00129   bool success = false;
00130   for( num_iter_ = 0; num_iter_ < max_iter(); ++num_iter_ ) {
00131 
00132     // Print out this iteration.
00133 
00134     value_type frac_phi = phi_k + eta() * (*alpha_k) * Dphi_k;
00135     if(out)
00136       *out  << setw(5)      << num_iter_
00137           << setw(w)      << *alpha_k
00138           << setw(w)      << *phi_kp1
00139           << setw(w)      << ((*phi_kp1)-frac_phi)  << endl;
00140     
00141     // Check that this is a valid number.
00142     if( RTOp_is_nan_inf( *phi_kp1 ) ) {
00143       // Cut back the step to min_frac * alpha_k
00144       *alpha_k = min_frac()*(*alpha_k);
00145       best_alpha = 0.0;
00146       best_phi = phi_k;
00147     }
00148     else {    
00149 
00150       // Armijo condition
00151       if( *phi_kp1 < frac_phi ) {
00152         // We have found an acceptable point
00153         success = true;
00154         if( !max_out_iter() || ( max_out_iter() && num_iter_ == max_iter() - 1 ) )
00155           break;  // get out of the loop
00156       }
00157 
00158       // Select a new alpha to try:
00159       //   alpha_k = ( min_frac*alpha_k <= quadratic interpolation <= max_frac*alpha_k )
00160 
00161       // Quadratic interpolation of alpha_k that minimizes phi.
00162       // We know the values of phi at the initail point and alpha_k and
00163       // the derivate of phi w.r.t alpha at the initial point and
00164       // that's enough information for a quandratic interpolation.
00165       
00166       value_type alpha_quad =   ( -0.5 * Dphi_k * (*alpha_k) * (*alpha_k) )
00167                     / ( (*phi_kp1) - phi_k - (*alpha_k) * Dphi_k );
00168 
00169       *alpha_k = min( max(min_frac()*(*alpha_k),alpha_quad), max_frac()*(*alpha_k) );
00170 
00171     }
00172     
00173     // Evaluate the point
00174 
00175     *phi_kp1 = phi(*alpha_k);
00176 
00177     // Save the best point found
00178     if(*phi_kp1 < best_phi) {
00179       best_phi = *phi_kp1;
00180       best_alpha = *alpha_k;
00181     }
00182 
00183   }
00184 
00185   // Be nice and reset the precision
00186   if(out) {
00187     out->precision(prec_saved);
00188   }
00189 
00190   if( success ) {
00191     return true;
00192   }
00193 
00194   // Line search failure.  Return the best point found and let the 
00195   // client decide what to do.
00196   *alpha_k = best_alpha;
00197   *phi_kp1 = phi(best_alpha); // Make this the last call to phi(x)
00198   return false; 
00199 }
00200 
00201 void ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::print_algorithm(
00202   std::ostream& out, const std::string& L) const
00203 {
00204   out
00205     << L << "*** start line search using the Armijo cord test and quadratic interpolation of alpha\n"
00206     << L << "default: max_ls_iter = " << max_iter() << std::endl
00207     << L << "         eta = " << eta() << std::endl
00208     << L << "         min_frac = " << min_frac() << std::endl
00209     << L << "         max_frac = " << max_frac() << std::endl
00210     << L << "         max_out_iter = " << max_out_iter() << std::endl
00211     << L << "Dphi_k = phi.deriv()\n"
00212     << L << "if Dphi_k >= 0\n"
00213     << L << "    throw not_descent_direction\n"
00214     << L << "    end line search\n"
00215     << L << "end\n"
00216     << L << "best_alpha = alpha_k\n"
00217     << L << "best_phi = phi_kp1\n"
00218     << L << "for num_iter = 0... max_ls_iter\n"
00219     << L << "    frac_phi = phi_k + eta * alpha_k * Dphi_k\n"
00220     << L << "    print iteration\n"
00221     << L << "    if phi_kp1 is not a valid number then\n"
00222     << L << "        *** Cut back the step so the NLP's functions\n"
00223     << L << "        *** will hopefully be defined.\n"
00224     << L << "        alpha_k = min_frac * alpha_k\n"
00225     << L << "        best_alpha = 0\n"
00226     << L << "        best_phi = phi_k\n"
00227     << L << "    else\n"
00228     << L << "        if ( phi_kp1 < frac_phi ) then\n"
00229     << L << "            *** We have found an acceptable point\n"
00230     << L << "            if( !max_out_iter || max_out_iter && num_iter == max_ls_iter - 1 ) )\n"
00231     << L << "                end line search\n"
00232     << L << "            end\n"
00233     << L << "        end\n"
00234     << L << "        *** Use a quadratic interpoation to minimize phi(alpha)\n"
00235     << L << "        alpha_quad = (-0.5 * Dphi_k * alpha_k^2) / ( phi_kp1 - phi_k - alpha_k*Dphi_k )\n"
00236     << L << "        alpha_k = min( max( min_frac*alpha_k, alpha_quad ), max_frac*alpha_k )\n"
00237     << L << "    end\n"
00238     << L << "    phi_kp1 = phi(alpha_k)\n"
00239     << L << "    if phi_kp1 < best_phi\n"
00240     << L << "        best_phi = phi_kp1\n"
00241     << L << "        best_alpha = alpha_k\n"
00242     << L << "    end\n"
00243     << L << "end\n"
00244     << L << "*** If you get there the line search failed.\n"
00245     << L << "alpha_k = best_alpha\n"
00246     << L << "phi_kp1 = phi(alpha_k)\n"
00247     << L << "linesearch_failure = true\n";
00248 }
00249 
00250 void ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::validate_parameters() const
00251 {
00252   if( eta() < 0.0 || 1.0 < eta() ) {
00253     std::ostringstream omsg;
00254     omsg
00255       << "DirectLineSearchArmQuad_Strategy::validate_parameters() : "
00256       << "Error, eta = " << eta() << " is not in the range [0, 1]";
00257     throw std::invalid_argument( omsg.str() );
00258   }
00259   if( min_frac() < 0.0 || 1.0 < min_frac() ) {
00260     std::ostringstream omsg;
00261     omsg
00262       << "DirectLineSearchArmQuad_Strategy::validate_parameters() : "
00263       << "Error, min_frac = " << min_frac() << " is not in the range [0, 1]";
00264     throw std::invalid_argument( omsg.str() );
00265   }
00266   if( max_frac() < 0.0 || ( !max_out_iter() && 1.0 < max_frac() ) ) {
00267     std::ostringstream omsg;
00268     omsg
00269       << "DirectLineSearchArmQuad_Strategy::validate_parameters() : "
00270       << "Error, max_frac = " << max_frac() << " is not in the range [0, 1]";
00271     throw std::invalid_argument( omsg.str() );
00272   }
00273   if( max_frac() < min_frac() ) {
00274     std::ostringstream omsg;
00275     omsg
00276       << "DirectLineSearchArmQuad_Strategy::validate_parameters() : "
00277       << "Error, max_frac = " << max_frac()
00278       << " < min_frac = " << min_frac();;
00279     throw std::invalid_argument( omsg.str() );
00280   }
00281 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:10:57 2011 for MOOCHO (Single Doxygen Collection) by  doxygen 1.6.3