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

Generated on Wed May 12 21:51:20 2010 for MoochoPack : Framework for Large-Scale Optimization Algorithms by  doxygen 1.4.7