MoochoPack_LineSearch2ndOrderCorrect_Step.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 <typeinfo>
00031 
00032 #include "MoochoPack_LineSearch2ndOrderCorrect_Step.hpp"
00033 #include "MoochoPack_moocho_algo_conversion.hpp"
00034 #include "IterationPack_print_algorithm_step.hpp"
00035 #include "ConstrainedOptPack_print_vector_change_stats.hpp"
00036 #include "ConstrainedOptPack_MeritFuncCalc1DQuadratic.hpp"
00037 #include "ConstrainedOptPack_MeritFuncCalcNLP.hpp"
00038 #include "ConstrainedOptPack_MeritFuncCalcNLE.hpp"
00039 #include "ConstrainedOptPack_MeritFuncNLESqrResid.hpp"
00040 #include "ConstrainedOptPack/src/VectorWithNorms.h"
00041 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixOp.hpp"
00042 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
00043 #include "AbstractLinAlgPack/src/max_near_feas_step.hpp"
00044 #include "DenseLinAlgPack_DVectorClass.hpp"
00045 #include "DenseLinAlgPack_DVectorOp.hpp"
00046 #include "DenseLinAlgPack_DVectorOut.hpp"
00047 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00048 
00049 namespace LinAlgOpPack {
00050   using AbstractLinAlgPack::Vp_StMtV;
00051 }
00052 
00053 namespace MoochoPack {
00054 
00055 LineSearch2ndOrderCorrect_Step::LineSearch2ndOrderCorrect_Step(
00056   const direct_ls_sqp_ptr_t&      direct_ls_sqp
00057   ,const merit_func_ptr_t&      merit_func
00058   ,const feasibility_step_ptr_t&    feasibility_step
00059   ,const direct_ls_newton_ptr_t&    direct_ls_newton
00060   ,value_type             eta
00061   ,ENewtonOutputLevel         newton_olevel
00062   ,value_type             constr_norm_threshold
00063   ,value_type             constr_incr_ratio
00064   ,int                after_k_iter
00065   ,EForcedConstrReduction       forced_constr_reduction
00066   ,value_type                         forced_reduct_ratio
00067   ,value_type             max_step_ratio
00068   ,int                max_newton_iter
00069   )
00070   :direct_ls_sqp_(direct_ls_sqp)
00071   ,merit_func_(merit_func)
00072   ,feasibility_step_(feasibility_step)
00073   ,direct_ls_newton_(direct_ls_newton)
00074   ,eta_(eta)
00075   ,newton_olevel_(newton_olevel)
00076   ,constr_norm_threshold_(constr_norm_threshold)
00077   ,constr_incr_ratio_(constr_incr_ratio)
00078   ,after_k_iter_(after_k_iter)
00079   ,forced_constr_reduction_(forced_constr_reduction)
00080   ,forced_reduct_ratio_(forced_reduct_ratio)
00081   ,max_step_ratio_(max_step_ratio)
00082   ,max_newton_iter_(max_newton_iter)
00083 {}
00084 
00085 bool LineSearch2ndOrderCorrect_Step::do_step(
00086   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00087   ,poss_type assoc_step_poss
00088   )
00089 {
00090 
00091 /* ToDo: Upate the below code!
00092 
00093   using std::setw;
00094 
00095   using DenseLinAlgPack::dot;
00096   using DenseLinAlgPack::norm_inf;
00097   using DenseLinAlgPack::V_VpV;
00098   using DenseLinAlgPack::V_VmV;
00099   using DenseLinAlgPack::Vp_StV;
00100   using DenseLinAlgPack::Vt_S;
00101 
00102   using LinAlgOpPack::Vp_V;
00103   using LinAlgOpPack::V_MtV;
00104 
00105   using AbstractLinAlgPack::max_near_feas_step;
00106 
00107   using ConstrainedOptPack::print_vector_change_stats;
00108 
00109   typedef LineSearch2ndOrderCorrect_Step  this_t;
00110 
00111   NLPAlgo &algo = rsqp_algo(_algo);
00112   NLPAlgoState  &s    = algo.rsqp_state();
00113   NLP     &nlp  = algo.nlp();
00114 
00115   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00116   std::ostream& out = algo.track().journal_out();
00117   out << std::boolalpha;
00118 
00119   // print step header.
00120   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00121     using IterationPack::print_algorithm_step;
00122     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00123   }
00124 
00125   // /////////////////////////////////////////
00126   // Set references to iteration quantities
00127   //
00128   // Set k+1 first then go back to get k to ensure
00129   // we have backward storage.
00130   
00131   DVector
00132     &x_kp1 = s.x().set_k(+1).v();
00133   value_type
00134     &f_kp1 = s.f().set_k(+1);
00135   DVector
00136     &c_kp1 = s.c().set_k(+1).v();
00137 
00138   const value_type
00139     &f_k = s.f().get_k(0);
00140   const DVector
00141     &c_k = s.c().get_k(0).v();
00142   const DVector
00143     &x_k = s.x().get_k(0).v();
00144   const DVector
00145     &d_k = s.d().get_k(0).v();
00146   value_type
00147     &alpha_k = s.alpha().get_k(0);
00148 
00149   // /////////////////////////////////////
00150   // Compute Dphi_k, phi_kp1 and phi_k
00151 
00152   // Dphi_k
00153   const value_type
00154     Dphi_k = merit_func().deriv();
00155   if( Dphi_k >= 0 ) {
00156     throw LineSearchFailure( "LineSearch2ndOrderCorrect_Step::do_step(...) : " 
00157       "Error, d_k is not a descent direction for the merit function " );
00158   }
00159 
00160   // ph_kp1
00161   value_type
00162     &phi_kp1 = s.phi().set_k(+1) = merit_func().value( f_kp1, c_kp1 );
00163 
00164   // Must compute phi(x) at the base point x_k since the penalty parameter may have changed.
00165   const value_type
00166     &phi_k = s.phi().set_k(0) = merit_func().value( f_k, c_k );
00167 
00168   // //////////////////////////////////////
00169   // Setup the calculation merit function
00170 
00171   // Here f_kp1, and c_kp1 are updated at the same time the
00172   // line search is being performed.
00173   nlp.set_f( &f_kp1 );
00174   nlp.set_c( &c_kp1 );
00175   MeritFuncCalcNLP
00176     phi_calc( &merit_func(), &nlp );
00177 
00178   // //////////////////////////////////////////////////
00179   // Concider 2nd order correction if near solution?
00180   
00181   bool considering_correction = false;
00182   {
00183     const value_type
00184       small_num = std::numeric_limits<value_type>::min(),
00185       nrm_c_k   = s.c().get_k(0).norm_inf(),
00186       nrm_c_kp1 = s.c().get_k(+1).norm_inf();
00187     const bool
00188       test_1 = nrm_c_k <= constr_norm_threshold(),
00189       test_2 = (nrm_c_kp1/(1.0 + nrm_c_k)) < constr_incr_ratio(),
00190       test_3 = s.k() >= after_k_iter();
00191     considering_correction = test_1 && test_2 && test_3;
00192     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00193       out << "\n||c_k||inf = " << nrm_c_k << (test_1 ? " <= " : " > " )
00194         << "constr_norm_threshold = " << constr_norm_threshold()
00195         << "\n||c_kp1||inf/(1.0+||c_k||inf) = "
00196         << "(" << nrm_c_kp1 << ")/(" << 1.0 << " + " << nrm_c_k << ") = "
00197         << ( nrm_c_kp1 / (1.0 + nrm_c_k ) ) << (test_2 ? " <= " : " > " )
00198         << "constr_incr_ratio = " << constr_incr_ratio()
00199         << "\nk = " << s.k() << (test_3 ? " >= " : " < ")
00200         << "after_k_iter = " << after_k_iter()
00201         << (considering_correction
00202           ? ("\nThe computation of a 2nd order correction for x_kp1 = x_k + alpha_k*d_k + alpha_k^2*w"
00203              " will be considered ...\n")
00204           : "\nThe critera for considering a 2nd order correction has not been met ...\n" );
00205     }
00206   }
00207 
00208   // //////////////////////////////
00209   // See if we can take a full step
00210 
00211   bool chose_point = false;
00212 
00213   const value_type frac_phi = phi_k + eta() * Dphi_k;
00214   const bool armijo_test = phi_kp1 <= frac_phi;
00215   if( armijo_test ) {
00216     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00217       out << "\nAccepting full step x_kp1 = x_k + d_k\n";
00218     }
00219     chose_point = true; // The point meets the Armijo test.
00220   }
00221 
00222   // This is storage for function and gradient evaluations for
00223   // the trial newton points and must be remembered for latter
00224   value_type f_xdww;
00225   DVector     c_xdww;
00226   DVector w(x_kp1.size()),    // Full correction after completed computation.
00227        xdww(x_kp1.size());  // Will be set to xdw + sum( w(newton_i), newton_i = 1... )
00228                 // where w(itr) is the local corrections for the current
00229                 // newton iteration.
00230   bool use_correction = false;
00231   bool newton_failed  = true;
00232   
00233   bool considered_correction = ( considering_correction && !chose_point );
00234   if( considered_correction ) {
00235 
00236     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00237       out << "\nConsidering whether to compute a 2nd order correction for\n"
00238           << "x_kp1 = x_k + alpha_k * d_k + alpha_k^2 * w ...\n";
00239     }
00240 
00241     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00242       const value_type obj_descent = dot( s.Gf().get_k(0)(), d_k() );
00243       out << "\nGf_k' * d_k = " << obj_descent << std::endl;
00244       if( obj_descent >= 0.0 ) {
00245         out << "\nWarning, this may not work well with Gf_k'*d_k >= 0.0\n";
00246       }
00247     }
00248 
00249     // Merit function for newton line searches
00250     ConstrainedOptPack::MeritFuncNLESqrResid
00251       phi_c;
00252 
00253     DVector
00254       xdw = x_kp1;  // Will be set to x + d + sum(w(i),i=1..itr-1)
00255               //     where w(i) are previous local corrections
00256     value_type
00257       phi_c_x    = phi_c.value( c_k() ),
00258       phi_c_xd   = phi_c.value( c_kp1() ),
00259       phi_c_xdw  = phi_c_xd,    // No correction is computed yet so w = 0
00260       phi_c_xdww = phi_c_xdw,
00261       nrm_d    = norm_inf( d_k() );
00262 
00263     // Merit function for newton line searches
00264     nlp.set_f( &(f_xdww = f_kp1) );
00265     nlp.set_c( &(c_xdww = c_kp1) );
00266     ConstrainedOptPack::MeritFuncCalcNLE
00267       phi_c_calc( &phi_c, &nlp );
00268 
00269     DVector wy(s.con_decomp().size());  // Range space wy (see latter).
00270 
00271     const bool sufficient_reduction =
00272       phi_c_xd < forced_reduct_ratio() * phi_c_x;
00273     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00274       out << "\nphi_c(c(x_k+d_k)) = " << phi_c_xd << (sufficient_reduction ? " <= " : " > ")
00275         << "forced_reduct_ratio* phi_c(c(x_k)) = " << forced_reduct_ratio() << " * " << phi_c_x
00276         << " = " << (forced_reduct_ratio()*phi_c_x)
00277         << (sufficient_reduction
00278           ? "\nNo need for a 2nd order correciton, perform regular line search ... \n"
00279           : "\nWe need to try to compute a correction w ...\n" );
00280     }
00281     if(sufficient_reduction) {
00282       use_correction = false;
00283     }
00284     else {
00285       // Try to compute a second order correction term.
00286 
00287       // Set print level newton iterations
00288       ENewtonOutputLevel  use_newton_olevel;
00289       if( newton_olevel() == PRINT_USE_DEFAULT ) {
00290         switch(olevel) {
00291             case PRINT_NOTHING:
00292             case PRINT_BASIC_ALGORITHM_INFO:
00293             use_newton_olevel = PRINT_NEWTON_NOTHING;
00294             break;
00295             case PRINT_ALGORITHM_STEPS:
00296             case PRINT_ACTIVE_SET:
00297             use_newton_olevel = PRINT_NEWTON_SUMMARY_INFO;
00298             break;
00299             case PRINT_VECTORS:
00300             use_newton_olevel = PRINT_NEWTON_STEPS;
00301             break;
00302             case PRINT_ITERATION_QUANTITIES:
00303             use_newton_olevel = PRINT_NEWTON_VECTORS;
00304             break;
00305             default:
00306             TEST_FOR_EXCEPT(true);
00307         }
00308       }
00309       else {
00310         use_newton_olevel = newton_olevel();
00311       }
00312       EJournalOutputLevel inner_olevel;
00313       switch(use_newton_olevel) {
00314           case PRINT_NEWTON_NOTHING:
00315           case PRINT_NEWTON_SUMMARY_INFO:
00316           inner_olevel = PRINT_NOTHING;
00317           break;
00318           case PRINT_NEWTON_STEPS:
00319           inner_olevel = PRINT_ALGORITHM_STEPS;
00320           break;
00321           case PRINT_NEWTON_VECTORS:
00322           if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES )
00323             inner_olevel = PRINT_ITERATION_QUANTITIES;
00324           else if( (int)olevel >= (int)PRINT_ACTIVE_SET )
00325             inner_olevel = PRINT_ACTIVE_SET;
00326           else
00327             inner_olevel = PRINT_VECTORS;
00328           break;
00329           default:
00330           TEST_FOR_EXCEPT(true);
00331       }
00332       
00333       // Print header for summary information
00334       const int dbl_min_w = 21;
00335       const int dbl_w = std::_MAX(dbl_min_w,int(out.precision()+8));
00336       if( use_newton_olevel == PRINT_NEWTON_SUMMARY_INFO ) {
00337         out << "\nStarting Newton iterations ...\n"
00338           << "\nphi_c_x   = " << phi_c_x 
00339           << "\nphi_c_xd  = " << phi_c_xd
00340           << "\n||d_k||nf = " << nrm_d << "\n\n"
00341           << setw(5)      << "it"
00342           << setw(dbl_w)    << "||w||inf"
00343           << setw(dbl_w)    << "u"
00344           << setw(dbl_w)    << "step_ratio"
00345           << setw(5)      << "lsit"
00346           << setw(dbl_w)    << "a"
00347           << setw(dbl_w)    << "phi_c_xdww"
00348           << setw(dbl_w)    << "phi_c_xdww-phi_c_x"
00349           << setw(dbl_w)    << "phi_c_xdww-phi_c_xd\n"
00350           << setw(5)      << "----"
00351           << setw(dbl_w)    << "--------------------"
00352           << setw(dbl_w)    << "-------------------"
00353           << setw(dbl_w)    << "-------------------"
00354           << setw(5)      << "----"
00355           << setw(dbl_w)    << "-------------------"
00356           << setw(dbl_w)    << "-------------------"
00357           << setw(dbl_w)    << "-------------------"
00358           << setw(dbl_w)    << "-------------------\n";
00359       }
00360 
00361       // Perform newton feasibility iterations
00362       int newton_i;
00363       for( newton_i = 1; newton_i <= max_newton_iter(); ++newton_i ) {
00364         
00365         if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS ) {
00366           out << "\n**** newton_i = " << newton_i << std::endl;
00367         }
00368 
00369         // Compute a feasibility step
00370         if(!feasibility_step().compute_feasibility_step(
00371           out,inner_olevel,&algo,&s,xdw,nlp.c()(),&w() ))
00372         {
00373           if( (int)use_newton_olevel == (int)PRINT_NEWTON_SUMMARY_INFO ) {
00374             out << "\nCould not compute feasible direction!\n";
00375           }
00376           break; // exit the newton iterations
00377         }
00378 
00379         value_type
00380           nrm_w = norm_inf(w());        
00381 
00382         if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS ) {
00383           out << "\n||w||inf = " << nrm_w << std::endl;
00384         }
00385 
00386         if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_VECTORS ) {
00387           out << "\nw = " << w();
00388         }
00389 
00390         // ////////////////////////////////
00391         // Cutting back w
00392 
00393         value_type a = 1.0; // This is the alpha for your linesearch
00394 
00395         // Cut back w to be in the relaxed bounds.
00396         std::pair<value_type,value_type>
00397           u_steps = max_near_feas_step( s.x().get_k(0)(), w()
00398             , algo.nlp().xl(), algo.nlp().xu()
00399             , algo.algo_cntr().max_var_bounds_viol() );
00400         const value_type u = u_steps.first;
00401 
00402         if( u < a ) {
00403           if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS ) {
00404             out << "\nCutting back w = (a=u) * w to be within relaxed bounds:\n"
00405               << "u = " << u << std::endl;
00406           }
00407           a = u;
00408         }
00409 
00410         // Cut back step so x+d+sum(w(i)) is not too far from x+d
00411         value_type
00412           step_ratio = nrm_w / nrm_d;
00413         if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS ) {
00414           out << "\nstep_ratio = ||w||inf/||d||inf = " << step_ratio
00415               << std::endl;
00416         }
00417         if( a * step_ratio > max_step_ratio() ) {
00418           const value_type aa = a*(max_step_ratio()/step_ratio);
00419           if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS ) {
00420             out << "\na*step_ratio = " << a*step_ratio
00421                 << " > max_step_ratio = " << max_step_ratio() << std::endl
00422               << "Cutting back a = a*max_step_ratio/step_ratio = "
00423                 << aa << std::endl;
00424           }
00425           a = aa;
00426         }
00427 
00428         // /////////////////////////////////////////////////
00429         // Perform a line search along xdww = xdw + a * w
00430         
00431         if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS ) {
00432           out << "\nPerform linesearch along xdww = xdw + a*w\n"
00433             << "starting from a = " << a << " ...\n";
00434         }
00435 
00436         xdww = xdw();                 // xdww = xdw + a*w
00437         Vp_StV( &xdww(), a, w() );
00438         phi_c.calc_deriv(nlp.c());  // Set the directional derivative at c(xdw)
00439         phi_c_xdww = phi_c_calc( xdww() );  // phi_c_xdww = phi(xdww)
00440         const DVectorSlice xdw_w[2] = { xdw(), w() };
00441         MeritFuncCalc1DQuadratic
00442           phi_c_calc_1d( phi_c_calc, 1 , xdw_w, &xdww() );
00443         bool ls_okay = false;
00444         try {
00445           ls_okay = direct_ls_newton().do_line_search(
00446             phi_c_calc_1d,phi_c_xdw
00447             ,&a,&phi_c_xdww
00448             , (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS 
00449             ? &out : 0                        
00450             );
00451         }
00452         catch(const DirectLineSearch_Strategy::NotDescentDirection& excpt ) {
00453           if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO ) {
00454             out << "\nThe line search object throw the exception:" << typeName(excpt) << ":\n"
00455               << excpt.what() << std::endl;
00456           }
00457           ls_okay = false;
00458         }
00459         // Note that the last value c(x) computed by the line search is for
00460         // xdw + a*w.
00461 
00462         // Print line for summary output
00463         if( use_newton_olevel == PRINT_NEWTON_SUMMARY_INFO ) {
00464           out << setw(5)      << newton_i
00465             << setw(dbl_w)    << nrm_w
00466             << setw(dbl_w)    << u
00467             << setw(dbl_w)    << step_ratio
00468             << setw(5)      << direct_ls_newton().num_iterations()
00469             << setw(dbl_w)    << a
00470             << setw(dbl_w)    << phi_c_xdww
00471             << setw(dbl_w)    << (phi_c_xdww-phi_c_x)
00472             << setw(dbl_w)    << (phi_c_xdww-phi_c_xd) << std::endl;
00473         }
00474 
00475         if(!ls_okay) {
00476           if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO ) {
00477             out << "\nThe line search failed so forget about computing a correction ...\n";
00478           }
00479           use_correction = false;
00480           newton_failed = true;
00481           break;
00482         }
00483 
00484         // See if this point is okay
00485         bool good_correction = false;
00486         switch( forced_constr_reduction() ) {
00487           case CONSTR_LESS_X_D: {
00488             good_correction = ( phi_c_xdww < forced_reduct_ratio()*phi_c_xd );
00489             if( good_correction
00490               && (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO )
00491             {
00492               out << "\nphi_c(c(x_k+d_k+w)) = " << phi_c_xdww
00493                 << " < forced_reduct_ratio * phi_c(c(x_k+d_k)) = "
00494                 << forced_reduct_ratio() << " * " << phi_c_xd
00495                 << " = " << (forced_reduct_ratio()*phi_c_xd) << std::endl;
00496             }
00497             break;
00498           }
00499           case CONSTR_LESS_X: {
00500             good_correction = ( phi_c_xdww < forced_reduct_ratio()*phi_c_x );
00501             if( good_correction
00502               && (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO )
00503             {
00504               out << "\nphi_c(c(x_k+d_k+w)) = " << phi_c_xdww
00505                 << " < forced_reduct_ratio * phi_c(c(x_k)) = "
00506                 << forced_reduct_ratio() << " * " << phi_c_x
00507                 << " = " << (forced_reduct_ratio()*phi_c_x) << std::endl;
00508             }
00509             break;
00510           }
00511           default:
00512             TEST_FOR_EXCEPT(true);
00513         }
00514 
00515         if(good_correction) {
00516           if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO ) {
00517             out << "\nAccept this point and compute our full correction w ... \n";
00518           }
00519           // Compute the full correction and do a curved linesearch
00520           // w = xdww - x_kp1
00521           V_VmV( &w(), xdww(), x_kp1() );
00522           if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO ) {
00523             out << "\n||w||inf = " << norm_inf(w()) << std::endl;
00524           }
00525           if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_VECTORS ) {
00526             out << "\nw = " << w();
00527           }
00528           use_correction = true;
00529           newton_failed  = false;
00530           break;
00531         }
00532 
00533         // Else perform another newton iteration.
00534         xdw       = xdww;
00535         phi_c_xdw = phi_c_xdww;
00536 
00537       } // end for
00538       if( !use_correction ) {
00539         newton_failed  = true;
00540         if( forced_constr_reduction() == CONSTR_LESS_X_D ) {
00541           if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00542             out << "\nDam! This is really bad!\n"
00543               << "We where only looking for point phi_c(c(x_k+d_k+w)"
00544               << " < phi_c(c(x_k+k_k) and we could not find it\n"
00545               << " in the aloted number of newton iterations!\n"
00546               << "Perhaps the Gc_k did not give us a descent direction?\n"
00547               << "Just perform a standard line search from here ...\n";
00548           }
00549           use_correction = false;
00550         }
00551         else {
00552           if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00553             out << "\nWe where looking for point phi_c(c(x_k+d_k+w))"
00554               << " < phi_c(c(x_k)) and we could not find it.\n";
00555           }
00556           if( phi_c_xdww < phi_c_xd ) {
00557             if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00558               out << "However, we did find a point less than phi_c(c(x_k+d_k))\n"
00559                 << "so lets use the correction anyway.\n";
00560             }
00561             // Compute the full correction and do a curved linesearch
00562             // w = xdww - x_kp1
00563             V_VmV( &w(), xdww(), x_kp1() );
00564             if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO ) {
00565               out << "\n||w||inf = " << norm_inf(w()) << std::endl;
00566             }
00567             if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_VECTORS ) {
00568               out << "\nw = " << w();
00569             }
00570             use_correction = true;
00571           }
00572           else {
00573             if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00574               out << "Dam! We did not even find a point less than phi_c(c(x_k+d_k))\n"
00575                 << "just perform a standard line search along x_k + alpha_k*d_k.\n";
00576             }
00577             use_correction = false;
00578           }
00579         }
00580       }
00581     } // end else from if phi_c_xdw > phi_c_x
00582   } // end considered_correction
00583 
00584   // //////////////////////////
00585   // Set up for the line search
00586 
00587   if( considered_correction ) {
00588     if( use_correction ) {
00589       // We are using the correction so setup the full step for the
00590       // NLP linesearch to come.
00591       Vp_V( &x_kp1(), w() );      // Set point to x_kp1 = x_k + d_k + w
00592       nlp.calc_f(x_kp1(),false);  // same as last call to calc_c(x)
00593       f_kp1 = nlp.f();            // Here f and c where computed at x_k+d_k+w
00594       c_kp1 = nlp.c()();
00595       phi_kp1 = merit_func().value( f_kp1, c_kp1 );
00596     }
00597     else {
00598       // Just pretend the second order correction never happened
00599       // and we don't need to do anything.
00600     }
00601     // Set back the base point
00602     nlp.set_f( &f_kp1 );
00603     nlp.set_c( &c_kp1 );
00604   }
00605 
00606   // //////////////////////
00607   // Do the line search
00608 
00609   if( !chose_point ) {
00610     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00611       if( use_correction ) {
00612         out << "\nPerform a curved linesearch along:\n"
00613           << "x_kp1 = x_k + alpha_k * d_k + alpha_k^2 * w ...\n";
00614       }
00615       else {
00616         out << "\nPerform a standard linesearch along:\n"
00617           << "x_kp1 = x_k + alpha_k * d_k ...\n";
00618       }
00619     }
00620     const DVectorSlice xdw[3] = { x_k(), d_k(), w() };
00621     MeritFuncCalc1DQuadratic
00622       phi_calc_1d( phi_calc, (use_correction?2:1) , xdw, &x_kp1() );
00623     if( !direct_ls_sqp().do_line_search( phi_calc_1d, phi_k, &alpha_k, &phi_kp1
00624       , static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ?
00625         &out : static_cast<std::ostream*>(0)  )   )
00626     {
00627       // If the line search failed but the value of the merit function is less than
00628       // the base point value then just accept it and move on.  This many be a
00629       // bad thing to do.
00630 
00631       const value_type
00632         scaled_ared   = (s.phi().get_k(0) - s.phi().get_k(+1))/s.phi().get_k(0),
00633         keep_on_frac  = 1.0e-10;  // Make adjustable?
00634       bool keep_on = scaled_ared < keep_on_frac;
00635 
00636       if( (int)olevel >= (int)PRINT_BASIC_ALGORITHM_INFO )
00637       {
00638         out
00639           << "\nThe maximum number of linesearch iterations has been exceeded "
00640           << "(k = " << algo.state().k() << ")\n"
00641           << "(phi_k - phi_kp1)/phi_k = " << scaled_ared;
00642 //        if(keep_on) {
00643 //          out
00644 //            << " < " << keep_on_frac
00645 //            << "\nso we will accept to step and move on.\n";
00646 //        }
00647 //        else {
00648 //          out
00649 //            << " > " << keep_on_frac
00650 //            << "\nso we will reject the step and declare a line search failure.\n";
00651 //        }
00652       }
00653 //
00654 //      if( keep_on ) return true;
00655       
00656       throw LineSearchFailure( "LineSearch2ndOrderCorrect_Step::do_step(): "
00657                    "Error, Line search failure" );
00658     }
00659   }
00660 
00661   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00662     out << "\nalpha_k      = "  << alpha_k        << std::endl
00663       << "\n||x_kp1||inf = "  << norm_inf( x_kp1 )  << std::endl
00664       << "\nf_kp1        = "  << f_kp1        << std::endl
00665       << "\n||c_kp1||inf = "  << norm_inf(c_kp1)    << std::endl
00666       << "\nphi_kp1      = "  << phi_kp1        << std::endl;
00667   }
00668 
00669   if( (int)olevel >= (int)PRINT_VECTORS ) {
00670     out << "\nx_kp1 =\n"  << x_kp1
00671       << "\nc_kp1 =\n"  << c_kp1;
00672   }
00673 
00674 */
00675   TEST_FOR_EXCEPT(true);
00676 
00677   return true;
00678 }
00679 
00680 void LineSearch2ndOrderCorrect_Step::print_step(
00681   const Algorithm& algo
00682   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00683   , std::ostream& out, const std::string& L ) const
00684 {
00685   out << L << "*** Calculate a second order correction when near solution.\n"
00686     << L << "*** If we can compute a correction w then perform a curved\n"
00687     << L << "*** line search along x_kp1 = x_k + alpha_k * d_k + alpha_k^2 * w.\n"
00688     << L << "default: eta                     = " << eta() << std::endl
00689     << L << "         constr_norm_threshold   = " << constr_norm_threshold() << std::endl
00690     << L << "         constr_incr_ratio       = " << constr_incr_ratio() << std::endl
00691     << L << "         after_k_iter            = " << after_k_iter() << std::endl
00692     << L << "         forced_constr_reduction = " << (forced_constr_reduction()== CONSTR_LESS_X_D
00693                               ? "CONSTR_LESS_X_D\n"
00694                               : "CONSTR_LESS_X\n" )
00695     << L << "         forced_reduct_ratio     = " << forced_reduct_ratio() << std::endl
00696     << L << "         max_step_ratio          = " << max_step_ratio() << std::endl
00697     << L << "         max_newton_iter         = " << max_newton_iter() << std::endl
00698     << L << "begin definition of NLP merit function phi.value(f(x),c(x)):\n";
00699   
00700   merit_func().print_merit_func( out, L + "  " );
00701   
00702   out << L << "end definition\n"
00703     << L << "Dphi_k = phi.deriv()\n"
00704     << L << "if Dphi_k >= 0 then\n"
00705     << L << "  throw line_search_failure\n"
00706     << L << "end\n"
00707     << L << "phi_kp1 = phi_k.value(f_kp1,c_kp1)\n"
00708     << L << "phi_k = phi.value(f_k,c_k)\n"
00709     << L << "if (norm(c_k,inf) <= constr_norm_threshold)\n"
00710     << L << "and (norm(c_kp1,inf)/(norm(c_k,inf)+small_num) <= constr_incr_ratio)\n"
00711     << L << "and (k >= after_k_iter) then\n"
00712     << L << "considering_correction = ( (norm(c_k,inf) <= constr_norm_threshold)\n"
00713     << L << "  and (norm(c_kp1,inf)/(1.0 + norm(c_k,inf)) <= constr_incr_ratio)\n"
00714     << L << "  and (k >= after_k_iter) )\n"
00715     << L << "chose_point = false\n"
00716     << L << "if phi_kp1 < phi_k + eta * Dphi_k then\n"
00717     << L << "  chose_point = true\n"
00718     << L << "else\n"
00719     << L << "if (considering_correction == true) and (chose_point == false) then\n"
00720     << L << "  considered_correction = true\n"
00721     << L << "  begin definition of c(x) merit function phi_c.value(c(x)):\n";
00722 
00723   ConstrainedOptPack::MeritFuncNLESqrResid().print_merit_func(
00724     out, L + "    " );
00725 
00726   out << L << "  end definition\n"
00727     << L << "  xdw = x_kp1;\n"
00728     << L << "  phi_c_x = phi_c.value(c_k);\n"
00729     << L << "  phi_c_xd = phi_c.value(c_kp1);\n"
00730     << L << "  phi_c_xdw = phi_c_xd;\n"
00731     << L << "  phi_c_xdww = phi_c_xdw;\n"
00732     << L << "  if phi_c_xd < forced_reduct_ratio * phi_c_x then\n"
00733     << L << "    *** There is no need to perform a correction.\n"
00734     << L << "    use_correction = false;\n"
00735     << L << "  else\n"
00736     << L << "    *** Lets try to compute a correction by performing\n"
00737     << L << "    *** a series of newton steps to compute local steps w\n"
00738     << L << "    for newton_i = 1...max_newton_itr\n"
00739     << L << "      begin feasibility step calculation: \"" << typeName(feasibility_step()) << "\"\n";
00740 
00741   feasibility_step().print_step( out, L + "        " );
00742 
00743   out << L << "      end feasibility step calculation\n"
00744     << L << "      Find the largest positive step u where the slightly\n"
00745     << L << "      relaxed variable bounds:\n"
00746     << L << "        xl - delta <= xdw + u * w <= xu + delta\n"
00747     << L << "      are strictly satisfied (where delta = max_var_bounds_viol).\n"
00748     << L << "      a = min(1,u);\n"
00749     << L << "      step_ratio = norm(w,inf)/norm(d,inf);\n"
00750     << L << "      a = min(a,max_step_ratio/step_ratio);\n"
00751     << L << "      Perform line search for phi_c.value(c(xdww = xdw+a*w))\n"
00752     << L << "      starting from a and compute:\n"
00753     << L << "        a,\n"
00754     << L << "        xdww = xdw + a * w,\n"
00755     << L << "        phi_c_xdww = phi_c.value(c(xdww))\n"
00756     << L << "      print summary information;\n"
00757     << L << "      if line search failed then\n"
00758     << L << "        use_correction = false;\n"
00759     << L << "        exit for loop;\n"
00760     << L << "      end\n"
00761     << L << "      *** Determine if this is sufficent reduction in c(x) error\n"
00762     << L << "      if forced_constr_reduction == CONSTR_LESS_X_D then\n"
00763     << L << "        good_correction = (phi_c.value(c(xdww)) < forced_reduct_ratio*phi_c_xd);\n"
00764     << L << "      else if forced_constr_reduction == CONSTR_LESS_X then\n"
00765     << L << "        good_correction = (phi_c.value(c(xdww)) < forced_reduct_ratio*phi_c_x);\n"
00766     << L << "      end\n"
00767     << L << "      if good_correction == true then\n"
00768     << L << "        w = xdww - (x_k+d_k);\n"
00769     << L << "        use_correction = true;\n"
00770     << L << "        exit for loop;\n"
00771     << L << "      end\n"
00772     << L << "      *** This is not sufficient reduction is c(x) error yet.\n"
00773     << L << "      xdw = xdww;\n"
00774     << L << "      phi_c_xdw = phi_c_xdww;\n"
00775     << L << "    end\n"
00776     << L << "    if use_correction == false then\n"
00777     << L << "      if forced_constr_reduction == CONSTR_LESS_X_D then\n"
00778     << L << "        *** Dam! We could not find a point phi_c_xdww < phi_c_xd.\n"
00779     << L << "        *** Perhaps Gc_k does not give a descent direction for phi_c!\n"
00780     << L << "      else if forced_constr_reduction == CONSTR_LESS_X then\n"
00781     << L << "        if phi_c_dww < phi_c_xd then\n"
00782     << L << "           *** Accept this correction anyway.\n"
00783     << L << "           use_correction = true\n"
00784     << L << "        else\n"
00785     << L << "          *** Dam! we could not find any reduction in phi_c so\n"
00786     << L << "          *** Perhaps Gc_k does not give a descent direction for phi_c!\n"
00787     << L << "      end\n"
00788     << L << "    end\n"
00789     << L << "  end\n"
00790     << L << "end\n"
00791     << L << "if chose_point == false then\n"
00792     << L << "  if use_correction == true then\n"
00793     << L << "    Perform line search along x_kp1 = x_k + alpha_k * d_k + alpha_k^2 * w\n"
00794     << L << "  else\n"
00795     << L << "    Perform line search along x_kp1 = x_k + alpha_k * d_k\n"
00796     << L << "  end\n"
00797     << L << "  begin direct line search : \"" << typeName(direct_ls_sqp()) << "\"\n";
00798 
00799   direct_ls_sqp().print_algorithm( out, L + "    " );
00800 
00801   out
00802     << L << "  end direct line search\n"
00803     << L << "  if maximum number of linesearch iterations are exceeded then\n"
00804     << L << "    throw line_search_failure\n"
00805     << L << "  end\n"
00806     << L << "end\n";
00807 }
00808 
00809 } // end namespace MoochoPack

Generated on Tue Oct 20 12:51:47 2009 for MOOCHO (Single Doxygen Collection) by doxygen 1.4.7