MOOCHO (Single Doxygen Collection) Version of the Day
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <ostream>
00043 #include <typeinfo>
00044 
00045 #include "MoochoPack_LineSearchDirect_Step.hpp"
00046 #include "MoochoPack_Exceptions.hpp"
00047 #include "MoochoPack_moocho_algo_conversion.hpp"
00048 #include "IterationPack_print_algorithm_step.hpp"
00049 #include "ConstrainedOptPack_MeritFuncCalc1DQuadratic.hpp"
00050 #include "ConstrainedOptPack_MeritFuncCalcNLP.hpp"
00051 #include "AbstractLinAlgPack_VectorMutable.hpp"
00052 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00053 #include "AbstractLinAlgPack_VectorOut.hpp"
00054 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00055 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00056 #include "Teuchos_Assert.hpp"
00057 
00058 namespace MoochoPack {
00059 
00060 LineSearchDirect_Step::LineSearchDirect_Step(
00061   const direct_line_search_ptr_t& direct_line_search
00062   )
00063   :direct_line_search_(direct_line_search)
00064 {}
00065 
00066 bool LineSearchDirect_Step::do_step(
00067   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00068   ,poss_type assoc_step_poss
00069   )
00070 {
00071   using AbstractLinAlgPack::Vp_StV;
00072   using LinAlgOpPack::V_VpV;
00073 
00074   NLPAlgo &algo = rsqp_algo(_algo);
00075   NLPAlgoState  &s = algo.rsqp_state();
00076   NLP &nlp = algo.nlp();
00077 
00078   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00079   std::ostream& out = algo.track().journal_out();
00080 
00081   // print step header.
00082   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00083     using IterationPack::print_algorithm_step;
00084     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00085   }
00086 
00087   Teuchos::VerboseObjectTempState<NLP>
00088     nlpOutputTempState(
00089       Teuchos::rcp(&nlp,false), Teuchos::getFancyOStream(Teuchos::rcp(&out,false)),
00090       Teuchos::VERB_DEFAULT );
00091   // Above, we don't want to litter the output with any trace from the NLP.
00092   // However, if the user forces the verbosity level to be increased, then we
00093   // want to set the stream so that it knows where to print to.
00094 
00095   const size_type
00096     m  = nlp.m();
00097 
00098   // /////////////////////////////////////////
00099   // Set references to iteration quantities
00100   //
00101   // Set k+1 first then go back to get k+0 to ensure
00102   // we have backward storage!
00103 
00104   IterQuantityAccess<value_type>
00105     &f_iq     = s.f(),
00106     &alpha_iq = s.alpha(),
00107     &phi_iq   = s.phi();
00108   IterQuantityAccess<VectorMutable>
00109     &x_iq    = s.x(),
00110     &d_iq    = s.d(),
00111     &c_iq    = s.c();
00112   
00113   VectorMutable        &x_kp1   = x_iq.get_k(+1);
00114   const Vector         &x_k     = x_iq.get_k(0);
00115   value_type           &f_kp1   = f_iq.get_k(+1);
00116   const value_type     &f_k     = f_iq.get_k(0);
00117   VectorMutable        *c_kp1   = m  ? &c_iq.get_k(+1) : NULL;
00118   const Vector         *c_k     = m  ? &c_iq.get_k(0)  : NULL;
00119   const Vector         &d_k     = d_iq.get_k(0);
00120   value_type           &alpha_k = alpha_iq.get_k(0);
00121 
00122   // /////////////////////////////////////
00123   // Compute Dphi_k, phi_kp1 and phi_k
00124 
00125   const MeritFuncNLP
00126     &merit_func_nlp_k = s.merit_func_nlp().get_k(0);
00127   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00128     out << "\nBegin definition of NLP merit function phi.value(f(x),c(x)):\n";
00129     merit_func_nlp_k.print_merit_func( out, "    " );
00130     out << "end definition of the NLP merit funciton\n";
00131   }
00132   // Dphi_k
00133   const value_type
00134     Dphi_k = merit_func_nlp_k.deriv();
00135   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00136     out << "\nDphi_k = "  << Dphi_k << std::endl;
00137   }
00138   TEUCHOS_TEST_FOR_EXCEPTION(
00139     Dphi_k >= 0, LineSearchFailure
00140     ,"LineSearch2ndOrderCorrect_Step::do_step(...) : " 
00141     "Error, d_k is not a descent direction for the merit function "
00142     "since Dphi_k = " << Dphi_k << " >= 0" );
00143   // ph_kp1
00144   value_type
00145     &phi_kp1 = phi_iq.set_k(+1) = merit_func_nlp_k.value(
00146       f_kp1
00147       ,c_kp1
00148       ,NULL     // h
00149       ,NULL     // hl
00150       ,NULL     // hu
00151       );
00152   // Must compute phi(x) at the base point x_k since the penalty parameter may have changed.
00153   const value_type
00154     &phi_k = phi_iq.set_k(0) = merit_func_nlp_k.value(
00155       f_k
00156       ,c_k
00157       ,NULL     // h
00158       ,NULL     // hl
00159       ,NULL     // hu
00160       );
00161   
00162   // //////////////////////////////////////
00163   // Setup the calculation merit function
00164 
00165   // Here f_kp1, c_kp1 and h_kp1 are updated at the same time the
00166   // line search is being performed!
00167   nlp.unset_quantities();
00168   nlp.set_f( &f_kp1 );
00169   if(m)  nlp.set_c( c_kp1 );
00170   MeritFuncCalcNLP
00171     phi_calc( &merit_func_nlp_k, &nlp );
00172 
00173   // //////////////////////
00174   // Do the line search
00175   
00176   const Vector* xd[2] = { &x_k, &d_k };
00177   MeritFuncCalc1DQuadratic
00178     phi_calc_1d( phi_calc, 1, xd, &x_kp1 );
00179   
00180   if( !direct_line_search().do_line_search(
00181       phi_calc_1d, phi_k, &alpha_k, &phi_kp1
00182       ,( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS)
00183          ? &out : static_cast<std::ostream*>(0) )
00184       )
00185     )
00186   {
00187     // The line search failed!
00188     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) )
00189       out
00190         << "\nThe maximum number of linesearch iterations has been exceeded "
00191         << "(k = " << algo.state().k() << ")\n"
00192         << "(phi_k - phi_kp1)/phi_k = " << ((phi_k - phi_kp1)/phi_k)
00193         << "\nso we will reject the step and declare a line search failure.\n";
00194     TEUCHOS_TEST_FOR_EXCEPTION(
00195       true, LineSearchFailure
00196       ,"LineSearchDirect_Step::do_step(): Line search failure" );
00197   }
00198   nlp.unset_quantities();
00199   
00200   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00201     out << "\nalpha_k                = " << alpha_k;
00202     out << "\n||x_kp1||inf           = " << x_kp1.norm_inf();
00203     out << "\nf_kp1                  = " << f_kp1;
00204     if(m)
00205       out << "\n||c_kp1||inf           = " << c_kp1->norm_inf();
00206     out << "\nphi_kp1                = " << phi_kp1;
00207     out << std::endl;
00208   }
00209   
00210   if( (int)olevel >= (int)PRINT_VECTORS ) {
00211     out << "\nx_kp1 =\n" << x_kp1;
00212     if(m)
00213       out <<"\nc_kp1 =\n" << *c_kp1;
00214   }
00215 
00216   return true;
00217 }
00218 
00219 void LineSearchDirect_Step::print_step(
00220   const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00221   ,std::ostream& out, const std::string& L
00222   ) const
00223 {
00224   out
00225     << L << "*** Preform a line search along the full space search direction d_k.\n"
00226     << L << "Dphi_k = merit_func_nlp_k.deriv()\n"
00227     << L << "if Dphi_k >= 0 then\n"
00228     << L << "    throw line_search_failure\n"
00229     << L << "end\n"
00230     << L << "phi_kp1 = merit_func_nlp_k.value(f_kp1,c_kp1,h_kp1,hl,hu)\n"
00231     << L << "phi_k = merit_func_nlp_k.value(f_k,c_k,h_k,hl,hu)\n"
00232     << L << "begin direct line search (where phi = merit_func_nlp_k): \"" << typeName(direct_line_search()) << "\"\n";
00233   direct_line_search().print_algorithm( out, L + "    " );
00234   out
00235     << L << "end direct line search\n"
00236     << L << "if maximum number of linesearch iterations are exceeded then\n"
00237     << L << "    throw line_search_failure\n"
00238     << L << "end\n";
00239 }
00240 
00241 } // end namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines