MOOCHO (Single Doxygen Collection) Version of the Day
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 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission.
00026 //
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //
00039 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00040 // 
00041 // ***********************************************************************
00042 // @HEADER
00043 
00044 #include <ostream>
00045 #include <typeinfo>
00046 
00047 #include "MoochoPack_MeritFunc_PenaltyParamsUpdateWithMult_AddedStep.hpp"
00048 #include "MoochoPack_moocho_algo_conversion.hpp"
00049 #include "IterationPack_print_algorithm_step.hpp"
00050 #include "ConstrainedOptPack_MeritFuncPenaltyParams.hpp"
00051 #include "ConstrainedOptPack_MeritFuncNLPDirecDeriv.hpp"
00052 #include "ConstrainedOptPack/src/VectorWithNorms.h"
00053 #include "DenseLinAlgPack_DVectorOp.hpp"
00054 #include "DenseLinAlgPack_DVectorClass.hpp"
00055 #include "DenseLinAlgPack_DVectorOut.hpp"
00056 
00057 namespace {
00058 
00059 typedef MoochoPack::value_type value_type;
00060 inline value_type max(value_type v1, value_type v2)
00061 { return (v1 > v2) ? v1 : v2; }
00062 
00063 }
00064 
00065 namespace MoochoPack {
00066 
00067 MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::MeritFunc_PenaltyParamsUpdateWithMult_AddedStep(
00068       const merit_func_ptr_t& merit_func, value_type small_mu, value_type min_mu_ratio
00069     , value_type mult_factor, value_type kkt_near_sol )
00070   : merit_func_(merit_func), near_solution_(false)
00071     , small_mu_(small_mu), min_mu_ratio_(min_mu_ratio), mult_factor_(mult_factor)
00072     , kkt_near_sol_(kkt_near_sol), norm_inf_mu_last_(0.0)
00073 {}
00074 
00075 bool MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::do_step(Algorithm& _algo
00076   , poss_type step_poss, IterationPack::EDoStepType type
00077   , poss_type assoc_step_poss)
00078 {
00079   using DenseLinAlgPack::norm_inf;
00080 
00081   NLPAlgo &algo = rsqp_algo(_algo);
00082   NLPAlgoState  &s    = algo.rsqp_state();
00083   
00084   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00085   std::ostream& out = algo.track().journal_out();
00086 
00087   // print step header.
00088   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00089     using IterationPack::print_algorithm_step;
00090     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00091   }
00092 
00093   MeritFuncPenaltyParams
00094     *params = dynamic_cast<MeritFuncPenaltyParams*>(&merit_func());
00095   if( !params ) {
00096     std::ostringstream omsg;
00097     omsg
00098       << "MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::do_step(...), Error "
00099       << "The class " << typeName(&merit_func()) << " does not support the "
00100       << "MeritFuncPenaltyParams iterface\n";
00101     out << omsg.str();
00102     throw std::logic_error( omsg.str() );
00103   }
00104 
00105   MeritFuncNLPDirecDeriv
00106     *direc_deriv = dynamic_cast<MeritFuncNLPDirecDeriv*>(&merit_func());
00107   if( !direc_deriv ) {
00108     std::ostringstream omsg;
00109     omsg
00110       << "MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::do_step(...), Error "
00111       << "The class " << typeName(&merit_func()) << " does not support the "
00112       << "MeritFuncNLPDirecDeriv iterface\n";
00113     out << omsg.str();
00114     throw std::logic_error( omsg.str() );
00115   }
00116 
00117   bool perform_update = true;
00118 
00119   if( s.mu().updated_k(0) ) {
00120     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00121       out << "\nmu_k is already updated by someone else?\n";
00122     }
00123     const value_type mu_k = s.mu().get_k(0);
00124     if( mu_k == norm_inf_mu_last_ ) {
00125       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00126         out << "\nmu_k " << mu_k << " == norm_inf_mu_last = " << norm_inf_mu_last_
00127           << "\nso we will take this as a signal to skip the update.\n";
00128       }
00129       perform_update = false;
00130     }
00131     else {
00132       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00133         out << "\nmu_k " << mu_k << " != norm_inf_mu_last = " << norm_inf_mu_last_
00134           << "\nso we will ignore this and perform the update anyway.\n";
00135       }
00136     }   
00137   }
00138   if(perform_update) {
00139 
00140     if ( s.lambda().updated_k(0) ) {
00141 
00142       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00143         out << "\nUpdate the penalty parameter...\n";
00144       }
00145 
00146       const DVector
00147         &lambda_k = s.lambda().get_k(0).cv();
00148 
00149       if( params->mu().size() != lambda_k.size() )
00150         params->resize( lambda_k.size() );
00151       DVectorSlice
00152         mu = params->mu();
00153 
00154       const value_type
00155         max_lambda  = norm_inf( lambda_k() ),
00156         mult_fact = (1.0 + mult_factor_);
00157 
00158       if(near_solution_) {
00159         if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00160           out << "\nNear solution, forcing mu(j) >= mu_old(j)...\n";
00161         }
00162         DVector::const_iterator lb_itr = lambda_k.begin();
00163         DVectorSlice::iterator  mu_itr = mu.begin();
00164         for( ; lb_itr != lambda_k.end(); ++mu_itr, ++ lb_itr )
00165           *mu_itr = max( max( *mu_itr, mult_fact * ::fabs(*lb_itr) ), small_mu_ );
00166       }
00167       else {
00168         if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00169           out << "\nNot near solution, allowing reduction in mu(j) ...\n";
00170         }
00171         DVector::const_iterator lb_itr = lambda_k.begin();
00172         DVectorSlice::iterator  mu_itr = mu.begin();
00173         for( ; lb_itr != lambda_k.end(); ++mu_itr, ++ lb_itr ) {
00174           const value_type lb_j = ::fabs(*lb_itr);
00175           *mu_itr = max(
00176                   (3.0 * (*mu_itr) + lb_j) / 4.0  
00177                 , max( mult_fact * lb_j, small_mu_ )
00178                 );
00179         }
00180         value_type kkt_error = s.opt_kkt_err().get_k(0) + s.feas_kkt_err().get_k(0);
00181         if(kkt_error <= kkt_near_sol_) {
00182           if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00183             out << "\nkkt_error = " << kkt_error << " <= kkt_near_sol = "
00184                 << kkt_near_sol_ << std::endl
00185               << "Switching to forcing mu_k >= mu_km1 in the future\n";
00186           }
00187           near_solution_ = true;
00188         }
00189       }
00190 
00191       // Force the ratio
00192       const value_type
00193           max_mu  = norm_inf( mu() ),
00194           min_mu  = min_mu_ratio_ * max_mu;
00195       for(DVectorSlice::iterator mu_itr = mu.begin(); mu_itr != mu.end(); ++mu_itr)
00196         *mu_itr = max( (*mu_itr), min_mu ); 
00197 
00198       s.mu().set_k(0) = norm_inf_mu_last_ = max_mu;
00199 
00200       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00201         out << "\nmax(|mu(j)|) = " << (*std::max_element( mu.begin(), mu.end() ))
00202           << "\nmin(|mu(j)|) = " << (*std::min_element( mu.begin(), mu.end() ))
00203             << std::endl;
00204       }
00205 
00206       if( (int)olevel >= (int)PRINT_VECTORS ) {
00207         out << "\nmu = \n" << mu;
00208       }
00209     }
00210     else {
00211       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00212         out << "\nDon't have the info to update penalty parameter so just use the last updated...\n";
00213       }
00214     }
00215   }
00216 
00217   // In addition also compute the directional derivative
00218   direc_deriv->calc_deriv( s.Gf().get_k(0)(), s.c().get_k(0)(), s.d().get_k(0)() );
00219 
00220   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00221     out << "\nmu_k = " << s.mu().get_k(0) << "\n";
00222   }
00223 
00224   return true;
00225 }
00226 
00227 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::print_step( const Algorithm& algo
00228   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00229   , std::ostream& out, const std::string& L ) const
00230 {
00231   out
00232     << L << "*** Update the penalty parameter for the merit function to ensure\n"
00233     << L << "*** a descent direction a directional derivatieve.\n"
00234     << L << "*** phi is a merit function object that uses the penalty parameter mu.\n"
00235     << L << "default: near_solution = false\n"
00236     << L << "         small_mu = " << small_mu_ << std::endl
00237     << L << "         min_mu_ratio = " << min_mu_ratio_ << std::endl
00238     << L << "         mult_factor = " << mult_factor_ << std::endl
00239     << L << "         kkt_near_sol = " << kkt_near_sol_ << std::endl
00240     << L << "perform_update = true\n"
00241     << L << "if mu_k is already updated then\n"
00242     << L << "    if mu_k == norm_inf_mu_last then\n"
00243     << L << "        *** We will use this as a signal to skip the update\n"
00244     << L << "        perform_update = false\n"
00245     << L << "    else\n"
00246     << L << "        *** We will perform the update anyway\n"
00247     << L << "    end\n"
00248     << L << "if perform_update == true then\n"
00249     << L << "    if lambda_k is updated then\n"
00250     << L << "        max_lambda = norm(lambda_k,inf)\n"
00251     << L << "        mult_fact = (1+mult_factor)\n"
00252     << L << "        mu = phi.mu()\n"
00253     << L << "        if near_solution == true\n"
00254     << L << "            for j = 1...m\n"
00255     << L << "                mu(j) = max(max(mu(j),mult_fact*abs(lambda_k(j))),small_mu)\n"
00256     << L << "            end\n"
00257     << L << "        else\n"
00258     << L << "            for j = 1...m\n"
00259     << L << "                mu(j) = max( ( 3.0 * mu(j) + abs(lambda_k(j)) ) / 4.0\n"
00260     << L << "                            , max( 1.001 * abs(lambda_k(j)) , small_mu ) )\n"
00261     << L << "            end\n"
00262     << L << "            kkt_error = opt_kkt_err_k + feas_kkt_err_k\n"
00263     << L << "            if kkt_error <= kkt_near_sol then\n"
00264     << L << "                near_solution = true\n"
00265     << L << "            end\n"
00266     << L << "        end\n"
00267     << L << "        min_mu = min_mu_ratio * norm(mu,inf)\n"
00268     << L << "        for j = 1...m\n"
00269     << L << "            mu(j) = max( mu(j), min_mu )\n"
00270     << L << "        end\n"
00271     << L << "    else\n"
00272     << L << "        *** Don't have the information to perform the update.\n"
00273     << L << "    end\n"
00274     << L << "end\n"
00275     << L << "phi.calc_deriv(Gf_k,c_k,d_k)\n";
00276 }
00277 
00278 // Overridden from MeritFunc_PenaltyParamUpdate_AddedStep
00279 
00280 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::small_mu( value_type small_mu )
00281 {
00282   small_mu_ = small_mu;
00283 }
00284 
00285 value_type MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::small_mu() const
00286 {
00287   return small_mu_;
00288 }
00289 
00290 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::min_mu_ratio( value_type min_mu_ratio )
00291 {
00292   min_mu_ratio_ = min_mu_ratio;
00293 }
00294 
00295 value_type MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::min_mu_ratio() const
00296 {
00297   return min_mu_ratio_;
00298 }
00299 
00300 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::mult_factor( value_type mult_factor )
00301 {
00302   mult_factor_ = mult_factor;
00303 }
00304 
00305 value_type MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::mult_factor() const
00306 {
00307   return mult_factor_;
00308 }
00309 
00310 void MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::kkt_near_sol( value_type kkt_near_sol )
00311 {
00312   kkt_near_sol_ = kkt_near_sol;
00313 }
00314 
00315 value_type MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::kkt_near_sol() const
00316 {
00317   return kkt_near_sol_;
00318 }
00319 
00320 
00321 } // end namespace MoochoPack
00322 
00323 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines