MoochoPack : Framework for Large-Scale Optimization Algorithms Version of the Day
MoochoPack_LineSearchFilter_Step.cpp
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 <math.h>
00043 
00044 #include <ostream>
00045 #include <fstream>
00046 #include <typeinfo>
00047 
00048 #include "MoochoPack_LineSearchFilter_Step.hpp"
00049 #include "MoochoPack_Exceptions.hpp"
00050 #include "MoochoPack_moocho_algo_conversion.hpp"
00051 #include "NLPInterfacePack_NLPObjGrad.hpp"
00052 #include "IterationPack_print_algorithm_step.hpp"
00053 #include "AbstractLinAlgPack_VectorOut.hpp"
00054 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00055 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00056 #include "AbstractLinAlgPack_VectorMutableSubView.hpp"
00057 #include "Teuchos_dyn_cast.hpp"
00058 #include "Teuchos_Assert.hpp"
00059 #include "Teuchos_as.hpp"
00060 
00061 //#define FILTER_DEBUG_OUT 1
00062 
00063 
00064 namespace {
00065 
00066 
00067 class RemoveFilterEntryPred {
00068 public:
00069 
00070   RemoveFilterEntryPred(
00071     const MoochoPack::value_type &f_with_boundary,
00072     const MoochoPack::value_type &theta_with_boundary,
00073     const Teuchos::RCP<std::ostream> &out
00074     )
00075     :f_with_boundary_(f_with_boundary),
00076      theta_with_boundary_(theta_with_boundary),
00077      out_(out)
00078     {}
00079 
00080   bool operator()(const MoochoPack::FilterEntry &entry)
00081     {
00082       if (
00083         entry.f >= f_with_boundary_
00084         &&
00085         entry.theta >= theta_with_boundary_
00086         )
00087       {
00088         if (!is_null(out_)) {
00089           *out_
00090             << "\nRemoving from the filter the redundant point:"
00091             << "\n  f_with_boundary     = " << entry.f
00092             << "\n  theta_with_boundary = " << entry.theta
00093             << "\n  iteration added     = " << entry.iter
00094             << std::endl;
00095         }
00096         return true;
00097       }
00098       return false;
00099     }
00100 
00101 private:
00102 
00103   const MoochoPack::value_type f_with_boundary_;
00104   const MoochoPack::value_type theta_with_boundary_;
00105   const Teuchos::RCP<std::ostream> out_;
00106 
00107 };
00108 
00109 
00110 } // namespace
00111 
00112 
00113 namespace MoochoPack {
00114 
00115 // This must exist somewhere already, ask Ross
00116 value_type MIN(value_type x, value_type y)
00117 { return (x < y) ? x : y; }
00118 
00119 value_type MAX(value_type x, value_type y)
00120 { return (x > y) ? x : y; }
00121 
00122 LineSearchFilter_Step::LineSearchFilter_Step( 
00123   Teuchos::RCP<NLPInterfacePack::NLP> nlp
00124   ,const std::string obj_iq_name
00125   ,const std::string grad_obj_iq_name
00126   ,const value_type &gamma_theta
00127   ,const value_type &gamma_f
00128   ,const value_type &f_min
00129   ,const value_type &gamma_alpha
00130   ,const value_type &delta        
00131   ,const value_type &s_f
00132   ,const value_type &s_theta
00133   ,const value_type &theta_small_fact
00134   ,const value_type &theta_max
00135   ,const value_type &eta_f
00136   ,const value_type &back_track_frac
00137   )
00138   :
00139   gamma_theta_(gamma_theta),
00140   gamma_f_(gamma_f),
00141   f_min_(f_min),
00142   gamma_alpha_(gamma_alpha),
00143   delta_(delta),
00144   s_f_(s_f),
00145   s_theta_(s_theta),
00146   theta_small_fact_(theta_small_fact),
00147   theta_max_(theta_max),
00148   eta_f_(eta_f),
00149   back_track_frac_(back_track_frac),
00150   filter_(FILTER_IQ_STRING),
00151   obj_f_(obj_iq_name),
00152   grad_obj_f_(grad_obj_iq_name),
00153   nlp_(nlp)
00154 {
00155   TEUCHOS_TEST_FOR_EXCEPTION(
00156     !nlp_.get(),
00157     std::logic_error,
00158     "Null nlp passed to LineSearchFilter_Step constructor"
00159     );
00160 
00161 #if defined(FILTER_DEBUG_OUT)
00162   std::ofstream fout("filter_out.xml", std::ofstream::out | std::ofstream::trunc);
00163   fout << "<FilterDebugDocument>" << std::endl;
00164   fout.close();
00165 #endif
00166 }
00167 
00168 LineSearchFilter_Step::~LineSearchFilter_Step()
00169 {
00170 #if defined(FILTER_DEBUG_OUT)
00171   std::ofstream fout("filter_out.xml", std::ofstream::out | std::ofstream::app);
00172   fout << "</FilterDebugDocument>" << std::endl;
00173   fout.close();
00174 #endif
00175 }
00176   
00177 bool LineSearchFilter_Step::do_step(
00178   Algorithm& _algo, poss_type step_poss, 
00179   IterationPack::EDoStepType type
00180   ,poss_type assoc_step_poss)
00181 {
00182 
00183   using Teuchos::dyn_cast;
00184   using IterationPack::print_algorithm_step;
00185   using LinAlgOpPack::Vp_StV;
00186   using std::setw;
00187     
00188   // Get Algorithm (cast), state, and problem
00189   NLPAlgo &algo = rsqp_algo(_algo);
00190   NLPAlgoState &s = algo.rsqp_state();
00191 
00192   EJournalOutputLevel olevel  = algo.algo_cntr().journal_output_level();
00193   std::ostream &out = algo.track().journal_out();
00194     
00195   // print step header
00196   if (static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS)) 
00197   { 
00198     using IterationPack::print_algorithm_step;
00199     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00200   }
00201 
00202   Teuchos::VerboseObjectTempState<NLP>
00203     nlpOutputTempState(
00204       nlp_, Teuchos::getFancyOStream(Teuchos::rcp(&out,false)),
00205       Teuchos::VERB_DEFAULT );
00206   // Above, we don't want to litter the output with any trace from the NLP.
00207   // However, if the user forces the verbosity level to be increased, then we
00208   // want to set the stream so that it knows where to print to.
00209     
00210   const size_type
00211     m  = nlp_->m();
00212     
00213   // Get the iteration quantity container objects
00214   IterQuantityAccess<value_type>
00215     &f_iq = obj_f_(s),
00216     &alpha_iq = s.alpha();
00217     
00218   IterQuantityAccess<VectorMutable>
00219     &x_iq   = s.x(),
00220     *c_iq   = m > 0 ? &s.c() : NULL,
00221     *h_iq   = NULL,                   // RAB: Replace latter!
00222     &Gf_iq  = grad_obj_f_(s);
00223 
00224   // check that all the pertinent information is known
00225   if (!s.d().updated_k(0) || !x_iq.updated_k(0))
00226   {
00227     // Dead in the water
00228     TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Error, d_k or x_k not updated." );     
00229     return false;
00230   }
00231     
00232   if (!alpha_iq.updated_k(0) || alpha_iq.get_k(0) > 1 || alpha_iq.get_k(0) <= 0)
00233   {
00234     // if alpha_k is not known then we would need to calculate all the new points
00235     TEUCHOS_TEST_FOR_EXCEPTION( true, std::out_of_range, "Error, alpha_k not updated or out of range [0, 1)." );    
00236     return false;
00237   }
00238 
00239   // Setup some necessary parameters
00240   // Assuming that f_iq, Gf_iq, c_iq, h_iq are updated for k
00241   const value_type Gf_t_dk =
00242     (
00243       Gf_iq.updated_k(0)
00244       ? Gf_iq.get_k(0).inner_product( s.d().get_k(0) )
00245       : dyn_cast<NLPInterfacePack::NLPObjGrad>(*nlp_).calc_Gf_prod(
00246         x_iq.get_k(0),s.d().get_k(0)
00247         )
00248       );
00249   const value_type theta_k = CalculateTheta_k( c_iq, h_iq, 0);
00250   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
00251     out << "\ntheta_k = ||c_k||1/c_k.dim() = " << theta_k << std::endl;
00252   const value_type gamma_f_used = CalculateGammaFUsed(f_iq,theta_k,olevel,out);
00253   const value_type theta_small = theta_small_fact_ * MAX(1.0,theta_k);
00254   const value_type alpha_min = CalculateAlphaMin( gamma_f_used, Gf_t_dk, theta_k, theta_small );
00255 
00256   value_type &alpha_k = alpha_iq.get_k(0);
00257   value_type theta_kp1 = 0.0;;
00258 
00259   // Print out some header/initial information
00260   int w = 15;
00261   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
00262   {
00263     out << "\nBeginning Filter line search method.\n"
00264         << "\nCurrent Filter\n"
00265         << "-----------------------------------------------------\n"
00266         << "|" << setw(25) << "f_with_boundary     " 
00267         << "|" << setw(25) << "theta_with_boundary    "
00268         << "|\n"
00269         << "-----------------------------------------------------\n";
00270 
00271     IterQuantityAccess<Filter_T>& filter_iq = filter_(s);
00272     
00273     if (filter_iq.updated_k(-1)) {
00274       Filter_T& filter = filter_iq.get_k(-1);
00275       if (!filter.empty()) {
00276         for (
00277           Filter_T::iterator entry = filter.begin();
00278           entry != filter.end();
00279           ++entry
00280           )
00281         { 
00282           out << "|" << setw(25) << entry->f
00283               << " " << setw(25) << entry->theta
00284               << "|\n";
00285         }
00286       }
00287       else {
00288         out << "Filter is empty.\n";
00289       }
00290     }
00291     else {
00292       out << "Filter is empty.\n";
00293     }
00294   
00295     // dump header
00296     out << "\n  Iteration Status\n";
00297     out << "----------------------------------------------------------------------------------------------------------\n"; 
00298 
00299     out << "|" << setw(w) << "alpha_k    " 
00300         << "|" << setw(w) << "f_kp1     "
00301         << "|" << setw(w) << "theta_kp1   "
00302         << "|" << setw(w) << "pt. status   "
00303         << "|" << setw(40) << "comment                "
00304         << "|" << std::endl;
00305   
00306     out << "----------------------------------------------------------------------------------------------------------" 
00307         << std::endl;
00308   }
00309 
00310   // Begin the line search
00311   bool augment_filter = false;
00312   bool accepted = false;
00313   while (alpha_k > alpha_min && !accepted)
00314   {
00315     accepted = true;
00316 
00317     // Check that point is safe for calculations (no nans, infs, etc)
00318     if (!ValidatePoint(x_iq, f_iq, c_iq, h_iq, false))
00319     {
00320       accepted = false;
00321       
00322       // Print out some point information
00323       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
00324       {
00325         int w = 15;
00326         // dump point
00327         out << "|" << setw(w) << " --- " 
00328             << " " << setw(w) << " --- "
00329             << " " << setw(w) << " --- "
00330             << " " << setw(w) << " failed "
00331             << " " << setw(40) << " nan_or_inf in calc"
00332             << " " << std::endl;
00333       }
00334   
00335       // Really, we do not need to throw an exception here, we can try and backtrack
00336       // alpha to get into an acceptable region
00337       TEUCHOS_TEST_FOR_EXCEPTION( true, std::out_of_range, "Point Not Valid." );  
00338     }
00339       
00340     // Check if point satisfies filter
00341     if (accepted)
00342     {
00343       theta_kp1 = CalculateTheta_k(c_iq, h_iq, +1);
00344 
00345       // Print out some point information
00346       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
00347       {
00348         // dump point
00349         out << "|" << setw(w) << alpha_k 
00350             << " " << setw(w) << f_iq.get_k(+1)
00351             << " " << setw(w) << theta_kp1;
00352       }
00353 
00354       accepted = CheckFilterAcceptability(f_iq.get_k(+1), theta_kp1, s);
00355 
00356       // Print out failure information
00357       if( !accepted && static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
00358       {
00359         out << " " << setw(w) << "failed"
00360             << " " << setw(40) << "Unacceptable to filter"
00361             << "|" << std::endl;
00362       }
00363     }
00364 
00365     // Check if point has theta less than theta_max
00366     if (accepted) {
00367       if (theta_kp1 > theta_max_) {
00368         accepted = false;
00369         // Print out failure information
00370         if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00371           out << " " << setw(w) << "failed"
00372               << " " << setw(40) << "theta_kp1 > theta_max"
00373               << "|" << std::endl;
00374         }
00375       }
00376     }
00377 
00378 
00379     // Check if point satisfies sufficient decrease (Armijo on f if switching cond holds)
00380     if (accepted)
00381     {
00382       // Check for switching condition
00383       if (ShouldSwitchToArmijo(Gf_t_dk, alpha_k, theta_k, theta_small))
00384       {
00385         accepted = CheckArmijo(Gf_t_dk, alpha_k, f_iq);
00386 
00387         // Print out point information
00388         if(static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
00389         {
00390           if (accepted)
00391           { out << " " << setw(w) << "accepted"; }
00392           else
00393           { out << " " << setw(w) << "failed"; }
00394 
00395           out << " " << setw(40) << "Switch Cond. Holds (Armijo)" << "|" << std::endl;
00396         }
00397       }
00398       else
00399       {
00400         accepted = CheckFractionalReduction(f_iq, gamma_f_used, theta_kp1, theta_k);
00401         if (accepted)
00402         { augment_filter = true; }
00403 
00404         // Print out point information
00405         if(static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
00406         {
00407           if (accepted)
00408           { out << " " << setw(w) << "accepted"; }
00409           else
00410           { out << " " << setw(w) << "failed"; }
00411 
00412           out << " " << setw(40) << "Fraction Reduction (! Switch Cond )" << "|" << std::endl;
00413         }
00414 
00415       }
00416     }
00417 
00418     // if the point fails any of the tests, then backtrack
00419     if (!accepted)
00420     {
00421       // try a smaller alpha_k
00422       alpha_k = alpha_k*back_track_frac_;
00423       UpdatePoint(s.d().get_k(0), alpha_k, x_iq, f_iq, c_iq, h_iq, *nlp_);
00424     }   
00425 
00426   } // end while
00427 
00428 
00429   if (accepted)
00430   {
00431     // Print status
00432     if(static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00433       if (augment_filter)
00434         out << "\nPoint was accepted - augmenting the filter ...\n";
00435       else
00436         out << "Point was accepted - did NOT augment filter.\n";
00437     }
00438 
00439     if (augment_filter) {
00440       AugmentFilter( gamma_f_used, f_iq.get_k(+1), theta_kp1, s, olevel, out );
00441     }
00442     else {
00443       // Just update the filter from the last iteration
00444       UpdateFilter(s);
00445     }
00446   }
00447   else {
00448     // Could not find an acceptable alpha_k, go to restoration
00449     // Print status
00450     //if(static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
00451     //  {
00452     //  out << "\nCould not find acceptable alpha_k - going to restoration phase.\n"; 
00453     //  }
00454     
00455     //TEUCHOS_TEST_FOR_EXCEPTION( true, std::out_of_range, "Tried to go to restoration phase." ); 
00456     
00457     TEUCHOS_TEST_FOR_EXCEPTION( true, LineSearchFailure
00458                         ,"FilterLineSearchFailure : Should go to restoration"
00459       );
00460   }
00461 
00462   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) 
00463   {
00464     out << "\nx_kp1 =\n" << x_iq.get_k(+1);
00465     if (c_iq)
00466     { out << "\nc_kp1 =\n" << c_iq->get_k(+1); }
00467     if (h_iq)
00468     { out << "\nh_kp1 =\n" << h_iq->get_k(+1); }
00469   }
00470 
00471 #if defined(FILTER_DEBUG_OUT)
00472   std::ofstream fout("filter_out.xml", std::ofstream::out | std::ostream::app);
00473   fout << "   <FilterIteration iter=\"" << s.k() << "\">" << std::endl;
00474   fout << "      <SelectedPoint alpha=\"" << alpha_k 
00475        << "\" f=\"" << f_iq.get_k(+1) 
00476        << "\" theta=\"" << theta_kp1
00477        << "\" />" << std::endl;
00478 
00479   // Output the filter
00480   fout << "      <Filter>" << std::endl;
00481     
00482   IterQuantityAccess<Filter_T>& filter_iq = filter_(s);
00483   if (filter_iq.updated_k(0))
00484   {    
00485     Filter_T& current_filter = filter_iq.get_k(0);
00486     for (
00487       Filter_T::iterator entry = current_filter.begin();
00488       entry != current_filter.end();
00489       ++entry
00490       )
00491     {
00492       fout << "         <FilterPoint iter=\"" << entry->iter 
00493            << "\" f=\"" << entry->f 
00494            << "\" theta=\"" << entry->theta << ">" << std::endl;
00495     }
00496   }
00497   else
00498   {
00499     fout << "         <FilterNotUpdated/>" << std::endl;
00500   }
00501 
00502   fout << "      </Filter>" << std::endl;
00503 
00504     
00505   // Output the alpha curve
00506   fout << "      <AlphaCurve>" << std::endl;
00507   value_type alpha_tmp = 1.0;
00508   for (int i = 0; i < 10 || alpha_tmp > alpha_k; ++i)
00509   {
00510     UpdatePoint(s.d().get_k(0), alpha_tmp, x_iq, f_iq, c_iq, h_iq, *nlp_);
00511     if (ValidatePoint(x_iq, f_iq, c_iq, h_iq, false))
00512     {
00513       value_type theta = CalculateTheta_k(c_iq, h_iq, +1);
00514       fout << "         <AlphaPoint "
00515            << "alpha=\"" << alpha_tmp << "\" "
00516            << "f=\"" << f_iq.get_k(+1) << "\" "
00517            << "theta=\"" << theta << ">" << std::endl;
00518     }
00519 
00520     alpha_tmp=alpha_tmp*back_track_frac_;
00521   }
00522 
00523   // restore alpha_k
00524   UpdatePoint(s.d().get_k(0), alpha_k, x_iq, f_iq, c_iq, h_iq, *nlp_);
00525 
00526   fout << "      </AlphaCurve>" << std::endl;
00527 
00528   fout << "   </FilterIteration" << std::endl;
00529 
00530   fout.close();
00531 
00532 #endif
00533     
00534   return true;
00535 }
00536   
00537 void LineSearchFilter_Step::print_step(
00538   const Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00539   ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
00540   ) const
00541 {
00542   out
00543   << L << "*** Filter line search method\n"
00544   << L << "# Assumes initial d_k & alpha_k (0-1) is known and\n"
00545   << L << "# x_kp1, f_kp1, c_kp1 and h_kp1 are calculated for that alpha_k\n"
00546   << L << "Gf_t_dk = <Gf,dk>\n"
00547   << L << "theta_k = norm_1(c_k)/c_k.dim()\n"
00548   << L << "theta_small = theta_small_fact*max(1.0,theta_k)\n"
00549   << L << "if f_min != F_MIN_UNBOUNDED then\n"
00550   << L << "  gamma_f_used = gamma_f * (f_k-f_min)/theta_k\n"
00551   << L << "else\n"
00552   << L
00553   << L << "  gamma_f_used = gamma_f\n"
00554   << L << "end\n"
00555   << L << "if Gf_t_dk < 0 then\n"
00556   << L << "  alpha_min = min(gamma_theta, gamma_f_used*theta_k/(-Gf_t_dk))\n"
00557   << L << "  if theta_k <= theta_small then\n"
00558   << L << "    alpha_min = min(alpha_min, delta_*(theta_k^s_theta)/((-Gf_t_dk)^s_f))\n"
00559   << L << "  end\n"
00560   << L << "else\n"
00561   << L << "  alpha_min = gamma_theta\n"
00562   << L << "end\n"
00563   << L << "alpha_min = alpha_min*gamma_alpha\n"
00564   << L << "# Start the line search\n"
00565   << L << "accepted = false\n"
00566   << L << "augment = false\n"
00567   << L << "while alpha > alpha_min and accepted = false then\n"
00568   << L << "  accepted = true"
00569   << L << "  if any values in x_kp1, f_kp1, c_kp1, h_kp1 are nan or inf then\n"
00570   << L << "    accepted = false\n"
00571   << L << "  end\n"
00572   << L << "  # Check filter\n"
00573   << L << "  if accepted = true then\n"
00574   << L << "    theta_kp1 = norm_1(c_kp1)/c_kp1.dim()\n"
00575   << L << "    for each pt in the filter do\n"
00576   << L << "      if theta_kp1 >= pt.theta and f_kp1 >= pt.f then\n"
00577   << L << "        accepted = false\n"
00578   << L << "        break for\n"
00579   << L << "      end\n"
00580   << L << "    next pt\n"
00581   << L << "  end\n"
00582   << L << "  #Check Sufficient Decrease\n"
00583   << L << "  if accepted = true then"
00584   << L << "    # if switching condition is satisfied, use Armijo on f\n"
00585   << L << "    if theta_k < theta_small and Gf_t_dk < 0 and\n"
00586   << L << "        ((-Gf_t_dk)^s_f)*alpha_k < delta*theta_k^s_theta then\n"
00587   << L << "      if f_kp1 <= f_k + eta_f*alpha_k*Gf_t_dk then\n"
00588   << L << "        accepted = true\n"
00589   << L << "      end\n"
00590   << L << "    else\n"
00591   << L << "      # Verify factional reduction\n"
00592   << L << "      if theta_kp1 < (1-gamma_theta)*theta_k or f_kp1 < f_k - gamma_f_used*theta_k then\n"
00593   << L << "        accepted = true\n"
00594   << L << "        augment = true\n"
00595   << L << "      end\n"
00596   << L << "    end\n"
00597   << L << "  end\n"
00598   << L << "  if accepted = false then\n"
00599   << L << "    alpha_k = alpha_k*back_track_frac\n"
00600   << L << "    x_kp1  = x_k + alpha_k * d_k\n"
00601   << L << "    f_kp1  = f(x_kp1)\n"
00602   << L << "    c_kp1  = c(x_kp1)\n"
00603   << L << "    h_kp1  = h(x_kp1)\n"
00604   << L << "  end\n"    
00605   << L << "end while\n"
00606   << L << "if accepted = true then\n"
00607   << L << "  if augment = true then\n"
00608   << L << "    Augment the filter (use f_with_boundary and theta_with_boundary)\n"
00609   << L << "  end\n"
00610   << L << "else\n"
00611   << L << "  goto the restoration phase\n"
00612   << L << "end\n";
00613 }
00614   
00615 bool LineSearchFilter_Step::ValidatePoint( 
00616   const IterQuantityAccess<VectorMutable>& x,
00617   const IterQuantityAccess<value_type>& f,
00618   const IterQuantityAccess<VectorMutable>* c,
00619   const IterQuantityAccess<VectorMutable>* h,
00620   const bool throw_excpt
00621   ) const
00622 {
00623 
00624   using AbstractLinAlgPack::assert_print_nan_inf;
00625   
00626   if (assert_print_nan_inf(x.get_k(+1), "x", throw_excpt, NULL) 
00627       || assert_print_nan_inf(f.get_k(+1), "f", throw_excpt, NULL)
00628       || (!c || assert_print_nan_inf(c->get_k(+1), "c", throw_excpt, NULL))
00629       || (!h || assert_print_nan_inf(h->get_k(+1), "c", throw_excpt, NULL)))
00630   {
00631     return true;
00632   }
00633   return false;
00634 }
00635   
00636 bool LineSearchFilter_Step::CheckFilterAcceptability( 
00637   const value_type f, 
00638   const value_type theta,
00639   const AlgorithmState& s
00640   ) const
00641 {
00642   bool accepted = true;
00643 
00644   const IterQuantityAccess<Filter_T>& filter_iq = filter_(s);
00645     
00646   if (filter_iq.updated_k(-1))
00647   {
00648     const Filter_T &current_filter = filter_iq.get_k(-1);
00649     
00650     for (
00651       Filter_T::const_iterator entry = current_filter.begin();
00652       entry != current_filter.end();
00653       ++entry
00654       )
00655     { 
00656       if (f >= entry->f && theta >= entry->theta)
00657       {
00658         accepted = false;
00659         break;
00660       }
00661     }
00662   }
00663 
00664   return accepted;
00665 }
00666 
00667 bool LineSearchFilter_Step::CheckArmijo( 
00668   const value_type Gf_t_dk, 
00669   const value_type alpha_k, 
00670   const IterQuantityAccess<value_type>& f_iq
00671   ) const
00672 {
00673   bool accepted = false;
00674 
00675   // Check Armijo on objective fn
00676   double f_kp1 = f_iq.get_k(+1);
00677   double f_k = f_iq.get_k(0);
00678   double lhs = f_k - f_kp1;
00679   double rhs = -eta_f_*alpha_k*Gf_t_dk;
00680   if ( lhs >= rhs )
00681   {
00682     // Accept pt, do NOT augment filter
00683     accepted = true;
00684   }
00685 
00686   return accepted;
00687 }
00688 
00689 bool LineSearchFilter_Step::CheckFractionalReduction( 
00690   const IterQuantityAccess<value_type>& f_iq,
00691   const value_type gamma_f_used,
00692   const value_type theta_kp1, 
00693   const value_type theta_k
00694   ) const
00695 {
00696   bool accepted = false;
00697   if (theta_kp1 <= (1-gamma_theta_)*theta_k
00698       || f_iq.get_k(+1) <= f_iq.get_k(0)-gamma_f_used*theta_k )
00699   {
00700     // Accept pt and augment filter
00701     accepted = true;
00702   }
00703   return accepted;
00704 }
00705 
00706 void LineSearchFilter_Step::UpdatePoint( 
00707   const VectorMutable& d,
00708   const value_type alpha, 
00709   IterQuantityAccess<VectorMutable> &x,
00710   IterQuantityAccess<value_type>& f,
00711   IterQuantityAccess<VectorMutable>* c,
00712   IterQuantityAccess<VectorMutable>* h,
00713   NLP& nlp ) const
00714 {  
00715   using LinAlgOpPack::Vp_StV;
00716   using AbstractLinAlgPack::assert_print_nan_inf;
00717   VectorMutable& x_kp1 = x.set_k(+1);
00718   x_kp1 = x.get_k(0);
00719   Vp_StV( &x_kp1, alpha, d);
00720 
00721   if (assert_print_nan_inf(x_kp1, "x", true, NULL))
00722   {
00723     // Calcuate f and c at the new point.
00724     nlp.unset_quantities();
00725     nlp.set_f( &f.set_k(+1) );
00726     if (c) nlp.set_c( &c->set_k(+1) );
00727     nlp.calc_f( x_kp1 ); 
00728     if (c) nlp.calc_c( x_kp1, false );
00729     nlp.unset_quantities();
00730   }
00731 }
00732 
00733 value_type LineSearchFilter_Step::CalculateAlphaMin( 
00734   const value_type gamma_f_used,
00735   const value_type Gf_t_dk,
00736   const value_type theta_k, 
00737   const value_type theta_small
00738   ) const
00739 {
00740   value_type alpha_min = 0;
00741     
00742   if (Gf_t_dk < 0)
00743   {
00744     alpha_min = MIN(gamma_theta_, gamma_f_used*theta_k/(-Gf_t_dk));
00745     if (theta_k <= theta_small)
00746     {
00747       value_type switch_bound = delta_*pow(theta_k, s_theta_)/pow(-Gf_t_dk,s_f_);
00748       alpha_min = MIN(alpha_min, switch_bound);
00749     }
00750   }
00751   else
00752   {
00753     alpha_min = gamma_theta_;
00754   }
00755   return alpha_min * gamma_alpha_;
00756 }
00757 
00758 value_type LineSearchFilter_Step::CalculateTheta_k( 
00759   const IterQuantityAccess<VectorMutable>* c,
00760   const IterQuantityAccess<VectorMutable>* h,
00761   int k
00762   ) const
00763 {
00764   value_type theta = 0.0;
00765   if (h) {
00766     TEUCHOS_TEST_FOR_EXCEPTION( true, std::out_of_range, "Error, do not support inequalities yet" );
00767   }
00768   if (c) {
00769     const Vector &c_k = c->get_k(k);
00770     theta += ( c_k.norm_1() / c_k.dim() );
00771   }
00772   return theta;
00773 }
00774 
00775 value_type LineSearchFilter_Step::CalculateGammaFUsed(
00776   const IterQuantityAccess<value_type> &f,
00777   const value_type theta_k,
00778   const EJournalOutputLevel olevel,
00779   std::ostream &out
00780   ) const
00781 {
00782   if( f_min() == F_MIN_UNBOUNDED) {
00783     if((int)olevel >= static_cast<int>(PRINT_ALGORITHM_STEPS))
00784       out << "\nf_min==F_MIN_UNBOUNDED:  Setting gamma_f_used = gamma_f = "<<gamma_f()<<".\n";
00785     return gamma_f();
00786   }
00787   const value_type f_k = f.get_k(0);
00788   if( f_k < f_min() ) {
00789     if((int)olevel >= static_cast<int>(PRINT_ALGORITHM_STEPS))
00790       out << "\nWarning!!! f_k = "<<f_k<<" < f_min = "<<f_min()<<"!  Setting gamma_f_used = gamma_f = "<<gamma_f()<<".\n";
00791     return gamma_f();
00792   }
00793   const value_type gamma_f_used = gamma_f() * ( f_k - f_min() ) / theta_k;
00794   if((int)olevel >= static_cast<int>(PRINT_ALGORITHM_STEPS))
00795     out << "\nf_min = "<<f_min()<<"!=F_MIN_UNBOUNDED: Setting gamma_f_used = gamma_f * (f_k - f_min)/theta_k = "
00796         <<gamma_f()<<" * ("<<f_k<<" - "<<f_min()<<")/"<<theta_k<<" = "<<gamma_f_used <<".\n";
00797   return gamma_f_used;
00798 }
00799 
00800 bool LineSearchFilter_Step::ShouldSwitchToArmijo( 
00801   const value_type Gf_t_dk,
00802   const value_type alpha_k,
00803   const value_type theta_k,
00804   const value_type theta_small
00805   ) const
00806 {
00807   if (theta_k < theta_small && Gf_t_dk < 0) {
00808     if (pow(-Gf_t_dk, s_f_)*alpha_k - delta_*pow(theta_k, s_theta_) > 0) {
00809       return true;
00810     }
00811   }
00812   return false;
00813 }
00814 
00815 
00816 void LineSearchFilter_Step::UpdateFilter( IterationPack::AlgorithmState& s ) const
00817 {
00818   
00819   IterQuantityAccess<Filter_T>& filter_iq = filter_(s);
00820     
00821   if (!filter_iq.updated_k(0))
00822   {
00823     if (filter_iq.updated_k(-1))
00824     {
00825       // initialize the filter from the last iteration
00826       filter_iq.set_k(0,-1);
00827     }
00828     else
00829     {
00830       // create an uninitialized filter
00831       filter_iq.set_k(0);
00832     }
00833   }
00834 
00835 }
00836 
00837 void LineSearchFilter_Step::AugmentFilter( 
00838   const value_type gamma_f_used,
00839   const value_type f,
00840   const value_type theta,
00841   IterationPack::AlgorithmState& s,
00842   const EJournalOutputLevel olevel,
00843   std::ostream &out
00844   ) const
00845 {
00846 
00847   using Teuchos::as;
00848 
00849   const value_type
00850     f_with_boundary = f-gamma_f_used*theta,
00851     theta_with_boundary = (1.0-gamma_theta())*theta;
00852 
00853   if(static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00854     out << "\nAugmenting the filter with the point:"
00855         << "\n  f_with_boundary = f_kp1 - gamma_f_used*theta_kp1 = "<<f<<" - "<<gamma_f_used<<"*"<<theta<< " = " <<f_with_boundary
00856         << "\n  theta_with_boundary = (1-gamma_theta)*theta_kp1 = (1-"<<gamma_theta()<<")*"<<theta<<" = "<<theta_with_boundary
00857         << std::endl;
00858   }
00859   
00860   UpdateFilter(s);
00861   Filter_T& current_filter = filter_(s).get_k(0);
00862 
00863   // Remove current filter entries that are not longer needed
00864   current_filter.remove_if(
00865     RemoveFilterEntryPred(
00866       f_with_boundary, theta_with_boundary,
00867       as<int>(olevel) >= as<int>(PRINT_ALGORITHM_STEPS)
00868       ? Teuchos::rcpFromRef(out) : Teuchos::null
00869       )
00870     );
00871 
00872   // Now append the current point
00873   current_filter.push_front(FilterEntry(f_with_boundary, theta_with_boundary, s.k()));
00874 
00875 }
00876 
00877 // static
00878 
00879 value_type LineSearchFilter_Step::F_MIN_UNBOUNDED = std::numeric_limits<value_type>::min();
00880 
00881 } // end namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends