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

Generated on Thu Sep 18 12:35:17 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1