MoochoPack_LineSearchDirect_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_LineSearchDirect_Step.hpp"
00033 #include "MoochoPack_Exceptions.hpp"
00034 #include "MoochoPack_moocho_algo_conversion.hpp"
00035 #include "IterationPack_print_algorithm_step.hpp"
00036 #include "ConstrainedOptPack_MeritFuncCalc1DQuadratic.hpp"
00037 #include "ConstrainedOptPack_MeritFuncCalcNLP.hpp"
00038 #include "AbstractLinAlgPack_VectorMutable.hpp"
00039 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00040 #include "AbstractLinAlgPack_VectorOut.hpp"
00041 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00042 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00043 #include "Teuchos_TestForException.hpp"
00044 
00045 namespace MoochoPack {
00046 
00047 LineSearchDirect_Step::LineSearchDirect_Step(
00048   const direct_line_search_ptr_t& direct_line_search
00049   )
00050   :direct_line_search_(direct_line_search)
00051 {}
00052 
00053 bool LineSearchDirect_Step::do_step(
00054   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00055   ,poss_type assoc_step_poss
00056   )
00057 {
00058   using AbstractLinAlgPack::Vp_StV;
00059   using LinAlgOpPack::V_VpV;
00060 
00061   NLPAlgo &algo = rsqp_algo(_algo);
00062   NLPAlgoState  &s    = algo.rsqp_state();
00063   NLP     &nlp  = algo.nlp();
00064 
00065   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00066   std::ostream& out = algo.track().journal_out();
00067 
00068   // print step header.
00069   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00070     using IterationPack::print_algorithm_step;
00071     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00072   }
00073 
00074   const size_type
00075     m  = nlp.m();
00076 
00077   // /////////////////////////////////////////
00078   // Set references to iteration quantities
00079   //
00080   // Set k+1 first then go back to get k+0 to ensure
00081   // we have backward storage!
00082 
00083   IterQuantityAccess<value_type>
00084     &f_iq     = s.f(),
00085     &alpha_iq = s.alpha(),
00086     &phi_iq   = s.phi();
00087   IterQuantityAccess<VectorMutable>
00088     &x_iq    = s.x(),
00089     &d_iq    = s.d(),
00090     &c_iq    = s.c();
00091   
00092   VectorMutable        &x_kp1   = x_iq.get_k(+1);
00093   const Vector         &x_k     = x_iq.get_k(0);
00094   value_type           &f_kp1   = f_iq.get_k(+1);
00095   const value_type     &f_k     = f_iq.get_k(0);
00096   VectorMutable        *c_kp1   = m  ? &c_iq.get_k(+1) : NULL;
00097   const Vector         *c_k     = m  ? &c_iq.get_k(0)  : NULL;
00098   const Vector         &d_k     = d_iq.get_k(0);
00099   value_type           &alpha_k = alpha_iq.get_k(0);
00100 
00101   // /////////////////////////////////////
00102   // Compute Dphi_k, phi_kp1 and phi_k
00103 
00104   const MeritFuncNLP
00105     &merit_func_nlp_k = s.merit_func_nlp().get_k(0);
00106   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00107     out << "\nBegin definition of NLP merit function phi.value(f(x),c(x)):\n";
00108     merit_func_nlp_k.print_merit_func( out, "    " );
00109     out << "end definition of the NLP merit funciton\n";
00110   }
00111   // Dphi_k
00112   const value_type
00113     Dphi_k = merit_func_nlp_k.deriv();
00114   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00115     out << "\nDphi_k = "  << Dphi_k << std::endl;
00116   }
00117   TEST_FOR_EXCEPTION(
00118     Dphi_k >= 0, LineSearchFailure
00119     ,"LineSearch2ndOrderCorrect_Step::do_step(...) : " 
00120     "Error, d_k is not a descent direction for the merit function "
00121     "since Dphi_k = " << Dphi_k << " >= 0" );
00122   // ph_kp1
00123   value_type
00124     &phi_kp1 = phi_iq.set_k(+1) = merit_func_nlp_k.value(
00125       f_kp1
00126       ,c_kp1
00127       ,NULL     // h
00128       ,NULL     // hl
00129       ,NULL     // hu
00130       );
00131   // Must compute phi(x) at the base point x_k since the penalty parameter may have changed.
00132   const value_type
00133     &phi_k = phi_iq.set_k(0) = merit_func_nlp_k.value(
00134       f_k
00135       ,c_k
00136       ,NULL     // h
00137       ,NULL     // hl
00138       ,NULL     // hu
00139       );
00140   
00141   // //////////////////////////////////////
00142   // Setup the calculation merit function
00143 
00144   // Here f_kp1, c_kp1 and h_kp1 are updated at the same time the
00145   // line search is being performed!
00146   nlp.unset_quantities();
00147   nlp.set_f( &f_kp1 );
00148   if(m)  nlp.set_c( c_kp1 );
00149   MeritFuncCalcNLP
00150     phi_calc( &merit_func_nlp_k, &nlp );
00151 
00152   // //////////////////////
00153   // Do the line search
00154   
00155   const Vector* xd[2] = { &x_k, &d_k };
00156   MeritFuncCalc1DQuadratic
00157     phi_calc_1d( phi_calc, 1, xd, &x_kp1 );
00158   
00159   if( !direct_line_search().do_line_search(
00160       phi_calc_1d, phi_k, &alpha_k, &phi_kp1
00161       ,( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS)
00162          ? &out : static_cast<std::ostream*>(0) )
00163       )
00164     )
00165   {
00166     // The line search failed!
00167     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) )
00168       out
00169         << "\nThe maximum number of linesearch iterations has been exceeded "
00170         << "(k = " << algo.state().k() << ")\n"
00171         << "(phi_k - phi_kp1)/phi_k = " << ((phi_k - phi_kp1)/phi_k)
00172         << "\nso we will reject the step and declare a line search failure.\n";
00173     TEST_FOR_EXCEPTION(
00174       true, LineSearchFailure
00175       ,"LineSearchDirect_Step::do_step(): Line search failure" );
00176   }
00177   nlp.unset_quantities();
00178   
00179   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00180     out << "\nalpha_k                = " << alpha_k;
00181     out << "\n||x_kp1||inf           = " << x_kp1.norm_inf();
00182     out << "\nf_kp1                  = " << f_kp1;
00183     if(m)
00184       out << "\n||c_kp1||inf           = " << c_kp1->norm_inf();
00185     out << "\nphi_kp1                = " << phi_kp1;
00186     out << std::endl;
00187   }
00188   
00189   if( (int)olevel >= (int)PRINT_VECTORS ) {
00190     out << "\nx_kp1 =\n" << x_kp1;
00191     if(m)
00192       out <<"\nc_kp1 =\n" << *c_kp1;
00193   }
00194 
00195   return true;
00196 }
00197 
00198 void LineSearchDirect_Step::print_step(
00199   const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00200   ,std::ostream& out, const std::string& L
00201   ) const
00202 {
00203   out
00204     << L << "*** Preform a line search along the full space search direction d_k.\n"
00205     << L << "Dphi_k = merit_func_nlp_k.deriv()\n"
00206     << L << "if Dphi_k >= 0 then\n"
00207     << L << "    throw line_search_failure\n"
00208     << L << "end\n"
00209     << L << "phi_kp1 = merit_func_nlp_k.value(f_kp1,c_kp1,h_kp1,hl,hu)\n"
00210     << L << "phi_k = merit_func_nlp_k.value(f_k,c_k,h_k,hl,hu)\n"
00211     << L << "begin direct line search (where phi = merit_func_nlp_k): \"" << typeid(direct_line_search()).name() << "\"\n";
00212   direct_line_search().print_algorithm( out, L + "    " );
00213   out
00214     << L << "end direct line search\n"
00215     << L << "if maximum number of linesearch iterations are exceeded then\n"
00216     << L << "    throw line_search_failure\n"
00217     << L << "end\n";
00218 }
00219 
00220 } // end namespace MoochoPack

Generated on Thu Sep 18 12:35:17 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1