ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
ConstrainedOptPack_DirectLineSearchArmQuad_Strategy.cpp
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <ostream>
00043 #include <iomanip>
00044 #include <sstream>
00045 
00046 #include "ConstrainedOptPack_DirectLineSearchArmQuad_Strategy.hpp"
00047 #include "ConstrainedOptPack_MeritFuncCalc1D.hpp"
00048 #include "check_nan_inf.h"
00049 
00050 namespace ConstrainedOptPack {
00051 inline value_type min(value_type v1, value_type v2) {
00052   return (v1 < v2) ? v1 : v2;
00053 }
00054 inline value_type max(value_type v1, value_type v2) {
00055   return (v1 > v2) ? v1 : v2;
00056 }
00057 }
00058 
00059 ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::DirectLineSearchArmQuad_Strategy(
00060   int           max_iter
00061   ,value_type   eta
00062   ,value_type   min_frac
00063   ,value_type   max_frac
00064   ,bool         max_out_iter
00065   )
00066   :max_iter_(max_iter)
00067   ,eta_(eta)
00068   ,min_frac_(min_frac)
00069   ,max_frac_(max_frac)
00070   ,max_out_iter_(max_out_iter)
00071 {}
00072 
00073 void ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::set_max_iter(int max_iter)
00074 {
00075   max_iter_ = max_iter;
00076 }
00077 
00078 int ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::max_iter() const
00079 {
00080   return max_iter_;
00081 }
00082 
00083 int ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::num_iterations() const
00084 {
00085   return num_iter_;
00086 }
00087 
00088 bool ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::do_line_search(
00089     const MeritFuncCalc1D& phi, value_type phi_k
00090   , value_type* alpha_k, value_type* phi_kp1
00091   , std::ostream* out )
00092 {
00093   using std::setw;
00094   using std::endl;
00095 
00096   if(*alpha_k < 0.0) {
00097     throw std::logic_error( "DirectLineSearchArmQuad_Strategy::do_line_search(): "
00098       "alpha_k can't start out less than 0.0" );
00099   }
00100 
00101   validate_parameters();
00102   
00103   int w = 20;
00104   int prec = 8;
00105   int prec_saved;
00106   if(out) {
00107     prec_saved = out->precision();
00108     *out  << std::setprecision(prec)
00109         << "\nStarting Armijo Quadratic interpolation linesearch ...\n";
00110   }
00111 
00112   // Loop initialization (technically the first iteration)
00113 
00114   const value_type Dphi_k = phi.deriv();
00115 
00116   // output header
00117   if(out) {
00118     if(max_out_iter())
00119       *out << "\nmax_out_iter == true, maxing out the number of iterations!"; 
00120     *out  << "\nDphi_k = "  << Dphi_k
00121         << "\nphi_k = "   << phi_k << "\n\n"
00122         << setw(5)      << "itr"
00123         << setw(w)      << "alpha_k"
00124         << setw(w)      << "phi_kp1"
00125         << setw(w)      << "phi_kp1-frac_phi\n"
00126         << setw(5)      << "----"
00127         << setw(w)      << "----------------"
00128         << setw(w)      << "----------------"
00129         << setw(w)      << "----------------\n";
00130   }
00131 
00132   // Check that this is a decent direction
00133   if(Dphi_k >= 0) throw DirectLineSearch_Strategy::NotDescentDirection(
00134     "DirectLineSearchArmQuad_Strategy::do_line_search(): "
00135     "The given d_k is not a descent direction for the given "
00136     "phi (phi.deriv() is >= 0)"     );
00137 
00138   // keep memory of the best value
00139   value_type  best_alpha = *alpha_k, best_phi = *phi_kp1;
00140 
00141   // Perform linesearch.
00142   bool success = false;
00143   for( num_iter_ = 0; num_iter_ < max_iter(); ++num_iter_ ) {
00144 
00145     // Print out this iteration.
00146 
00147     value_type frac_phi = phi_k + eta() * (*alpha_k) * Dphi_k;
00148     if(out)
00149       *out  << setw(5)      << num_iter_
00150           << setw(w)      << *alpha_k
00151           << setw(w)      << *phi_kp1
00152           << setw(w)      << ((*phi_kp1)-frac_phi)  << endl;
00153     
00154     // Check that this is a valid number.
00155     if( RTOp_is_nan_inf( *phi_kp1 ) ) {
00156       // Cut back the step to min_frac * alpha_k
00157       *alpha_k = min_frac()*(*alpha_k);
00158       best_alpha = 0.0;
00159       best_phi = phi_k;
00160     }
00161     else {    
00162 
00163       // Armijo condition
00164       if( *phi_kp1 < frac_phi ) {
00165         // We have found an acceptable point
00166         success = true;
00167         if( !max_out_iter() || ( max_out_iter() && num_iter_ == max_iter() - 1 ) )
00168           break;  // get out of the loop
00169       }
00170 
00171       // Select a new alpha to try:
00172       //   alpha_k = ( min_frac*alpha_k <= quadratic interpolation <= max_frac*alpha_k )
00173 
00174       // Quadratic interpolation of alpha_k that minimizes phi.
00175       // We know the values of phi at the initail point and alpha_k and
00176       // the derivate of phi w.r.t alpha at the initial point and
00177       // that's enough information for a quandratic interpolation.
00178       
00179       value_type alpha_quad =   ( -0.5 * Dphi_k * (*alpha_k) * (*alpha_k) )
00180                     / ( (*phi_kp1) - phi_k - (*alpha_k) * Dphi_k );
00181 
00182       *alpha_k = min( max(min_frac()*(*alpha_k),alpha_quad), max_frac()*(*alpha_k) );
00183 
00184     }
00185     
00186     // Evaluate the point
00187 
00188     *phi_kp1 = phi(*alpha_k);
00189 
00190     // Save the best point found
00191     if(*phi_kp1 < best_phi) {
00192       best_phi = *phi_kp1;
00193       best_alpha = *alpha_k;
00194     }
00195 
00196   }
00197 
00198   // Be nice and reset the precision
00199   if(out) {
00200     out->precision(prec_saved);
00201   }
00202 
00203   if( success ) {
00204     return true;
00205   }
00206 
00207   // Line search failure.  Return the best point found and let the 
00208   // client decide what to do.
00209   *alpha_k = best_alpha;
00210   *phi_kp1 = phi(best_alpha); // Make this the last call to phi(x)
00211   return false; 
00212 }
00213 
00214 void ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::print_algorithm(
00215   std::ostream& out, const std::string& L) const
00216 {
00217   out
00218     << L << "*** start line search using the Armijo cord test and quadratic interpolation of alpha\n"
00219     << L << "default: max_ls_iter = " << max_iter() << std::endl
00220     << L << "         eta = " << eta() << std::endl
00221     << L << "         min_frac = " << min_frac() << std::endl
00222     << L << "         max_frac = " << max_frac() << std::endl
00223     << L << "         max_out_iter = " << max_out_iter() << std::endl
00224     << L << "Dphi_k = phi.deriv()\n"
00225     << L << "if Dphi_k >= 0\n"
00226     << L << "    throw not_descent_direction\n"
00227     << L << "    end line search\n"
00228     << L << "end\n"
00229     << L << "best_alpha = alpha_k\n"
00230     << L << "best_phi = phi_kp1\n"
00231     << L << "for num_iter = 0... max_ls_iter\n"
00232     << L << "    frac_phi = phi_k + eta * alpha_k * Dphi_k\n"
00233     << L << "    print iteration\n"
00234     << L << "    if phi_kp1 is not a valid number then\n"
00235     << L << "        *** Cut back the step so the NLP's functions\n"
00236     << L << "        *** will hopefully be defined.\n"
00237     << L << "        alpha_k = min_frac * alpha_k\n"
00238     << L << "        best_alpha = 0\n"
00239     << L << "        best_phi = phi_k\n"
00240     << L << "    else\n"
00241     << L << "        if ( phi_kp1 < frac_phi ) then\n"
00242     << L << "            *** We have found an acceptable point\n"
00243     << L << "            if( !max_out_iter || max_out_iter && num_iter == max_ls_iter - 1 ) )\n"
00244     << L << "                end line search\n"
00245     << L << "            end\n"
00246     << L << "        end\n"
00247     << L << "        *** Use a quadratic interpoation to minimize phi(alpha)\n"
00248     << L << "        alpha_quad = (-0.5 * Dphi_k * alpha_k^2) / ( phi_kp1 - phi_k - alpha_k*Dphi_k )\n"
00249     << L << "        alpha_k = min( max( min_frac*alpha_k, alpha_quad ), max_frac*alpha_k )\n"
00250     << L << "    end\n"
00251     << L << "    phi_kp1 = phi(alpha_k)\n"
00252     << L << "    if phi_kp1 < best_phi\n"
00253     << L << "        best_phi = phi_kp1\n"
00254     << L << "        best_alpha = alpha_k\n"
00255     << L << "    end\n"
00256     << L << "end\n"
00257     << L << "*** If you get there the line search failed.\n"
00258     << L << "alpha_k = best_alpha\n"
00259     << L << "phi_kp1 = phi(alpha_k)\n"
00260     << L << "linesearch_failure = true\n";
00261 }
00262 
00263 void ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::validate_parameters() const
00264 {
00265   if( eta() < 0.0 || 1.0 < eta() ) {
00266     std::ostringstream omsg;
00267     omsg
00268       << "DirectLineSearchArmQuad_Strategy::validate_parameters() : "
00269       << "Error, eta = " << eta() << " is not in the range [0, 1]";
00270     throw std::invalid_argument( omsg.str() );
00271   }
00272   if( min_frac() < 0.0 || 1.0 < min_frac() ) {
00273     std::ostringstream omsg;
00274     omsg
00275       << "DirectLineSearchArmQuad_Strategy::validate_parameters() : "
00276       << "Error, min_frac = " << min_frac() << " is not in the range [0, 1]";
00277     throw std::invalid_argument( omsg.str() );
00278   }
00279   if( max_frac() < 0.0 || ( !max_out_iter() && 1.0 < max_frac() ) ) {
00280     std::ostringstream omsg;
00281     omsg
00282       << "DirectLineSearchArmQuad_Strategy::validate_parameters() : "
00283       << "Error, max_frac = " << max_frac() << " is not in the range [0, 1]";
00284     throw std::invalid_argument( omsg.str() );
00285   }
00286   if( max_frac() < min_frac() ) {
00287     std::ostringstream omsg;
00288     omsg
00289       << "DirectLineSearchArmQuad_Strategy::validate_parameters() : "
00290       << "Error, max_frac = " << max_frac()
00291       << " < min_frac = " << min_frac();;
00292     throw std::invalid_argument( omsg.str() );
00293   }
00294 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends