MoochoPack_UpdateBarrierParameter_Step.cpp

// @HEADER
// ***********************************************************************
// 
// Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
//                  Copyright (2003) Sandia Corporation
// 
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// 
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//  
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
// 
// ***********************************************************************
// @HEADER

#include <ostream>
#include <typeinfo>
#include <iostream>
#include <math.h>

#include "MoochoPack_UpdateBarrierParameter_Step.hpp"
#include "MoochoPack_moocho_algo_conversion.hpp"
#include "IterationPack_print_algorithm_step.hpp"
#include "Teuchos_dyn_cast.hpp"
#include "MoochoPack_IpState.hpp"
#include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
#include "AbstractLinAlgPack_VectorStdOps.hpp"
#include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
#include "Teuchos_TestForException.hpp"

#include "OptionsFromStreamPack_StringToBool.hpp"

#define min(a,b) ( (a < b) ? a : b )
#define max(a,b) ( (a > b) ? a : b )

namespace MoochoPack {

UpdateBarrierParameter_Step::UpdateBarrierParameter_Step(
  const value_type init_barrier_parameter,
  const value_type tau_mu,
  const value_type theta_mu,
  const value_type tau_epsilon,
  const value_type theta_epsilon,
  const value_type e_tol_max
  )
  :
  init_barrier_parameter_(init_barrier_parameter),
  tau_mu_(tau_mu),
  theta_mu_(theta_mu),
  tau_epsilon_(tau_epsilon),
  theta_epsilon_(theta_epsilon),
  e_tol_max_(e_tol_max)
  {}
  

bool UpdateBarrierParameter_Step::do_step(
  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
  ,poss_type assoc_step_poss
  )
  {
  using Teuchos::dyn_cast;
  using IterationPack::print_algorithm_step;

  NLPAlgo            &algo   = dyn_cast<NLPAlgo>(_algo);
  IpState             &s      = dyn_cast<IpState>(_algo.state());
  NLP                 &nlp    = algo.nlp();

  if (!nlp.is_initialized())
    {
    nlp.initialize(false);
    }
  
  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
  std::ostream& out = algo.track().journal_out();
  
  // print step header.
  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
    {
    using IterationPack::print_algorithm_step;
    print_algorithm_step( _algo, step_poss, type, assoc_step_poss, out );
    }

  
  // Get iteration quantities
  IterQuantityAccess<value_type>& e_tol_iq = s.e_tol();
  IterQuantityAccess<value_type>& mu_iq = s.barrier_parameter();

  // Check values and initialize, if necessary
  /*  if (mu_iq.last_updated() == IterQuantity::NONE_UPDATED)
    {
    mu_iq.set_k(0) = init_barrier_parameter_;
    e_tol_iq.set_k(0) = Calculate_e_tol(mu_iq.get_k(0));

    if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
      {
      out << "\nInitializing barrier parameter (mu) and sub problem tolerance (e_tol) ...\n";
      }
    }
    else*/
    {
    // if tolerance is satisfied, calculate new barrier_parameter 
    //  and e_tol, otherwise update mu and e_tol from last iter
    const value_type opt_err = s.opt_kkt_err().get_k(0);
    const value_type feas_err = s.feas_kkt_err().get_k(0);
    const value_type comp_err_mu = s.comp_err_mu().get_k(0);

    const value_type mu_km1 = mu_iq.get_k(-1);
    if (e_tol_iq.last_updated() == IterQuantity::NONE_UPDATED)
      {
      // First time through, don't let mu change
      mu_iq.set_k(0,-1);
      e_tol_iq.set_k(0) = Calculate_e_tol(mu_iq.get_k(0));
          }
    else
      {
      const value_type e_tol_km1 = e_tol_iq.get_k(-1);
      bool sub_prob_converged = (opt_err < e_tol_km1 && feas_err < e_tol_km1 && comp_err_mu < e_tol_km1); 
      if (sub_prob_converged)
        {
        // Calculate new mu and e_tol
        value_type& mu_k = mu_iq.set_k(0);
        mu_k = min(tau_mu_*mu_km1, pow(mu_km1, theta_mu_));
        e_tol_iq.set_k(0) = Calculate_e_tol(mu_k);
        
        if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
          {
          out << "\nSub-problem converged!\n"
            << " Updating barrier parameter (mu) and sub problem tolerance (e_tol) ...\n";
          }
        }
      else
        {
        if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
          {
          out << "\nSub-problem not-converged!\n"
            << " Keeping existing barrier parameter (mu) and sub problem tolerance (e_tol) ...\n";
          }
        mu_iq.set_k(0,-1);
        e_tol_iq.set_k(0,-1);
        }
      }

    if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
      {
      out << "\nbarrier_parameter (mu) = " << mu_iq.get_k(0)
        << "\nsub problem tolerance (e_tol) = " << e_tol_iq.get_k(0)  << std::endl;
      }
    }

  return true;
  }


void UpdateBarrierParameter_Step::print_step(
  const Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
  ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
  ) const
  {
  //const NLPAlgo   &algo = rsqp_algo(_algo);
  //const NLPAlgoState  &s    = algo.rsqp_state();
  out << L << "# Update the interior point barrier parameter (mu)\n"
    << L << "if (KKTerror < e_tol) then\n"
    << L << "   mu_kp1 = min(tau_mu*mu_k, mu_k^theta_mu)\n"
    << L << "   e_tol_kp1 = min(e_tol_max, tau_epsilon*min(mu_k, mu_k^theta_epsilon))\n"
    << L << "else\n"
    << L << "   mu_kp1 = mu_k\n"
    << L << "   e_tol_kp1 = e_tol_k\n"
    << L << "end;\n";
  }

value_type UpdateBarrierParameter_Step::Calculate_e_tol(value_type mu)
  { 
  value_type e_tol = tau_epsilon_*min(mu, pow(mu, theta_epsilon_));
  e_tol = min(e_tol_max_, e_tol);

  return e_tol;
  }


namespace {

const int local_num_options = 5;

enum local_EOptions 
  {
    TAU_MU,
    THETA_MU,
    TAU_EPSILON,
    THETA_EPSILON,
    E_TOL_MAX
  };

const char* local_SOptions[local_num_options] = 
  {
    "tau_mu",
    "theta_mu",
    "tau_epsilon",
    "theta_epsilon",
    "e_tol_max"
  };

}

 
UpdateBarrierParameter_StepSetOptions::UpdateBarrierParameter_StepSetOptions(
  UpdateBarrierParameter_Step* target
  , const char opt_grp_name[] )
  :
  OptionsFromStreamPack::SetOptionsFromStreamNode(
    opt_grp_name, local_num_options, local_SOptions ),
  OptionsFromStreamPack::SetOptionsToTargetBase< UpdateBarrierParameter_Step >( target )
  {
  }

void UpdateBarrierParameter_StepSetOptions::setOption( 
  int option_num, const std::string& option_value )
  {
  using OptionsFromStreamPack::StringToBool;
  
  typedef UpdateBarrierParameter_Step target_t;
  switch( (local_EOptions)option_num ) 
    {
    case TAU_MU:
      target().tau_mu(std::atof(option_value.c_str()));
      break;
    case THETA_MU:
      target().theta_mu(std::atof(option_value.c_str()));
      break;
    case TAU_EPSILON:
      target().tau_epsilon(std::atof(option_value.c_str()));
      break;
    case THETA_EPSILON:
      target().theta_epsilon(std::atof(option_value.c_str()));
      break;
    case E_TOL_MAX:
      target().e_tol_max(std::atof(option_value.c_str()));
      break;
    default:
      TEST_FOR_EXCEPT(true);  // Local error only?
    }
  }

} // end namespace MoochoPack