MoochoPack_MeritFunc_PenaltyParamsUpdateWithMult_AddedStep.cpp

Go to the documentation of this file.
00001 #if 0
00002 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00007 //                  Copyright (2003) Sandia Corporation
00008 // 
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 // 
00012 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00027 // 
00028 // ***********************************************************************
00029 // @HEADER
00030 
00031 #include <ostream>
00032 #include <typeinfo>
00033 
00034 #include "MoochoPack_MeritFunc_PenaltyParamsUpdateWithMult_AddedStep.hpp"
00035 #include "MoochoPack_moocho_algo_conversion.hpp"
00036 #include "IterationPack_print_algorithm_step.hpp"
00037 #include "ConstrainedOptPack_MeritFuncPenaltyParams.hpp"
00038 #include "ConstrainedOptPack_MeritFuncNLPDirecDeriv.hpp"
00039 #include "ConstrainedOptPack/src/VectorWithNorms.h"
00040 #include "DenseLinAlgPack_DVectorOp.hpp"
00041 #include "DenseLinAlgPack_DVectorClass.hpp"
00042 #include "DenseLinAlgPack_DVectorOut.hpp"
00043 
00044 namespace {
00045 
00046 typedef MoochoPack::value_type value_type;
00047 inline value_type max(value_type v1, value_type v2)
00048 { return (v1 > v2) ? v1 : v2; }
00049 
00050 }
00051 
00052 namespace MoochoPack {
00053 
00054 MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::MeritFunc_PenaltyParamsUpdateWithMult_AddedStep(
00055       const merit_func_ptr_t& merit_func, value_type small_mu, value_type min_mu_ratio
00056     , value_type mult_factor, value_type kkt_near_sol )
00057   : merit_func_(merit_func), near_solution_(false)
00058     , small_mu_(small_mu), min_mu_ratio_(min_mu_ratio), mult_factor_(mult_factor)
00059     , kkt_near_sol_(kkt_near_sol), norm_inf_mu_last_(0.0)
00060 {}
00061 
00062 bool MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::do_step(Algorithm& _algo
00063   , poss_type step_poss, IterationPack::EDoStepType type
00064   , poss_type assoc_step_poss)
00065 {
00066   using DenseLinAlgPack::norm_inf;
00067 
00068   NLPAlgo &algo = rsqp_algo(_algo);
00069   NLPAlgoState  &s    = algo.rsqp_state();
00070   
00071   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00072   std::ostream& out = algo.track().journal_out();
00073 
00074   // print step header.
00075   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00076     using IterationPack::print_algorithm_step;
00077     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00078   }
00079 
00080   MeritFuncPenaltyParams
00081     *params = dynamic_cast<MeritFuncPenaltyParams*>(&merit_func());
00082   if( !params ) {
00083     std::ostringstream omsg;
00084     omsg
00085       << "MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::do_step(...), Error "
00086       << "The class " << typeName(&merit_func()) << " does not support the "
00087       << "MeritFuncPenaltyParams iterface\n";
00088     out << omsg.str();
00089     throw std::logic_error( omsg.str() );
00090   }
00091 
00092   MeritFuncNLPDirecDeriv
00093     *direc_deriv = dynamic_cast<MeritFuncNLPDirecDeriv*>(&merit_func());
00094   if( !direc_deriv ) {
00095     std::ostringstream omsg;
00096     omsg
00097       << "MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::do_step(...), Error "
00098       << "The class " << typeName(&merit_func()) << " does not support the "
00099       << "MeritFuncNLPDirecDeriv iterface\n";
00100     out << omsg.str();
00101     throw std::logic_error( omsg.str() );
00102   }
00103 
00104   bool perform_update = true;
00105 
00106   if( s.mu().updated_k(0) ) {
00107     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00108       out << "\nmu_k is already updated by someone else?\n";
00109     }
00110     const value_type mu_k = s.mu().get_k(0);
00111     if( mu_k == norm_inf_mu_last_ ) {
00112       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00113         out << "\nmu_k " << mu_k << " == norm_inf_mu_last = " << norm_inf_mu_last_
00114           << "\nso we will take this as a signal to skip the update.\n";
00115       }
00116       perform_update = false;
00117     }
00118     else {
00119       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00120         out << "\nmu_k " << mu_k << " != norm_inf_mu_last = " << norm_inf_mu_last_
00121           << "\nso we will ignore this and perform the update anyway.\n";
00122       }
00123     }   
00124   }
00125   if(perform_update) {
00126 
00127     if ( s.lambda().updated_k(0) ) {
00128 
00129       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00130         out << "\nUpdate the penalty parameter...\n";
00131       }
00132 
00133       const DVector
00134         &lambda_k = s.lambda().get_k(0).cv();
00135 
00136       if( params->mu().size() != lambda_k.size() )
00137         params->resize( lambda_k.size() );
00138       DVectorSlice
00139         mu = params->mu();
00140 
00141       const value_type
00142         max_lambda  = norm_inf( lambda_k() ),
00143         mult_fact = (1.0 + mult_factor_);
00144 
00145       if(near_solution_) {
00146         if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00147           out << "\nNear solution, forcing mu(j) >= mu_old(j)...\n";
00148         }
00149         DVector::const_iterator lb_itr = lambda_k.begin();
00150         DVectorSlice::iterator  mu_itr = mu.begin();
00151         for( ; lb_itr != lambda_k.end(); ++mu_itr, ++ lb_itr )
00152           *mu_itr = max( max( *mu_itr, mult_fact * ::fabs(*lb_itr) ), small_mu_ );
00153       }
00154       else {
00155         if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00156           out << "\nNot near solution, allowing reduction in mu(j) ...\n";
00157         }
00158         DVector::const_iterator lb_itr = lambda_k.begin();
00159         DVectorSlice::iterator  mu_itr = mu.begin();
00160         for( ; lb_itr != lambda_k.end(); ++mu_itr, ++ lb_itr ) {
00161           const value_type lb_j = ::fabs(*lb_itr);
00162           *mu_itr = max(
00163                   (3.0 * (*mu_itr) + lb_j) / 4.0  
00164                 , max( mult_fact * lb_j, small_mu_ )
00165                 );
00166         }
00167         value_type kkt_error = s.opt_kkt_err().get_k(0) + s.feas_kkt_err().get_k(0);
00168         if(kkt_error <= kkt_near_sol_) {
00169           if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00170             out << "\nkkt_error = " << kkt_error << " <= kkt_near_sol = "
00171                 << kkt_near_sol_ << std::endl
00172               << "Switching to forcing mu_k >= mu_km1 in the future\n";
00173           }
00174           near_solution_ = true;
00175         }
00176       }
00177 
00178       // Force the ratio
00179       const value_type
00180           max_mu  = norm_inf( mu() ),
00181           min_mu  = min_mu_ratio_ * max_mu;
00182       for(DVectorSlice::iterator mu_itr = mu.begin(); mu_itr != mu.end(); ++mu_itr)
00183         *mu_itr = max( (*mu_itr), min_mu ); 
00184 
00185       s.mu().set_k(0) = norm_inf_mu_last_ = max_mu;
00186 
00187       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00188         out << "\nmax(|mu(j)|) = " << (*std::max_element( mu.begin(), mu.end() ))
00189           << "\nmin(|mu(j)|) = " << (*std::min_element( mu.begin(), mu.end() ))
00190             << std::endl;
00191       }
00192 
00193       if( (int)olevel >= (int)PRINT_VECTORS ) {
00194         out << "\nmu = \n" << mu;
00195       }
00196     }
00197     else {
00198       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00199         out << "\nDon't have the info to update penalty parameter so just use the last updated...\n";
00200       }
00201     }
00202   }
00203 
00204   // In addition also compute the directional derivative
00205   direc_deriv->calc_deriv( s.Gf().get_k(0)(), s.c().get_k(0)(), s.d().get_k(0)() );
00206 
00207   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00208     out << "\nmu_k = " << s.mu().get_k(0) << "\n";
00209   }
00210 
00211   return true;
00212 }
00213 
00214 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::print_step( const Algorithm& algo
00215   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00216   , std::ostream& out, const std::string& L ) const
00217 {
00218   out
00219     << L << "*** Update the penalty parameter for the merit function to ensure\n"
00220     << L << "*** a descent direction a directional derivatieve.\n"
00221     << L << "*** phi is a merit function object that uses the penalty parameter mu.\n"
00222     << L << "default: near_solution = false\n"
00223     << L << "         small_mu = " << small_mu_ << std::endl
00224     << L << "         min_mu_ratio = " << min_mu_ratio_ << std::endl
00225     << L << "         mult_factor = " << mult_factor_ << std::endl
00226     << L << "         kkt_near_sol = " << kkt_near_sol_ << std::endl
00227     << L << "perform_update = true\n"
00228     << L << "if mu_k is already updated then\n"
00229     << L << "    if mu_k == norm_inf_mu_last then\n"
00230     << L << "        *** We will use this as a signal to skip the update\n"
00231     << L << "        perform_update = false\n"
00232     << L << "    else\n"
00233     << L << "        *** We will perform the update anyway\n"
00234     << L << "    end\n"
00235     << L << "if perform_update == true then\n"
00236     << L << "    if lambda_k is updated then\n"
00237     << L << "        max_lambda = norm(lambda_k,inf)\n"
00238     << L << "        mult_fact = (1+mult_factor)\n"
00239     << L << "        mu = phi.mu()\n"
00240     << L << "        if near_solution == true\n"
00241     << L << "            for j = 1...m\n"
00242     << L << "                mu(j) = max(max(mu(j),mult_fact*abs(lambda_k(j))),small_mu)\n"
00243     << L << "            end\n"
00244     << L << "        else\n"
00245     << L << "            for j = 1...m\n"
00246     << L << "                mu(j) = max( ( 3.0 * mu(j) + abs(lambda_k(j)) ) / 4.0\n"
00247     << L << "                            , max( 1.001 * abs(lambda_k(j)) , small_mu ) )\n"
00248     << L << "            end\n"
00249     << L << "            kkt_error = opt_kkt_err_k + feas_kkt_err_k\n"
00250     << L << "            if kkt_error <= kkt_near_sol then\n"
00251     << L << "                near_solution = true\n"
00252     << L << "            end\n"
00253     << L << "        end\n"
00254     << L << "        min_mu = min_mu_ratio * norm(mu,inf)\n"
00255     << L << "        for j = 1...m\n"
00256     << L << "            mu(j) = max( mu(j), min_mu )\n"
00257     << L << "        end\n"
00258     << L << "    else\n"
00259     << L << "        *** Don't have the information to perform the update.\n"
00260     << L << "    end\n"
00261     << L << "end\n"
00262     << L << "phi.calc_deriv(Gf_k,c_k,d_k)\n";
00263 }
00264 
00265 // Overridden from MeritFunc_PenaltyParamUpdate_AddedStep
00266 
00267 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::small_mu( value_type small_mu )
00268 {
00269   small_mu_ = small_mu;
00270 }
00271 
00272 value_type MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::small_mu() const
00273 {
00274   return small_mu_;
00275 }
00276 
00277 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::min_mu_ratio( value_type min_mu_ratio )
00278 {
00279   min_mu_ratio_ = min_mu_ratio;
00280 }
00281 
00282 value_type MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::min_mu_ratio() const
00283 {
00284   return min_mu_ratio_;
00285 }
00286 
00287 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::mult_factor( value_type mult_factor )
00288 {
00289   mult_factor_ = mult_factor;
00290 }
00291 
00292 value_type MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::mult_factor() const
00293 {
00294   return mult_factor_;
00295 }
00296 
00297 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::kkt_near_sol( value_type kkt_near_sol )
00298 {
00299   kkt_near_sol_ = kkt_near_sol;
00300 }
00301 
00302 value_type MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::kkt_near_sol() const
00303 {
00304   return kkt_near_sol_;
00305 }
00306 
00307 
00308 } // end namespace MoochoPack
00309 
00310 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:10:59 2011 for MOOCHO (Single Doxygen Collection) by  doxygen 1.6.3