MoochoPack_LineSearch2ndOrderCorrect_Step.cpp

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

Generated on Wed May 12 21:51:20 2010 for MoochoPack : Framework for Large-Scale Optimization Algorithms by  doxygen 1.4.7