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

Generated on Tue Oct 20 12:51:47 2009 for MOOCHO (Single Doxygen Collection) by doxygen 1.4.7