MoochoPack_MeritFunc_PenaltyParamsUpdateWithMult_AddedStep.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 <ostream>
00030 #include <typeinfo>
00031 
00032 #include "../std/MoochoPack_MeritFunc_PenaltyParamsUpdateWithMult_AddedStep.hpp"
00033 #include "../MoochoPack_moocho_algo_conversion.hpp"
00034 #include "IterationPack_print_algorithm_step.hpp"
00035 #include "ConstrainedOptPack_MeritFuncPenaltyParams.hpp"
00036 #include "ConstrainedOptPack_MeritFuncNLPDirecDeriv.hpp"
00037 #include "ConstrainedOptPack/src/VectorWithNorms.h"
00038 #include "DenseLinAlgPack_DVectorOp.hpp"
00039 #include "DenseLinAlgPack_DVectorClass.hpp"
00040 #include "DenseLinAlgPack_DVectorOut.hpp"
00041 
00042 namespace {
00043 
00044 typedef MoochoPack::value_type value_type;
00045 inline value_type max(value_type v1, value_type v2)
00046 { return (v1 > v2) ? v1 : v2; }
00047 
00048 }
00049 
00050 namespace MoochoPack {
00051 
00052 MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::MeritFunc_PenaltyParamsUpdateWithMult_AddedStep(
00053       const merit_func_ptr_t& merit_func, value_type small_mu, value_type min_mu_ratio
00054     , value_type mult_factor, value_type kkt_near_sol )
00055   : merit_func_(merit_func), near_solution_(false)
00056     , small_mu_(small_mu), min_mu_ratio_(min_mu_ratio), mult_factor_(mult_factor)
00057     , kkt_near_sol_(kkt_near_sol), norm_inf_mu_last_(0.0)
00058 {}
00059 
00060 bool MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::do_step(Algorithm& _algo
00061   , poss_type step_poss, IterationPack::EDoStepType type
00062   , poss_type assoc_step_poss)
00063 {
00064   using DenseLinAlgPack::norm_inf;
00065 
00066   NLPAlgo &algo = rsqp_algo(_algo);
00067   NLPAlgoState  &s    = algo.rsqp_state();
00068   
00069   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00070   std::ostream& out = algo.track().journal_out();
00071 
00072   // print step header.
00073   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00074     using IterationPack::print_algorithm_step;
00075     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00076   }
00077 
00078   MeritFuncPenaltyParams
00079     *params = dynamic_cast<MeritFuncPenaltyParams*>(&merit_func());
00080   if( !params ) {
00081     std::ostringstream omsg;
00082     omsg
00083       << "MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::do_step(...), Error "
00084       << "The class " << typeid(&merit_func()).name() << " does not support the "
00085       << "MeritFuncPenaltyParams iterface\n";
00086     out << omsg.str();
00087     throw std::logic_error( omsg.str() );
00088   }
00089 
00090   MeritFuncNLPDirecDeriv
00091     *direc_deriv = dynamic_cast<MeritFuncNLPDirecDeriv*>(&merit_func());
00092   if( !direc_deriv ) {
00093     std::ostringstream omsg;
00094     omsg
00095       << "MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::do_step(...), Error "
00096       << "The class " << typeid(&merit_func()).name() << " does not support the "
00097       << "MeritFuncNLPDirecDeriv iterface\n";
00098     out << omsg.str();
00099     throw std::logic_error( omsg.str() );
00100   }
00101 
00102   bool perform_update = true;
00103 
00104   if( s.mu().updated_k(0) ) {
00105     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00106       out << "\nmu_k is already updated by someone else?\n";
00107     }
00108     const value_type mu_k = s.mu().get_k(0);
00109     if( mu_k == norm_inf_mu_last_ ) {
00110       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00111         out << "\nmu_k " << mu_k << " == norm_inf_mu_last = " << norm_inf_mu_last_
00112           << "\nso we will take this as a signal to skip the update.\n";
00113       }
00114       perform_update = false;
00115     }
00116     else {
00117       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00118         out << "\nmu_k " << mu_k << " != norm_inf_mu_last = " << norm_inf_mu_last_
00119           << "\nso we will ignore this and perform the update anyway.\n";
00120       }
00121     }   
00122   }
00123   if(perform_update) {
00124 
00125     if ( s.lambda().updated_k(0) ) {
00126 
00127       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00128         out << "\nUpdate the penalty parameter...\n";
00129       }
00130 
00131       const DVector
00132         &lambda_k = s.lambda().get_k(0).cv();
00133 
00134       if( params->mu().size() != lambda_k.size() )
00135         params->resize( lambda_k.size() );
00136       DVectorSlice
00137         mu = params->mu();
00138 
00139       const value_type
00140         max_lambda  = norm_inf( lambda_k() ),
00141         mult_fact = (1.0 + mult_factor_);
00142 
00143       if(near_solution_) {
00144         if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00145           out << "\nNear solution, forcing mu(j) >= mu_old(j)...\n";
00146         }
00147         DVector::const_iterator lb_itr = lambda_k.begin();
00148         DVectorSlice::iterator  mu_itr = mu.begin();
00149         for( ; lb_itr != lambda_k.end(); ++mu_itr, ++ lb_itr )
00150           *mu_itr = max( max( *mu_itr, mult_fact * ::fabs(*lb_itr) ), small_mu_ );
00151       }
00152       else {
00153         if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00154           out << "\nNot near solution, allowing reduction in mu(j) ...\n";
00155         }
00156         DVector::const_iterator lb_itr = lambda_k.begin();
00157         DVectorSlice::iterator  mu_itr = mu.begin();
00158         for( ; lb_itr != lambda_k.end(); ++mu_itr, ++ lb_itr ) {
00159           const value_type lb_j = ::fabs(*lb_itr);
00160           *mu_itr = max(
00161                   (3.0 * (*mu_itr) + lb_j) / 4.0  
00162                 , max( mult_fact * lb_j, small_mu_ )
00163                 );
00164         }
00165         value_type kkt_error = s.opt_kkt_err().get_k(0) + s.feas_kkt_err().get_k(0);
00166         if(kkt_error <= kkt_near_sol_) {
00167           if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00168             out << "\nkkt_error = " << kkt_error << " <= kkt_near_sol = "
00169                 << kkt_near_sol_ << std::endl
00170               << "Switching to forcing mu_k >= mu_km1 in the future\n";
00171           }
00172           near_solution_ = true;
00173         }
00174       }
00175 
00176       // Force the ratio
00177       const value_type
00178           max_mu  = norm_inf( mu() ),
00179           min_mu  = min_mu_ratio_ * max_mu;
00180       for(DVectorSlice::iterator mu_itr = mu.begin(); mu_itr != mu.end(); ++mu_itr)
00181         *mu_itr = max( (*mu_itr), min_mu ); 
00182 
00183       s.mu().set_k(0) = norm_inf_mu_last_ = max_mu;
00184 
00185       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00186         out << "\nmax(|mu(j)|) = " << (*std::max_element( mu.begin(), mu.end() ))
00187           << "\nmin(|mu(j)|) = " << (*std::min_element( mu.begin(), mu.end() ))
00188             << std::endl;
00189       }
00190 
00191       if( (int)olevel >= (int)PRINT_VECTORS ) {
00192         out << "\nmu = \n" << mu;
00193       }
00194     }
00195     else {
00196       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00197         out << "\nDon't have the info to update penalty parameter so just use the last updated...\n";
00198       }
00199     }
00200   }
00201 
00202   // In addition also compute the directional derivative
00203   direc_deriv->calc_deriv( s.Gf().get_k(0)(), s.c().get_k(0)(), s.d().get_k(0)() );
00204 
00205   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00206     out << "\nmu_k = " << s.mu().get_k(0) << "\n";
00207   }
00208 
00209   return true;
00210 }
00211 
00212 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::print_step( const Algorithm& algo
00213   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00214   , std::ostream& out, const std::string& L ) const
00215 {
00216   out
00217     << L << "*** Update the penalty parameter for the merit function to ensure\n"
00218     << L << "*** a descent direction a directional derivatieve.\n"
00219     << L << "*** phi is a merit function object that uses the penalty parameter mu.\n"
00220     << L << "default: near_solution = false\n"
00221     << L << "         small_mu = " << small_mu_ << std::endl
00222     << L << "         min_mu_ratio = " << min_mu_ratio_ << std::endl
00223     << L << "         mult_factor = " << mult_factor_ << std::endl
00224     << L << "         kkt_near_sol = " << kkt_near_sol_ << std::endl
00225     << L << "perform_update = true\n"
00226     << L << "if mu_k is already updated then\n"
00227     << L << "    if mu_k == norm_inf_mu_last then\n"
00228     << L << "        *** We will use this as a signal to skip the update\n"
00229     << L << "        perform_update = false\n"
00230     << L << "    else\n"
00231     << L << "        *** We will perform the update anyway\n"
00232     << L << "    end\n"
00233     << L << "if perform_update == true then\n"
00234     << L << "    if lambda_k is updated then\n"
00235     << L << "        max_lambda = norm(lambda_k,inf)\n"
00236     << L << "        mult_fact = (1+mult_factor)\n"
00237     << L << "        mu = phi.mu()\n"
00238     << L << "        if near_solution == true\n"
00239     << L << "            for j = 1...m\n"
00240     << L << "                mu(j) = max(max(mu(j),mult_fact*abs(lambda_k(j))),small_mu)\n"
00241     << L << "            end\n"
00242     << L << "        else\n"
00243     << L << "            for j = 1...m\n"
00244     << L << "                mu(j) = max( ( 3.0 * mu(j) + abs(lambda_k(j)) ) / 4.0\n"
00245     << L << "                            , max( 1.001 * abs(lambda_k(j)) , small_mu ) )\n"
00246     << L << "            end\n"
00247     << L << "            kkt_error = opt_kkt_err_k + feas_kkt_err_k\n"
00248     << L << "            if kkt_error <= kkt_near_sol then\n"
00249     << L << "                near_solution = true\n"
00250     << L << "            end\n"
00251     << L << "        end\n"
00252     << L << "        min_mu = min_mu_ratio * norm(mu,inf)\n"
00253     << L << "        for j = 1...m\n"
00254     << L << "            mu(j) = max( mu(j), min_mu )\n"
00255     << L << "        end\n"
00256     << L << "    else\n"
00257     << L << "        *** Don't have the information to perform the update.\n"
00258     << L << "    end\n"
00259     << L << "end\n"
00260     << L << "phi.calc_deriv(Gf_k,c_k,d_k)\n";
00261 }
00262 
00263 // Overridden from MeritFunc_PenaltyParamUpdate_AddedStep
00264 
00265 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::small_mu( value_type small_mu )
00266 {
00267   small_mu_ = small_mu;
00268 }
00269 
00270 value_type MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::small_mu() const
00271 {
00272   return small_mu_;
00273 }
00274 
00275 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::min_mu_ratio( value_type min_mu_ratio )
00276 {
00277   min_mu_ratio_ = min_mu_ratio;
00278 }
00279 
00280 value_type MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::min_mu_ratio() const
00281 {
00282   return min_mu_ratio_;
00283 }
00284 
00285 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::mult_factor( value_type mult_factor )
00286 {
00287   mult_factor_ = mult_factor;
00288 }
00289 
00290 value_type MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::mult_factor() const
00291 {
00292   return mult_factor_;
00293 }
00294 
00295 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::kkt_near_sol( value_type kkt_near_sol )
00296 {
00297   kkt_near_sol_ = kkt_near_sol;
00298 }
00299 
00300 value_type MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::kkt_near_sol() const
00301 {
00302   return kkt_near_sol_;
00303 }
00304 
00305 
00306 } // end namespace MoochoPack

Generated on Thu Sep 18 12:34:28 2008 for MoochoPack : Framework for Large-Scale Optimization Algorithms by doxygen 1.3.9.1