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