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   Teuchos::VerboseObjectTempState<NLP>
00075     nlpOutputTempState(
00076       Teuchos::rcp(&nlp,false), Teuchos::getFancyOStream(Teuchos::rcp(&out,false)),
00077       Teuchos::VERB_DEFAULT );
00078   // Above, we don't want to litter the output with any trace from the NLP.
00079   // However, if the user forces the verbosity level to be increased, then we
00080   // want to set the stream so that it knows where to print to.
00081 
00082   const size_type
00083     m  = nlp.m();
00084 
00085   // /////////////////////////////////////////
00086   // Set references to iteration quantities
00087   //
00088   // Set k+1 first then go back to get k+0 to ensure
00089   // we have backward storage!
00090 
00091   IterQuantityAccess<value_type>
00092     &f_iq     = s.f(),
00093     &alpha_iq = s.alpha(),
00094     &phi_iq   = s.phi();
00095   IterQuantityAccess<VectorMutable>
00096     &x_iq    = s.x(),
00097     &d_iq    = s.d(),
00098     &c_iq    = s.c();
00099   
00100   VectorMutable        &x_kp1   = x_iq.get_k(+1);
00101   const Vector         &x_k     = x_iq.get_k(0);
00102   value_type           &f_kp1   = f_iq.get_k(+1);
00103   const value_type     &f_k     = f_iq.get_k(0);
00104   VectorMutable        *c_kp1   = m  ? &c_iq.get_k(+1) : NULL;
00105   const Vector         *c_k     = m  ? &c_iq.get_k(0)  : NULL;
00106   const Vector         &d_k     = d_iq.get_k(0);
00107   value_type           &alpha_k = alpha_iq.get_k(0);
00108 
00109   // /////////////////////////////////////
00110   // Compute Dphi_k, phi_kp1 and phi_k
00111 
00112   const MeritFuncNLP
00113     &merit_func_nlp_k = s.merit_func_nlp().get_k(0);
00114   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00115     out << "\nBegin definition of NLP merit function phi.value(f(x),c(x)):\n";
00116     merit_func_nlp_k.print_merit_func( out, "    " );
00117     out << "end definition of the NLP merit funciton\n";
00118   }
00119   // Dphi_k
00120   const value_type
00121     Dphi_k = merit_func_nlp_k.deriv();
00122   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00123     out << "\nDphi_k = "  << Dphi_k << std::endl;
00124   }
00125   TEST_FOR_EXCEPTION(
00126     Dphi_k >= 0, LineSearchFailure
00127     ,"LineSearch2ndOrderCorrect_Step::do_step(...) : " 
00128     "Error, d_k is not a descent direction for the merit function "
00129     "since Dphi_k = " << Dphi_k << " >= 0" );
00130   // ph_kp1
00131   value_type
00132     &phi_kp1 = phi_iq.set_k(+1) = merit_func_nlp_k.value(
00133       f_kp1
00134       ,c_kp1
00135       ,NULL     // h
00136       ,NULL     // hl
00137       ,NULL     // hu
00138       );
00139   // Must compute phi(x) at the base point x_k since the penalty parameter may have changed.
00140   const value_type
00141     &phi_k = phi_iq.set_k(0) = merit_func_nlp_k.value(
00142       f_k
00143       ,c_k
00144       ,NULL     // h
00145       ,NULL     // hl
00146       ,NULL     // hu
00147       );
00148   
00149   // //////////////////////////////////////
00150   // Setup the calculation merit function
00151 
00152   // Here f_kp1, c_kp1 and h_kp1 are updated at the same time the
00153   // line search is being performed!
00154   nlp.unset_quantities();
00155   nlp.set_f( &f_kp1 );
00156   if(m)  nlp.set_c( c_kp1 );
00157   MeritFuncCalcNLP
00158     phi_calc( &merit_func_nlp_k, &nlp );
00159 
00160   // //////////////////////
00161   // Do the line search
00162   
00163   const Vector* xd[2] = { &x_k, &d_k };
00164   MeritFuncCalc1DQuadratic
00165     phi_calc_1d( phi_calc, 1, xd, &x_kp1 );
00166   
00167   if( !direct_line_search().do_line_search(
00168       phi_calc_1d, phi_k, &alpha_k, &phi_kp1
00169       ,( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS)
00170          ? &out : static_cast<std::ostream*>(0) )
00171       )
00172     )
00173   {
00174     // The line search failed!
00175     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) )
00176       out
00177         << "\nThe maximum number of linesearch iterations has been exceeded "
00178         << "(k = " << algo.state().k() << ")\n"
00179         << "(phi_k - phi_kp1)/phi_k = " << ((phi_k - phi_kp1)/phi_k)
00180         << "\nso we will reject the step and declare a line search failure.\n";
00181     TEST_FOR_EXCEPTION(
00182       true, LineSearchFailure
00183       ,"LineSearchDirect_Step::do_step(): Line search failure" );
00184   }
00185   nlp.unset_quantities();
00186   
00187   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00188     out << "\nalpha_k                = " << alpha_k;
00189     out << "\n||x_kp1||inf           = " << x_kp1.norm_inf();
00190     out << "\nf_kp1                  = " << f_kp1;
00191     if(m)
00192       out << "\n||c_kp1||inf           = " << c_kp1->norm_inf();
00193     out << "\nphi_kp1                = " << phi_kp1;
00194     out << std::endl;
00195   }
00196   
00197   if( (int)olevel >= (int)PRINT_VECTORS ) {
00198     out << "\nx_kp1 =\n" << x_kp1;
00199     if(m)
00200       out <<"\nc_kp1 =\n" << *c_kp1;
00201   }
00202 
00203   return true;
00204 }
00205 
00206 void LineSearchDirect_Step::print_step(
00207   const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00208   ,std::ostream& out, const std::string& L
00209   ) const
00210 {
00211   out
00212     << L << "*** Preform a line search along the full space search direction d_k.\n"
00213     << L << "Dphi_k = merit_func_nlp_k.deriv()\n"
00214     << L << "if Dphi_k >= 0 then\n"
00215     << L << "    throw line_search_failure\n"
00216     << L << "end\n"
00217     << L << "phi_kp1 = merit_func_nlp_k.value(f_kp1,c_kp1,h_kp1,hl,hu)\n"
00218     << L << "phi_k = merit_func_nlp_k.value(f_k,c_k,h_k,hl,hu)\n"
00219     << L << "begin direct line search (where phi = merit_func_nlp_k): \"" << typeName(direct_line_search()) << "\"\n";
00220   direct_line_search().print_algorithm( out, L + "    " );
00221   out
00222     << L << "end direct line search\n"
00223     << L << "if maximum number of linesearch iterations are exceeded then\n"
00224     << L << "    throw line_search_failure\n"
00225     << L << "end\n";
00226 }
00227 
00228 } // end namespace MoochoPack

Generated on Wed May 12 21:52:30 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7