Changeset 542


Ignore:
Timestamp:
Oct 13, 2005 6:43:08 PM (14 years ago)
Author:
andreasw
Message:
  • cleaned up line search to allow for alternative globalization scheme
Location:
branches/dev/Algorithm
Files:
1 added
8 edited
2 copied
2 moved

Legend:

Unmodified
Added
Removed
  • branches/dev/Algorithm/IpAdaptiveMuUpdate.cpp

    r529 r542  
    374374
    375375      linesearch_->Reset();
    376       // Uncomment the next line if the filter should not switch to
     376      // Uncomment the next line if the line search should not switch to
    377377      // the restoration phase in the free mode
     378
    378379      // linesearch_->SetRigorousLineSearch(false);
    379380    }
  • branches/dev/Algorithm/IpAlgBuilder.cpp

    r517 r542  
    1414#include "IpPDPerturbationHandler.hpp"
    1515#include "IpOptErrorConvCheck.hpp"
    16 #include "IpFilterLineSearch.hpp"
     16#include "IpBacktrackingLineSearch.hpp"
     17#include "IpFilterLSAcceptor.hpp"
    1718#include "IpMonotoneMuUpdate.hpp"
    1819#include "IpAdaptiveMuUpdate.hpp"
     
    207208    SmartPtr<RestoRestorationPhase> resto_resto =
    208209      new RestoRestorationPhase();
    209     SmartPtr<FilterLineSearch> resto_LineSearch =
    210       new FilterLineSearch(GetRawPtr(resto_resto), GetRawPtr(resto_PDSolver),
    211                            GetRawPtr(resto_convCheck));
     210    SmartPtr<FilterLSAcceptor> resto_filterLSacceptor =
     211      new FilterLSAcceptor(GetRawPtr(resto_PDSolver));
     212    SmartPtr<LineSearch> resto_LineSearch =
     213      new BacktrackingLineSearch(GetRawPtr(resto_filterLSacceptor),
     214                                 GetRawPtr(resto_resto), GetRawPtr(resto_convCheck));
    212215
    213216    // Create the mu update that will be used by the restoration phase
     
    282285
    283286    // Create the line search to be used by the main algorithm
    284     SmartPtr<FilterLineSearch> lineSearch =
    285       new FilterLineSearch(GetRawPtr(resto_phase), GetRawPtr(PDSolver),
    286                            convCheck);
     287    SmartPtr<FilterLSAcceptor> filterLSacceptor =
     288      new FilterLSAcceptor(GetRawPtr(PDSolver));
     289    SmartPtr<LineSearch> lineSearch =
     290      new BacktrackingLineSearch(GetRawPtr(filterLSacceptor),
     291                                 GetRawPtr(resto_phase), convCheck);
    287292
    288293    // The following cross reference is not good: We have to store a
     
    290295    // non-SmartPtr to make sure that things are properly deleted when
    291296    // the IpoptAlgorithm return by the Builder is destructed.
    292     resto_convCheck->SetOrigFilterLineSearch(*lineSearch);
     297    resto_convCheck->SetOrigFilterLSAcceptor(*filterLSacceptor);
    293298
    294299    // Create the mu update that will be used by the main algorithm
  • branches/dev/Algorithm/IpAlgorithmRegOp.cpp

    r507 r542  
    1313#include "IpAlgBuilder.hpp"
    1414#include "IpDefaultIterateInitializer.hpp"
    15 #include "IpFilterLineSearch.hpp"
     15#include "IpBacktrackingLineSearch.hpp"
     16#include "IpFilterLSAcceptor.hpp"
    1617#include "IpGradientScaling.hpp"
    1718#include "IpIpoptAlg.hpp"
     
    4546    AlgorithmBuilder::RegisterOptions(roptions);
    4647    roptions->SetRegisteringCategory("Line Search");
    47     FilterLineSearch::RegisterOptions(roptions);
     48    BacktrackingLineSearch::RegisterOptions(roptions);
     49    FilterLSAcceptor::RegisterOptions(roptions);
    4850    roptions->SetRegisteringCategory("NLP Scaling");
    4951    StandardScalingBase::RegisterOptions(roptions);
  • branches/dev/Algorithm/IpBacktrackingLineSearch.cpp

    r541 r542  
    66//
    77// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
    8 
    9 #include "IpFilterLineSearch.hpp"
     8//           Andreas Waechter                 IBM    2005-10-13
     9//               derived file from IpFilterLineSearch.cpp
     10
     11#include "IpBacktrackingLineSearch.hpp"
    1012#include "IpJournalist.hpp"
    1113#include "IpRestoPhase.hpp"
     
    3941#endif
    4042
    41   FilterLineSearch::FilterLineSearch(const SmartPtr<RestorationPhase>& resto_phase,
    42                                      const SmartPtr<PDSystemSolver>& pd_solver,
    43                                      const SmartPtr<ConvergenceCheck>& conv_check)
     43  BacktrackingLineSearch::BacktrackingLineSearch(
     44    const SmartPtr<BacktrackingLSAcceptor>& acceptor,
     45    const SmartPtr<RestorationPhase>& resto_phase,
     46    const SmartPtr<ConvergenceCheck>& conv_check)
    4447      :
    4548      LineSearch(),
    46       filter_(2),
     49      acceptor_(acceptor),
    4750      resto_phase_(resto_phase),
    48       pd_solver_(pd_solver),
    4951      conv_check_(conv_check)
    5052  {
    51     DBG_START_FUN("FilterLineSearch::FilterLineSearch",
     53    DBG_START_FUN("BacktrackingLineSearch::BacktrackingLineSearch",
    5254                  dbg_verbosity);
    53   }
    54 
    55   FilterLineSearch::~FilterLineSearch()
    56   {
    57     DBG_START_FUN("FilterLineSearch::~FilterLineSearch()",
     55    DBG_ASSERT(IsValid(acceptor_));
     56  }
     57
     58  BacktrackingLineSearch::~BacktrackingLineSearch()
     59  {
     60    DBG_START_FUN("BacktrackingLineSearch::~BacktrackingLineSearch()",
    5861                  dbg_verbosity);
    5962  }
    6063
    61   void FilterLineSearch::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
    62   {
    63     roptions->AddLowerBoundedNumberOption(
    64       "theta_max_fact",
    65       "Determines upper bound for constraint violation in the filter.",
    66       0.0, true, 1e4,
    67       "The algorithmic parameter theta_max is determined as theta_max_fact "
    68       "times the maximum of 1 and the constraint violation at initial point.  "
    69       "Any point with a constraint violation larger than theta_max is "
    70       "unacceptable to the filter (see Eqn. (21) in implementation paper).");
    71     roptions->AddLowerBoundedNumberOption(
    72       "theta_min_fact",
    73       "Determines constraint violation threshold in the switching rule.",
    74       0.0, true, 1e-4,
    75       "The algorithmic parameter theta_min is determined as theta_min_fact "
    76       "times the maximum of 1 and the constraint violation at initial point.  "
    77       "The switching rules treats an iteration as an h-type iteration whenever "
    78       "the current constraint violation is larger than theta_min (see "
    79       "paragraph before Eqn. (19) in implementation paper).");
    80     roptions->AddBoundedNumberOption(
    81       "eta_phi",
    82       "Relaxation factor in the Armijo condition.",
    83       0.0, true, 0.5, true, 1e-8,
    84       "(See Eqn. (20) in implementation paper)");
    85     roptions->AddLowerBoundedNumberOption(
    86       "delta", "Multiplier for constraint violation in the switching rule.",
    87       0.0, true, 1.0,
    88       "(See Eqn. (19) in implementation paper.)");
    89     roptions->AddLowerBoundedNumberOption(
    90       "s_phi",
    91       "Exponent for linear barrier function model in the switching rule.",
    92       1.0, true, 2.3,
    93       "(See Eqn. (19) in implementation paper.)");
    94     roptions->AddLowerBoundedNumberOption(
    95       "s_theta",
    96       "Exponent for current constraint violation in the switching rule.",
    97       1.0, true, 1.1,
    98       "(See Eqn. (19) in implementation paper.)");
    99     roptions->AddBoundedNumberOption(
    100       "gamma_phi",
    101       "Relaxation factor in the filter margin for the barrier function.",
    102       0.0, true, 1.0, true, 1e-8,
    103       "(See Eqn. (18a) in implementation paper.)");
    104     roptions->AddBoundedNumberOption(
    105       "gamma_theta",
    106       "Relaxation factor in the filter margin for the constraint violation.",
    107       0.0, true, 1.0, true, 1e-5,
    108       "(See Eqn. (18b) in implementation paper.)");
    109     roptions->AddBoundedNumberOption(
    110       "alpha_min_frac",
    111       "Safety factor for the minimal step size (before switching to restoration phase).",
    112       0.0, true, 1.0, true, 0.05,
    113       "(This is gamma_alpha in Eqn. (20) in implementation paper.)");
     64  void BacktrackingLineSearch::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
     65  {
    11466    roptions->AddBoundedNumberOption(
    11567      "alpha_red_factor",
     
    11870      "At every step of the backtracking line search, the trial step size is "
    11971      "reduced by this factor.");
    120     roptions->AddLowerBoundedIntegerOption(
    121       "max_soc",
    122       "Maximum number of second order correction trial steps at each iteration.",
    123       0, 4,
    124       "Choosing 0 disables the second order "
    125       "corrections. (This is p^{max} of Step A-5.9 of "
    126       "Algorithm A in implementation paper.)");
    127     roptions->AddLowerBoundedNumberOption(
    128       "kappa_soc",
    129       "Factor in the sufficient reduction rule for second order correction.",
    130       0.0, true, 0.99,
    131       "This option determines how much a second order correction step must reduce the "
    132       "constraint violation so that further correction steps are attempted.  "
    133       "(See Step A-5.9 of Algorithm A in implementation paper.)");
    134     roptions->AddLowerBoundedNumberOption(
    135       "obj_max_inc",
    136       "Determines the upper bound on the acceptable increase of barrier objective function.",
    137       1.0, true, 5.0,
    138       "Trial points are rejected if they lead to an increase in the "
    139       "barrier objective function by more than obj_max_inc orders "
    140       "of magnitude.");
    14172
    14273    std::string prev_category = roptions->RegisteringCategory();
     
    15081      "DOESN'T REALLY WORK YET!");
    15182    roptions->SetRegisteringCategory(prev_category);
    152     roptions->AddStringOption3(
    153       "corrector_type",
    154       "The type of corrector steps that should be taken.",
    155       "none",
    156       "none", "no corrector",
    157       "affine", "corrector step towards mu=0",
    158       "primal-dual", "corrector step towards current mu",
    159       "If \"mu_strategy\" is \"adaptive\", this option determines "
    160       "what kind of corrector steps should be tried.");
    161 
    162     roptions->AddStringOption2(
    163       "skip_corr_if_neg_curv",
    164       "Skip the corrector step in negative curvature iteration.",
    165       "yes",
    166       "no", "don't skip",
    167       "yes", "skip",
    168       "The corrector step is not tried if negative curvature has been "
    169       "encountered during the computation of the search direction in "
    170       "the current iteration. This option is only used if \"mu_strategy\" is "
    171       "\"adaptive\".");
    172 
    173     roptions->AddStringOption2(
    174       "skip_corr_in_monotone_mode",
    175       "Skip the corrector step during monotone barrier parameter mode.",
    176       "yes",
    177       "no", "don't skip",
    178       "yes", "skip",
    179       "The corrector step is not tried if the algorithm is currently in the "
    180       "monotone mode (see also option \"barrier_strategy\")."
    181       "This option is only used if \"mu_strategy\" is \"adaptive\".");
    18283
    18384    roptions->AddStringOption2(
     
    203104      "This option determines how the step size (alpha_y) will be calculated when updating the "
    204105      "constraint multipliers.");
    205 
    206     roptions->AddLowerBoundedNumberOption(
    207       "corrector_compl_avrg_red_fact",
    208       "Complementarity tolerance factor for accepting corrector step.",
    209       0.0, true, 1.0,
    210       "This option determines the factor by which complementarity is allowed to increase "
    211       "for a corrector step to be accepted.");
    212106
    213107    roptions->AddStringOption2(
     
    231125      "search are disabled. If the problem is square, this options is set to "
    232126      "0.");
    233     roptions->AddLowerBoundedNumberOption(
    234       "soft_resto_pderror_reduction_factor",
    235       "Required reduction in primal-dual error in the soft restoration phase.",
    236       0.0, false, (1.0 - 1e-4),
    237       "The soft restoration phase attempts to reduce the "
    238       "primal-dual error with regular steps. If the damped "
    239       "primal-dual step (damped only to satisfy the "
    240       "fraction-to-the-boundary rule) is not decreasing the primal-dual error "
    241       "by at least this factor, then the regular restoration phase is called. "
    242       "Choosing \"0\" here disables the soft "
    243       "restoration phase.");
    244127    roptions->AddStringOption2(
    245128      "start_with_resto",
     
    275158      "allowed before the watchdog "
    276159      "procedure is aborted and the algorithm returns to the stored point.");
    277   }
    278 
    279   bool FilterLineSearch::InitializeImpl(const OptionsList& options,
    280                                         const std::string& prefix)
    281   {
    282     options.GetNumericValue("theta_max_fact", theta_max_fact_, prefix);
    283     options.GetNumericValue("theta_min_fact", theta_min_fact_, prefix);
    284     ASSERT_EXCEPTION(theta_min_fact_ < theta_max_fact_, OPTION_INVALID,
    285                      "Option \"theta_min_fact\": This value must be larger than 0 and less than theta_max_fact.");
    286     options.GetNumericValue("eta_phi", eta_phi_, prefix);
    287     options.GetNumericValue("delta", delta_, prefix);
    288     options.GetNumericValue("s_phi", s_phi_, prefix);
    289     options.GetNumericValue("s_theta", s_theta_, prefix);
    290     options.GetNumericValue("gamma_phi", gamma_phi_, prefix);
    291     options.GetNumericValue("gamma_theta", gamma_theta_, prefix);
    292     options.GetNumericValue("alpha_min_frac", alpha_min_frac_, prefix);
     160
     161    roptions->AddLowerBoundedNumberOption(
     162      "soft_resto_pderror_reduction_factor",
     163      "Required reduction in primal-dual error in the soft restoration phase.",
     164      0.0, false, (1.0 - 1e-4),
     165      "The soft restoration phase attempts to reduce the "
     166      "primal-dual error with regular steps. If the damped "
     167      "primal-dual step (damped only to satisfy the "
     168      "fraction-to-the-boundary rule) is not decreasing the primal-dual error "
     169      "by at least this factor, then the regular restoration phase is called. "
     170      "Choosing \"0\" here disables the soft "
     171      "restoration phase.");
     172  }
     173
     174  bool BacktrackingLineSearch::InitializeImpl(const OptionsList& options,
     175      const std::string& prefix)
     176  {
    293177    options.GetNumericValue("alpha_red_factor", alpha_red_factor_, prefix);
    294     options.GetIntegerValue("max_soc", max_soc_, prefix);
    295     if (max_soc_>0) {
    296       ASSERT_EXCEPTION(IsValid(pd_solver_), OPTION_INVALID,
    297                        "Option \"max_soc\": This option is non-negative, but no linear solver for computing the SOC given to FilterLineSearch object.");
    298     }
    299     options.GetNumericValue("kappa_soc", kappa_soc_, prefix);
    300     options.GetNumericValue("obj_max_inc", obj_max_inc_, prefix);
    301178    options.GetBoolValue("magic_steps", magic_steps_, prefix);
     179    options.GetBoolValue("accept_every_trial_step", accept_every_trial_step_, prefix);
    302180    Index enum_int;
    303     options.GetEnumValue("corrector_type", enum_int, prefix);
    304     corrector_type_ = CorrectorTypeEnum(enum_int);
    305     options.GetBoolValue("skip_corr_if_neg_curv", skip_corr_if_neg_curv_, prefix);
    306     options.GetBoolValue("skip_corr_in_monotone_mode", skip_corr_in_monotone_mode_, prefix);
    307     options.GetBoolValue("accept_every_trial_step", accept_every_trial_step_, prefix);
    308181    options.GetEnumValue("alpha_for_y", enum_int, prefix);
    309182    alpha_for_y_ = AlphaForYEnum(enum_int);
    310     options.GetNumericValue("corrector_compl_avrg_red_fact", corrector_compl_avrg_red_fact_, prefix);
    311183    options.GetNumericValue("expect_infeasible_problem_ctol", expect_infeasible_problem_ctol_, prefix);
    312184    options.GetBoolValue("expect_infeasible_problem", expect_infeasible_problem_, prefix);
     
    314186    options.GetBoolValue("start_with_resto", start_with_resto_, prefix);
    315187
    316     bool retvalue = true;
    317     if (IsValid(resto_phase_)) {
    318       retvalue = resto_phase_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
    319                                           options, prefix);
    320     }
    321 
    322     options.GetNumericValue("soft_resto_pderror_reduction_factor",
    323                             soft_resto_pderror_reduction_factor_, prefix);
    324188    options.GetNumericValue("tiny_step_tol", tiny_step_tol_, prefix);
    325189    options.GetIntegerValue("watchdog_trial_iter_max", watchdog_trial_iter_max_, prefix);
    326190    options.GetIntegerValue("watchdog_shortened_iter_trigger", watchdog_shortened_iter_trigger_, prefix);
     191    options.GetNumericValue("soft_resto_pderror_reduction_factor",
     192                            soft_resto_pderror_reduction_factor_, prefix);
     193
     194    bool retvalue = true;
     195    if (IsValid(resto_phase_)) {
     196      if (!resto_phase_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
     197                                    options, prefix)) {
     198        return false;
     199      }
     200      if (!acceptor_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
     201                                 options, prefix)) {
     202        return false;
     203      }
     204    }
    327205
    328206    rigorous_ = true;
    329207    skipped_line_search_ = false;
    330208    tiny_step_last_iteration_ = false;
    331     theta_min_ = -1.;
    332     theta_max_ = -1.;
    333209
    334210    Reset();
     
    341217  }
    342218
    343   void FilterLineSearch::FindAcceptableTrialPoint()
    344   {
    345     DBG_START_METH("FilterLineSearch::FindAcceptableTrialPoint",
     219  void BacktrackingLineSearch::FindAcceptableTrialPoint()
     220  {
     221    DBG_START_METH("BacktrackingLineSearch::FindAcceptableTrialPoint",
    346222                   dbg_verbosity);
    347223    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
     
    368244    skipped_line_search_ = false;
    369245
    370     // Set the values for the reference point
    371     if (!in_watchdog_) {
    372       reference_theta_ = IpCq().curr_constraint_violation();
    373       reference_barr_ = IpCq().curr_barrier_obj();
    374       reference_gradBarrTDelta_ = IpCq().curr_gradBarrTDelta();
    375     }
    376     else {
    377       reference_theta_ = watchdog_theta_;
    378       reference_barr_ = watchdog_barr_;
    379       reference_gradBarrTDelta_ = watchdog_gradBarrTDelta_;
    380     }
     246    // Initialize the acceptor for this backtracking line search
     247    acceptor_->InitThisLineSearch(in_watchdog_);
    381248
    382249    // Get the search directions (this will store the actual search
    383250    // direction, possibly including higher order corrections)
    384     SmartPtr<IteratesVector> actual_delta = IpData().delta()->MakeNewContainer();
     251    SmartPtr<IteratesVector> actual_delta =
     252      IpData().delta()->MakeNewContainer();
    385253
    386254    bool goto_resto = false;
     
    393261
    394262    if (start_with_resto_) {
    395       // If the use requested to start with the restoration phase,
    396       // skip the lin e search and do exactly that.  Reset the flag so
     263      // If the user requested to start with the restoration phase,
     264      // skip the line search and do exactly that.  Reset the flag so
    397265      // that this happens only once.
    398266      goto_resto = true;
     
    407275
    408276    // Check if search direction becomes too small
    409     // ToDo: move this into place independent of this particular line search?
    410277    bool tiny_step = (!goto_resto && DetectTinyStep());
    411278
     
    427294    }
    428295
     296    // Handle the situation of a tiny step
    429297    if (tiny_step) {
    430298      alpha_primal =
     
    454322
    455323      if (in_soft_resto_phase_) {
    456         bool satisfies_original_filter = false;
     324        // If we are currently in the soft restoration phase, continue
     325        // that way, and switch back if enough progress is made to the
     326        // original criterion (e.g., the filter)
     327        bool satisfies_original_criterion = false;
    457328        // ToDo use tiny_step in TrySoftRestoStep?
    458329        accept = TrySoftRestoStep(actual_delta,
    459                                   satisfies_original_filter);
     330                                  satisfies_original_criterion);
    460331        if (accept) {
    461332          IpData().Set_info_alpha_primal_char('s');
    462           if (satisfies_original_filter) {
     333          if (satisfies_original_criterion) {
    463334            in_soft_resto_phase_ = false;
    464335            IpData().Set_info_alpha_primal_char('S');
     
    467338      }
    468339      else {
     340        // Start the backtracking line search
    469341        bool done = false;
    470342        bool skip_first_trial_point = false;
     
    507379    } /* if (!goto_resto && !tiny_step) { */
    508380
    509     // If line search has been aborted because the step size becomes too small,
    510     // go to the restoration phase
     381    // If line search has been aborted because the step size becomes
     382    // too small, go to the restoration phase or continue with soft
     383    // restoration phase
    511384    if (!accept) {
    512385      // If we are not asked to do a rigorous line search, do no call
     
    518391      else {
    519392        // Check if we should start the soft restoration phase
    520         if (!in_soft_resto_phase_ && soft_resto_pderror_reduction_factor_>0.
     393        if (!in_soft_resto_phase_
    521394            && !goto_resto && !expect_infeasible_problem_) {
    522395          Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    523396                         "--> Starting soft restoration phase <--\n");
    524           // Augment the filter with the current point
    525           AugmentFilter();
     397          // Prepare the restoration phase, e.g., augment the filter
     398          // with the current point.
     399          acceptor_->PrepareRestoPhaseStart();
    526400
    527401          // Try the current search direction for the soft restoration phase
    528           bool satisfies_original_filter;
     402          bool satisfies_original_criterion;
    529403          accept = TrySoftRestoStep(actual_delta,
    530                                     satisfies_original_filter);
    531           // If it has been accepted: If the original filter is also
     404                                    satisfies_original_criterion);
     405          // If it has been accepted: If the original criterion is also
    532406          // satisfied, we can just take that step and continue with
    533407          // the regular algorithm, otherwise we stay in the soft
    534408          // restoration phase
    535409          if (accept) {
    536             if (satisfies_original_filter) {
     410            if (satisfies_original_criterion) {
    537411              IpData().Set_info_alpha_primal_char('S');
    538412            }
     
    545419
    546420        if (!accept) {
     421          // Go to the restoration phase
    547422          if (!in_soft_resto_phase_) {
    548             // Augment the filter with the current point if we are
    549             // already in the soft restoration phase, this has been
    550             // done earlier
    551             AugmentFilter();
     423            // Prepare the restoration phase, e.g., augment the filter
     424            // with the current point. If we are already in the soft
     425            // restoration phase, this has been done earlier
     426            acceptor_->PrepareRestoPhaseStart();
    552427          }
    553428          if (CurrentIsAcceptable()) {
     
    633508  }
    634509
    635   bool FilterLineSearch::DoBacktrackingLineSearch(bool skip_first_trial_point,
     510  bool BacktrackingLineSearch::DoBacktrackingLineSearch(bool skip_first_trial_point,
    636511      Number& alpha_primal,
    637512      bool& corr_taken,
     
    644519    bool accept = false;
    645520
    646     DBG_START_METH("FilterLineSearch::DoBacktrackingLineSearch",
     521    DBG_START_METH("BacktrackingLineSearch::DoBacktrackingLineSearch",
    647522                   dbg_verbosity);
    648523
     
    656531    Number alpha_min = alpha_primal_max;
    657532    if (!in_watchdog_) {
    658       alpha_min = CalculateAlphaMin();
     533      alpha_min = acceptor_->CalculateAlphaMin();
    659534    }
    660535    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
     
    674549    }
    675550
    676     filter_.Print(Jnlst());
    677 
    678     if (corrector_type_!=NO_CORRECTOR && !skip_first_trial_point &&
    679         (!skip_corr_if_neg_curv_ || IpData().info_regu_x()==0.) &&
    680         (!skip_corr_in_monotone_mode_ || IpData().FreeMuMode()) ) {
     551    if (!skip_first_trial_point) {
    681552      // Before we do the actual backtracking line search for the
    682553      // regular primal-dual search direction, let's see if a step
    683554      // including a higher-order correctior is already acceptable
    684       accept = TryCorrector(alpha_primal_test,
    685                             alpha_primal,
    686                             actual_delta);
     555      accept = acceptor_->TryCorrector(alpha_primal_test,
     556                                       alpha_primal,
     557                                       actual_delta);
    687558    }
    688559    if (accept) {
     
    711582          // If it is acceptable, stop the search
    712583          alpha_primal_test = alpha_primal;
    713           accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test);
     584          if (accept_every_trial_step_) {
     585            // We call the evaluation at the trial point here, so that an
     586            // exception will the thrown if there are problem during the
     587            // evaluation of the functions (in that case, we want to further
     588            // reduce the step size
     589            IpCq().trial_barrier_obj();
     590            IpCq().trial_constraint_violation();
     591            accept = true;
     592          }
     593          else {
     594            accept = acceptor_->CheckAcceptabilityOfTrialPoint(alpha_primal_test);
     595          }
    714596        }
    715597        catch(IpoptNLP::Eval_Error& e) {
     
    745627          Number theta_trial = IpCq().trial_constraint_violation();
    746628          if (alpha_primal==alpha_primal_max &&       // i.e. first trial point
    747               theta_curr<=theta_trial && max_soc_>0) {
     629              theta_curr<=theta_trial) {
    748630            // Try second order correction
    749             accept = TrySecondOrderCorrection(alpha_primal_test,
    750                                               alpha_primal,
    751                                               actual_delta);
     631            accept = acceptor_->TrySecondOrderCorrection(alpha_primal_test,
     632                     alpha_primal,
     633                     actual_delta);
    752634          }
    753635          if (accept) {
     
    767649      info_alpha_primal_char = 'w';
    768650    }
    769     else {
    770       // Augment the filter if required
    771       if (!IsFtype(alpha_primal_test) ||
    772           !ArmijoHolds(alpha_primal_test)) {
    773         AugmentFilter();
    774         info_alpha_primal_char = 'h';
    775       }
    776       else {
    777         info_alpha_primal_char = 'f';
    778       }
     651    else if (accept) {
     652      info_alpha_primal_char =
     653        acceptor_->UpdateForNextIteration(alpha_primal_test);
    779654    }
    780655    if (soc_taken) {
     
    790665  }
    791666
    792   bool FilterLineSearch::IsFtype(Number alpha_primal_test)
    793   {
    794     DBG_START_METH("FilterLineSearch::IsFtype",
    795                    dbg_verbosity);
    796     DBG_ASSERT(reference_theta_>0. || reference_gradBarrTDelta_ < 0.0);
    797     return (reference_gradBarrTDelta_ < 0.0 &&
    798             alpha_primal_test*pow(-reference_gradBarrTDelta_,s_phi_) >
    799             delta_*pow(reference_theta_,s_theta_));
    800   }
    801 
    802   void FilterLineSearch::AugmentFilter()
    803   {
    804     DBG_START_METH("FilterLineSearch::AugmentFilter",
    805                    dbg_verbosity);
    806 
    807     Number phi_add = reference_barr_ - gamma_phi_*reference_theta_;
    808     Number theta_add = (1.-gamma_theta_)*reference_theta_;
    809 
    810     filter_.AddEntry(phi_add, theta_add, IpData().iter_count());
    811   }
    812 
    813   bool
    814   FilterLineSearch::CheckAcceptabilityOfTrialPoint(Number alpha_primal_test)
    815   {
    816     DBG_START_METH("FilterLineSearch::CheckAcceptabilityOfTrialPoint",
    817                    dbg_verbosity);
    818 
    819     if (accept_every_trial_step_) {
    820 
    821       // We call the evaluation at the trial point here, so that an
    822       // exception will the thrown if there are problem during the
    823       // evaluation of the functions (in that case, we want to further
    824       // reduce the step size
    825       /* Number trial_barr = */ IpCq().trial_barrier_obj();
    826       /* Number trial_theta = */
    827       IpCq().trial_constraint_violation();
    828       return true;
    829     }
    830 
    831     bool accept;
    832 
    833     // First compute the barrier function and constraint violation at the
    834     // current iterate and the trial point
    835 
    836     Number trial_theta = IpCq().trial_constraint_violation();
    837     // Check if constraint violation is becoming too large
    838     if (theta_max_ < 0.0) {
    839       theta_max_ = theta_max_fact_*Max(1.0, reference_theta_);
    840     }
    841     if (theta_min_ < 0.0) {
    842       theta_min_ = theta_min_fact_*Max(1.0, reference_theta_);
    843     }
    844 
    845     if (theta_max_>0 && trial_theta>theta_max_) {
    846       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    847                      "trial_theta = %e is larger than theta_max = %e\n",
    848                      trial_theta, theta_max_);
    849       IpData().Append_info_string("Tmax");
    850       return false;
    851     }
    852 
    853     Number trial_barr = IpCq().trial_barrier_obj();
    854     DBG_ASSERT(IsFiniteNumber(trial_barr));
    855 
    856     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    857                    "Checking acceptability for trial step size alpha_primal_test=%13.6e:\n", alpha_primal_test);
    858     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    859                    "  New values of barrier function     = %23.16e  (reference %23.16e):\n", trial_barr, reference_barr_);
    860     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    861                    "  New values of constraint violation = %23.16e  (reference %23.16e):\n", trial_theta, reference_theta_);
    862 
    863     // Check if point is acceptable w.r.t current iterate
    864     if (IsFtype(alpha_primal_test) && reference_theta_ <= theta_min_) {
    865       // Armijo condition for the barrier function has to be satisfied
    866       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Checking Armijo Condition...\n");
    867       accept = ArmijoHolds(alpha_primal_test);
    868     }
    869     else {
    870       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Checking sufficient reduction...\n");
    871       accept = IsAcceptableToCurrentIterate(trial_barr, trial_theta);
    872     }
    873 
    874     if (!accept) {
    875       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Failed...\n");
    876       return accept;
    877     }
    878     else {
    879       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Succeeded...\n");
    880     }
    881 
    882     // Now check if that pair is acceptable to the filter
    883     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Checking filter acceptability...\n");
    884     accept = IsAcceptableToCurrentFilter(trial_barr, trial_theta);
    885     if (!accept) {
    886       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Failed...\n");
    887       return accept;
    888     }
    889     else {
    890       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Succeeded...\n");
    891     }
    892 
    893     return accept;
    894   }
    895 
    896   bool FilterLineSearch::ArmijoHolds(Number alpha_primal_test)
    897   {
    898     /*
    899     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    900                    "ArmijoHolds test with trial_barr = %25.16e reference_barr = %25.16e\n        alpha_primal_test = %25.16e reference_gradBarrTDelta = %25.16e\n", IpCq().trial_barrier_obj(), reference_barr_,alpha_primal_test,reference_gradBarrTDelta_);
    901     */
    902     return Compare_le(IpCq().trial_barrier_obj()-reference_barr_,
    903                       eta_phi_*alpha_primal_test*reference_gradBarrTDelta_,
    904                       reference_barr_);
    905   }
    906 
    907   Number FilterLineSearch::CalculateAlphaMin()
    908   {
    909     Number gBD = IpCq().curr_gradBarrTDelta();
    910     Number curr_theta = IpCq().curr_constraint_violation();
    911     Number alpha_min = gamma_theta_;
    912 
    913     if (gBD < 0) {
    914       alpha_min = Min( gamma_theta_,
    915                        gamma_phi_*curr_theta/(-gBD));
    916       if (curr_theta <= theta_min_) {
    917         alpha_min = Min( alpha_min,
    918                          delta_*pow(curr_theta,s_theta_)/pow(-gBD,s_phi_)
    919                        );
    920       }
    921     }
    922 
    923     return alpha_min_frac_*alpha_min;
    924   }
    925 
    926   bool FilterLineSearch::IsAcceptableToCurrentIterate(Number trial_barr,
    927       Number trial_theta,
    928       bool called_from_restoration /*=false*/) const
    929   {
    930     DBG_START_METH("FilterLineSearch::IsAcceptableToCurrentIterate",
    931                    dbg_verbosity);
    932 
    933     // Check if the barrier objective function is increasing to
    934     // rapidly (according to option obj_max_inc)
    935     if (!called_from_restoration && trial_barr > reference_barr_) {
    936       Number basval = 1.;
    937       if (fabs(reference_barr_)>10.) {
    938         basval = log10(fabs(reference_barr_));
    939       }
    940       if (log10(trial_barr-reference_barr_)>obj_max_inc_+basval) {
    941         Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Rejecting trial point because barrier objective function increasing too rapidly (from %27.15e to %27.15e)\n",reference_barr_,trial_barr);
    942         return false;
    943       }
    944     }
    945 
    946     DBG_PRINT((1,"trial_barr  = %e reference_barr  = %e\n", trial_barr, reference_barr_));
    947     DBG_PRINT((1,"trial_theta = %e reference_theta = %e\n", trial_theta, reference_theta_));
    948     return (Compare_le(trial_theta, (1.-gamma_theta_)*reference_theta_, reference_theta_)
    949             || Compare_le(trial_barr-reference_barr_, -gamma_phi_*reference_theta_, reference_barr_));
    950   }
    951 
    952   bool FilterLineSearch::IsAcceptableToCurrentFilter(Number trial_barr, Number trial_theta) const
    953   {
    954     return filter_.Acceptable(trial_barr, trial_theta);
    955   }
    956 
    957   bool FilterLineSearch::Compare_le(Number lhs, Number rhs, Number BasVal)
    958   {
    959     DBG_START_FUN("FilterLineSearch::Compare_le",
    960                   dbg_verbosity);
    961     DBG_PRINT((1,"lhs = %27.16e rhs = %27.16e  BasVal = %27.16e\n",lhs,rhs,BasVal));
    962     // ToDo: Comparison based on machine precision
    963     return (lhs - rhs <= 1e-15*fabs(BasVal));
    964   }
    965 
    966   bool FilterLineSearch::TrySoftRestoStep(SmartPtr<IteratesVector>& actual_delta,
    967                                           bool &satisfies_original_filter)
    968   {
    969     DBG_START_FUN("FilterLineSearch::TrySoftRestoStep", dbg_verbosity);
    970 
    971     satisfies_original_filter = false;
    972 
    973     // ToDo: Need to decide if we want to try a corrector step first
    974 
    975     // Compute the maximal step sizes (we use identical step sizes for
    976     // primal and dual variables
    977     Number alpha_primal_max =
    978       IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
    979                                       *actual_delta->x(),
    980                                       *actual_delta->s());
    981     Number alpha_dual_max =
    982       IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
    983                                     *actual_delta->z_L(),
    984                                     *actual_delta->z_U(),
    985                                     *actual_delta->v_L(),
    986                                     *actual_delta->v_U());
    987     Number alpha_max =  Min(alpha_primal_max, alpha_dual_max);
    988 
    989     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    990                    "Trying soft restoration phase step with step length %13.6e\n",
    991                    alpha_max);
    992 
    993     // Set the trial point
    994     IpData().SetTrialPrimalVariablesFromStep(alpha_max, *actual_delta->x(), *actual_delta->s());
    995     PerformDualStep(alpha_max, alpha_max,
    996                     actual_delta);
    997 
    998     // Check if that point is acceptable with respect to the current
    999     // original filter
    1000 
    1001     Number trial_barr;
    1002     Number trial_theta;
    1003     try {
    1004       trial_barr = IpCq().trial_barrier_obj();
    1005       trial_theta = IpCq().trial_constraint_violation();
    1006     }
    1007     catch(IpoptNLP::Eval_Error& e) {
    1008       e.ReportException(Jnlst(), J_DETAILED);
    1009       Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
    1010                      "Warning: Evaluation error during soft restoration phase step.\n");
    1011       IpData().Append_info_string("e");
    1012       return false;
    1013     }
    1014     if (theta_max_<=0 || trial_theta<=theta_max_) {
    1015       if (IsAcceptableToCurrentIterate(trial_barr, trial_theta)) {
    1016         Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1017                        "  Trial step acceptable with respect to original filter.\n");
    1018         satisfies_original_filter = true;
    1019         return true;
    1020       }
    1021     }
    1022 
    1023     // Evaluate the optimality error at the new point
    1024     Number mu = .0;
    1025     if (!IpData().FreeMuMode()) {
    1026       mu = IpData().curr_mu();
    1027     }
    1028     Number trial_pderror;
    1029     Number curr_pderror;
    1030     try {
    1031       trial_pderror = IpCq().trial_primal_dual_system_error(mu);
    1032       curr_pderror = IpCq().curr_primal_dual_system_error(mu);
    1033     }
    1034     catch(IpoptNLP::Eval_Error& e) {
    1035       e.ReportException(Jnlst(), J_DETAILED);
    1036       Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
    1037                      "Warning: Evaluation error during soft restoration phase step.\n");
    1038       IpData().Append_info_string("e");
    1039       return false;
    1040     }
    1041 
    1042     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1043                    "  Primal-dual error at current point:  %23.16e\n", curr_pderror);
    1044     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1045                    "  Primal-dual error at trial point  :  %23.16e\n", trial_pderror);
    1046     // Check if there is sufficient reduction in the optimality error
    1047     if (trial_pderror <= soft_resto_pderror_reduction_factor_*curr_pderror) {
    1048       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1049                      "  Trial step accepted.\n");
    1050       return true;
    1051     }
    1052 
    1053     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1054                    "  Trial step rejected.\n");
    1055     return false;
    1056   }
    1057 
    1058   void FilterLineSearch::StartWatchDog()
    1059   {
    1060     DBG_START_FUN("FilterLineSearch::StartWatchDog", dbg_verbosity);
     667  void BacktrackingLineSearch::StartWatchDog()
     668  {
     669    DBG_START_FUN("BacktrackingLineSearch::StartWatchDog", dbg_verbosity);
    1061670
    1062671    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
     
    1069678    watchdog_alpha_primal_test_ =
    1070679      IpCq().curr_primal_frac_to_the_bound(IpData().curr_tau());
    1071     watchdog_theta_ = IpCq().curr_constraint_violation();
    1072     watchdog_barr_ = IpCq().curr_barrier_obj();
    1073     watchdog_gradBarrTDelta_ = IpCq().curr_gradBarrTDelta();
    1074   }
    1075 
    1076   void FilterLineSearch::StopWatchDog(SmartPtr<IteratesVector>& actual_delta)
    1077   {
    1078     DBG_START_FUN("FilterLineSearch::StopWatchDog", dbg_verbosity);
     680
     681    acceptor_->StartWatchDog();
     682  }
     683
     684  void BacktrackingLineSearch::StopWatchDog(SmartPtr<IteratesVector>& actual_delta)
     685  {
     686    DBG_START_FUN("BacktrackingLineSearch::StopWatchDog", dbg_verbosity);
    1079687
    1080688    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
     
    1098706    watchdog_shortened_iter_ = 0;
    1099707
    1100     reference_theta_ = watchdog_theta_;
    1101     reference_barr_ = watchdog_barr_;
    1102     reference_gradBarrTDelta_ = watchdog_gradBarrTDelta_;
    1103   }
    1104 
    1105   void FilterLineSearch::Reset()
    1106   {
    1107     DBG_START_FUN("FilterLineSearch::Reset", dbg_verbosity);
     708    acceptor_->StopWatchDog();
     709  }
     710
     711  void BacktrackingLineSearch::Reset()
     712  {
     713    DBG_START_FUN("BacktrackingLineSearch::Reset", dbg_verbosity);
    1108714    in_soft_resto_phase_ = false;
    1109715
     
    1114720    watchdog_shortened_iter_ = 0;
    1115721
    1116     filter_.Clear();
    1117   }
    1118 
    1119   void FilterLineSearch::PerformDualStep(Number alpha_primal,
    1120                                          Number alpha_dual,
    1121                                          SmartPtr<IteratesVector>& delta)
    1122   {
    1123     DBG_START_FUN("FilterLineSearch::PerformDualStep", dbg_verbosity);
     722    acceptor_->Reset();
     723  }
     724
     725  void BacktrackingLineSearch::PerformDualStep(Number alpha_primal,
     726      Number alpha_dual,
     727      SmartPtr<IteratesVector>& delta)
     728  {
     729    DBG_START_FUN("BacktrackingLineSearch::PerformDualStep", dbg_verbosity);
    1124730
    1125731    // set the bound multipliers from the step
     
    1191797  }
    1192798
    1193   bool
    1194   FilterLineSearch::TrySecondOrderCorrection(
    1195     Number alpha_primal_test,
    1196     Number& alpha_primal,
    1197     SmartPtr<IteratesVector>& actual_delta)
    1198   {
    1199     DBG_START_METH("FilterLineSearch::TrySecondOrderCorrection",
    1200                    dbg_verbosity);
    1201 
    1202     bool accept = false;
    1203     Index count_soc = 0;
    1204 
    1205     Number theta_soc_old = 0.;
    1206     Number theta_trial = IpCq().trial_constraint_violation();
    1207     Number alpha_primal_soc = alpha_primal;
    1208 
    1209     SmartPtr<Vector> c_soc = IpCq().curr_c()->MakeNew();
    1210     SmartPtr<Vector> dms_soc = IpCq().curr_d_minus_s()->MakeNew();
    1211     c_soc->Copy(*IpCq().curr_c());
    1212     dms_soc->Copy(*IpCq().curr_d_minus_s());
    1213     while (count_soc<max_soc_ && !accept &&
    1214            (count_soc==0 || theta_trial<=kappa_soc_*theta_soc_old) ) {
    1215       theta_soc_old = theta_trial;
    1216 
    1217       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1218                      "Trying second order correction number %d\n",
    1219                      count_soc+1);
    1220 
    1221       // Compute SOC constraint violation
    1222       c_soc->AddOneVector(1.0, *IpCq().trial_c(), alpha_primal_soc);
    1223       dms_soc->AddOneVector(1.0, *IpCq().trial_d_minus_s(), alpha_primal_soc);
    1224 
    1225       // Compute the SOC search direction
    1226       SmartPtr<IteratesVector> delta_soc = actual_delta->MakeNewIteratesVector(true);
    1227       SmartPtr<IteratesVector> rhs = actual_delta->MakeNewContainer();
    1228       rhs->Set_x(*IpCq().curr_grad_lag_with_damping_x());
    1229       rhs->Set_s(*IpCq().curr_grad_lag_with_damping_s());
    1230       rhs->Set_y_c(*c_soc);
    1231       rhs->Set_y_d(*dms_soc);
    1232       rhs->Set_z_L(*IpCq().curr_relaxed_compl_x_L());
    1233       rhs->Set_z_U(*IpCq().curr_relaxed_compl_x_U());
    1234       rhs->Set_v_L(*IpCq().curr_relaxed_compl_s_L());
    1235       rhs->Set_v_U(*IpCq().curr_relaxed_compl_s_U());
    1236       pd_solver_->Solve(-1.0, 0.0, *rhs, *delta_soc, true);
    1237 
    1238       // Compute step size
    1239       alpha_primal_soc =
    1240         IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
    1241                                         *delta_soc->x(),
    1242                                         *delta_soc->s());
    1243 
    1244       // Check if trial point is acceptable
    1245       try {
    1246         // Compute the primal trial point
    1247         IpData().SetTrialPrimalVariablesFromStep(alpha_primal_soc, *delta_soc->x(), *delta_soc->s());
    1248 
    1249         // in acceptance tests, use original step size!
    1250         accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test);
    1251       }
    1252       catch(IpoptNLP::Eval_Error& e) {
    1253         e.ReportException(Jnlst(), J_DETAILED);
    1254         Jnlst().Printf(J_WARNING, J_MAIN, "Warning: SOC step rejected due to evaluation error\n");
    1255         IpData().Append_info_string("e");
    1256         accept = false;
    1257         // There is no point in continuing SOC procedure
    1258         break;
    1259       }
    1260 
    1261       if (accept) {
    1262         Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Second order correction step accepted with %d corrections.\n", count_soc+1);
    1263         // Accept all SOC quantities
    1264         alpha_primal = alpha_primal_soc;
    1265         actual_delta = delta_soc;
    1266       }
    1267       else {
    1268         count_soc++;
    1269         theta_trial = IpCq().trial_constraint_violation();
    1270       }
    1271     }
    1272     return accept;
    1273   }
    1274 
    1275   bool
    1276   FilterLineSearch::TryCorrector(
    1277     Number alpha_primal_test,
    1278     Number& alpha_primal,
    1279     SmartPtr<IteratesVector>& actual_delta)
    1280   {
    1281     DBG_START_METH("FilterLineSearch::TryCorrector",
    1282                    dbg_verbosity);
    1283 
    1284     Index n_bounds = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim()
    1285                      + IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim();
    1286     if (n_bounds==0) {
    1287       // Nothing to be done
    1288       return false;
    1289     }
    1290 
    1291     IpData().TimingStats().TryCorrector().Start();
    1292 
    1293     bool accept = false;
    1294 
    1295     // Compute the corrector step based on corrector_type parameter
    1296     // create a new iterates vector and allocate space for all the entries
    1297     SmartPtr<IteratesVector> delta_corr = actual_delta->MakeNewIteratesVector(true);
    1298 
    1299     switch (corrector_type_) {
    1300       case AFFINE_CORRECTOR : {
    1301         // 1: Standard MPC corrector
    1302 
    1303         if (!IpData().HaveAffineDeltas()) {
    1304           Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1305                          "Solving the Primal Dual System for the affine step\n");
    1306           // First get the right hand side
    1307           SmartPtr<IteratesVector> rhs_aff = delta_corr->MakeNewContainer();
    1308 
    1309           rhs_aff->Set_x(*IpCq().curr_grad_lag_x());
    1310           rhs_aff->Set_s(*IpCq().curr_grad_lag_s());
    1311           rhs_aff->Set_y_c(*IpCq().curr_c());
    1312           rhs_aff->Set_y_d(*IpCq().curr_d_minus_s());
    1313           rhs_aff->Set_z_L(*IpCq().curr_compl_x_L());
    1314           rhs_aff->Set_z_U(*IpCq().curr_compl_x_U());
    1315           rhs_aff->Set_v_L(*IpCq().curr_compl_s_L());
    1316           rhs_aff->Set_v_U(*IpCq().curr_compl_s_U());
    1317 
    1318           // create a new iterates vector (with allocated space)
    1319           // for the affine scaling step
    1320           SmartPtr<IteratesVector> step_aff = delta_corr->MakeNewIteratesVector(true);
    1321 
    1322           // Now solve the primal-dual system to get the step
    1323           pd_solver_->Solve(-1.0, 0.0, *rhs_aff, *step_aff, false);
    1324 
    1325           DBG_PRINT_VECTOR(2, "step_aff", *step_aff);
    1326 
    1327           IpData().set_delta_aff(step_aff);
    1328           IpData().SetHaveAffineDeltas(true);
    1329         }
    1330 
    1331         DBG_ASSERT(IpData().HaveAffineDeltas());
    1332 
    1333         const SmartPtr<const IteratesVector> delta_aff = IpData().delta_aff();
    1334 
    1335         delta_corr->Copy(*actual_delta);
    1336 
    1337         // create a rhs vector and allocate entries
    1338         SmartPtr<IteratesVector> rhs = actual_delta->MakeNewIteratesVector(true);
    1339 
    1340         rhs->x_NonConst()->Set(0.);
    1341         rhs->s_NonConst()->Set(0.);
    1342         rhs->y_c_NonConst()->Set(0.);
    1343         rhs->y_d_NonConst()->Set(0.);
    1344         IpNLP().Px_L()->TransMultVector(-1., *delta_aff->x(), 0., *rhs->z_L_NonConst());
    1345         rhs->z_L_NonConst()->ElementWiseMultiply(*delta_aff->z_L());
    1346         IpNLP().Px_U()->TransMultVector(1., *delta_aff->x(), 0., *rhs->z_U_NonConst());
    1347         rhs->z_U_NonConst()->ElementWiseMultiply(*delta_aff->z_U());
    1348         IpNLP().Pd_L()->TransMultVector(-1., *delta_aff->s(), 0., *rhs->v_L_NonConst());
    1349         rhs->v_L_NonConst()->ElementWiseMultiply(*delta_aff->v_L());
    1350         IpNLP().Pd_U()->TransMultVector(1., *delta_aff->s(), 0., *rhs->v_U_NonConst());
    1351         rhs->v_U_NonConst()->ElementWiseMultiply(*delta_aff->v_U());
    1352 
    1353         pd_solver_->Solve(1.0, 1.0, *rhs, *delta_corr, true);
    1354 
    1355         DBG_PRINT_VECTOR(2, "delta_corr", *delta_corr);
    1356       }
    1357       break;
    1358       case PRIMAL_DUAL_CORRECTOR : {
    1359         // 2: Second order correction for primal-dual step to
    1360         // primal-dual mu
    1361 
    1362         delta_corr->Copy(*actual_delta);
    1363 
    1364         // allocate space for the rhs
    1365         SmartPtr<IteratesVector> rhs = actual_delta->MakeNewIteratesVector(true);
    1366 
    1367         rhs->x_NonConst()->Set(0.);
    1368         rhs->s_NonConst()->Set(0.);
    1369         rhs->y_c_NonConst()->Set(0.);
    1370         rhs->y_d_NonConst()->Set(0.);
    1371 
    1372         Number mu = IpData().curr_mu();
    1373         SmartPtr<Vector> tmp;
    1374 
    1375         rhs->z_L_NonConst()->Copy(*IpCq().curr_slack_x_L());
    1376         IpNLP().Px_L()->TransMultVector(-1., *actual_delta->x(),
    1377                                         -1., *rhs->z_L_NonConst());
    1378         tmp = actual_delta->z_L()->MakeNew();
    1379         tmp->AddTwoVectors(1., *IpData().curr()->z_L(), 1., *actual_delta->z_L(), 0.);
    1380         rhs->z_L_NonConst()->ElementWiseMultiply(*tmp);
    1381         rhs->z_L_NonConst()->AddScalar(mu);
    1382 
    1383         rhs->z_U_NonConst()->Copy(*IpCq().curr_slack_x_U());
    1384         IpNLP().Px_U()->TransMultVector(1., *actual_delta->x(),
    1385                                         -1., *rhs->z_U_NonConst());
    1386         tmp = actual_delta->z_U()->MakeNew();
    1387         tmp->AddTwoVectors(1., *IpData().curr()->z_U(), 1., *actual_delta->z_U(), 0.);
    1388         rhs->z_U_NonConst()->ElementWiseMultiply(*tmp);
    1389         rhs->z_U_NonConst()->AddScalar(mu);
    1390 
    1391         rhs->v_L_NonConst()->Copy(*IpCq().curr_slack_s_L());
    1392         IpNLP().Pd_L()->TransMultVector(-1., *actual_delta->s(),
    1393                                         -1., *rhs->v_L_NonConst());
    1394         tmp = actual_delta->v_L()->MakeNew();
    1395         tmp->AddTwoVectors(1., *IpData().curr()->v_L(), 1., *actual_delta->v_L(), 0.);
    1396         rhs->v_L_NonConst()->ElementWiseMultiply(*tmp);
    1397         rhs->v_L_NonConst()->AddScalar(mu);
    1398 
    1399         rhs->v_U_NonConst()->Copy(*IpCq().curr_slack_s_U());
    1400         IpNLP().Pd_U()->TransMultVector(1., *actual_delta->s(),
    1401                                         -1., *rhs->v_U_NonConst());
    1402         tmp = actual_delta->v_U()->MakeNew();
    1403         tmp->AddTwoVectors(1., *IpData().curr()->v_U(), 1., *actual_delta->v_U(), 0.);
    1404         rhs->v_U_NonConst()->ElementWiseMultiply(*tmp);
    1405         rhs->v_U_NonConst()->AddScalar(mu);
    1406 
    1407         DBG_PRINT_VECTOR(2, "rhs", *rhs);
    1408 
    1409         pd_solver_->Solve(1.0, 1.0, *rhs, *delta_corr, true);
    1410 
    1411         DBG_PRINT_VECTOR(2, "delta_corr", *delta_corr);
    1412       }
    1413       break;
    1414       default:
    1415       DBG_ASSERT("Unknown corrector_type value.");
    1416     }
    1417 
    1418     // Compute step size
    1419     Number alpha_primal_corr =
    1420       IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
    1421                                       *delta_corr->x(),
    1422                                       *delta_corr->s());
    1423     // Set the primal trial point
    1424     IpData().SetTrialPrimalVariablesFromStep(alpha_primal_corr, *delta_corr->x(), *delta_corr->s());
    1425 
    1426     // Check if we want to not even try the filter criterion
    1427     Number alpha_dual_max =
    1428       IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
    1429                                     *delta_corr->z_L(), *delta_corr->z_U(),
    1430                                     *delta_corr->v_L(), *delta_corr->v_U());
    1431 
    1432     IpData().SetTrialBoundMultipliersFromStep(alpha_dual_max, *delta_corr->z_L(), *delta_corr->z_U(), *delta_corr->v_L(), *delta_corr->v_U());
    1433 
    1434     Number trial_avrg_compl = IpCq().trial_avrg_compl();
    1435     Number curr_avrg_compl = IpCq().curr_avrg_compl();
    1436     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1437                    "avrg_compl(curr) = %e, avrg_compl(trial) = %e\n",
    1438                    curr_avrg_compl, trial_avrg_compl);
    1439     if (corrector_type_==AFFINE_CORRECTOR &&
    1440         trial_avrg_compl>=corrector_compl_avrg_red_fact_*curr_avrg_compl) {
    1441       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1442                      "Rejecting corrector step, because trial complementarity is too large.\n" );
    1443       IpData().TimingStats().TryCorrector().End();
    1444       return false;
    1445     }
    1446 
    1447     // Check if trial point is acceptable
    1448     try {
    1449       // in acceptance tests, use original step size!
    1450       accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test);
    1451     }
    1452     catch(IpoptNLP::Eval_Error& e) {
    1453       e.ReportException(Jnlst(), J_DETAILED);
    1454       Jnlst().Printf(J_WARNING, J_MAIN,
    1455                      "Warning: Corrector step rejected due to evaluation error\n");
    1456       IpData().Append_info_string("e");
    1457       accept = false;
    1458     }
    1459 
    1460     if (accept) {
    1461       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1462                      "Corrector step accepted with alpha_primal = %e\n",
    1463                      alpha_primal_corr);
    1464       // Accept all SOC quantities
    1465       alpha_primal = alpha_primal_corr;
    1466       actual_delta = delta_corr;
    1467 
    1468       if (Jnlst().ProduceOutput(J_MOREVECTOR, J_MAIN)) {
    1469         Jnlst().Printf(J_MOREVECTOR, J_MAIN,
    1470                        "*** Accepted corrector for Iteration: %d\n",
    1471                        IpData().iter_count());
    1472         delta_corr->Print(Jnlst(), J_MOREVECTOR, J_MAIN, "delta_corr");
    1473       }
    1474     }
    1475 
    1476     IpData().TimingStats().TryCorrector().End();
    1477     return accept;
    1478   }
    1479 
    1480799  void
    1481   FilterLineSearch::PerformMagicStep()
    1482   {
    1483     DBG_START_METH("FilterLineSearch::PerformMagicStep",
     800  BacktrackingLineSearch::PerformMagicStep()
     801  {
     802    DBG_START_METH("BacktrackingLineSearch::PerformMagicStep",
    1484803                   2);//dbg_verbosity);
    1485804
     
    1578897  }
    1579898
     899  bool BacktrackingLineSearch::TrySoftRestoStep(SmartPtr<IteratesVector>& actual_delta,
     900      bool &satisfies_original_criterion)
     901  {
     902    DBG_START_FUN("FilterLSAcceptor::TrySoftRestoStep", dbg_verbosity);
     903
     904    if (soft_resto_pderror_reduction_factor_==0.) {
     905      return false;
     906    }
     907
     908    satisfies_original_criterion = false;
     909
     910    // ToDo: Need to decide if we want to try a corrector step first
     911
     912    // Compute the maximal step sizes (we use identical step sizes for
     913    // primal and dual variables
     914    Number alpha_primal_max =
     915      IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
     916                                      *actual_delta->x(),
     917                                      *actual_delta->s());
     918    Number alpha_dual_max =
     919      IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
     920                                    *actual_delta->z_L(),
     921                                    *actual_delta->z_U(),
     922                                    *actual_delta->v_L(),
     923                                    *actual_delta->v_U());
     924    Number alpha_max =  Min(alpha_primal_max, alpha_dual_max);
     925
     926    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
     927                   "Trying soft restoration phase step with step length %13.6e\n",
     928                   alpha_max);
     929
     930    // Set the trial point
     931    IpData().SetTrialPrimalVariablesFromStep(alpha_max, *actual_delta->x(), *actual_delta->s());
     932    PerformDualStep(alpha_max, alpha_max,
     933                    actual_delta);
     934
     935    // Check if that point is acceptable with respect to the current
     936    // original filter
     937
     938    Number trial_barr;
     939    Number trial_theta;
     940    try {
     941      trial_barr = IpCq().trial_barrier_obj();
     942      trial_theta = IpCq().trial_constraint_violation();
     943    }
     944    catch(IpoptNLP::Eval_Error& e) {
     945      e.ReportException(Jnlst(), J_DETAILED);
     946      Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
     947                     "Warning: Evaluation error during soft restoration phase step.\n");
     948      IpData().Append_info_string("e");
     949      return false;
     950    }
     951    if (acceptor_->CheckAcceptabilityOfTrialPoint(0.)) {
     952      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
     953                     "  Trial step acceptable with respect to original backtracking globalization.\n");
     954      satisfies_original_criterion = true;
     955      return true;
     956    }
     957
     958    // Evaluate the optimality error at the new point
     959    Number mu = .0;
     960    if (!IpData().FreeMuMode()) {
     961      mu = IpData().curr_mu();
     962    }
     963    Number trial_pderror;
     964    Number curr_pderror;
     965    try {
     966      trial_pderror = IpCq().trial_primal_dual_system_error(mu);
     967      curr_pderror = IpCq().curr_primal_dual_system_error(mu);
     968    }
     969    catch(IpoptNLP::Eval_Error& e) {
     970      e.ReportException(Jnlst(), J_DETAILED);
     971      Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
     972                     "Warning: Evaluation error during soft restoration phase step.\n");
     973      IpData().Append_info_string("e");
     974      return false;
     975    }
     976
     977    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
     978                   "  Primal-dual error at current point:  %23.16e\n", curr_pderror);
     979    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
     980                   "  Primal-dual error at trial point  :  %23.16e\n", trial_pderror);
     981    // Check if there is sufficient reduction in the optimality error
     982    if (trial_pderror <= soft_resto_pderror_reduction_factor_*curr_pderror) {
     983      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
     984                     "  Trial step accepted.\n");
     985      return true;
     986    }
     987
     988    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
     989                   "  Trial step rejected.\n");
     990    return false;
     991  }
     992
    1580993  bool
    1581   FilterLineSearch::DetectTinyStep()
    1582   {
    1583     DBG_START_METH("FilterLineSearch::DetectTinyStep",
     994  BacktrackingLineSearch::DetectTinyStep()
     995  {
     996    DBG_START_METH("BacktrackingLineSearch::DetectTinyStep",
    1584997                   dbg_verbosity);
    1585998
     
    16291042  }
    16301043
    1631   bool FilterLineSearch::CurrentIsAcceptable()
     1044  bool BacktrackingLineSearch::CurrentIsAcceptable()
    16321045  {
    16331046    return (IsValid(conv_check_) &&
     
    16351048  }
    16361049
    1637   void FilterLineSearch::StoreAcceptablePoint()
    1638   {
    1639     DBG_START_METH("FilterLineSearch::StoreAcceptablePoint",
     1050  void BacktrackingLineSearch::StoreAcceptablePoint()
     1051  {
     1052    DBG_START_METH("BacktrackingLineSearch::StoreAcceptablePoint",
    16401053                   dbg_verbosity);
    16411054
     
    16431056  }
    16441057
    1645   bool FilterLineSearch::RestoreAcceptablePoint()
    1646   {
    1647     DBG_START_METH("FilterLineSearch::RestoreAcceptablePoint",
     1058  bool BacktrackingLineSearch::RestoreAcceptablePoint()
     1059  {
     1060    DBG_START_METH("BacktrackingLineSearch::RestoreAcceptablePoint",
    16481061                   dbg_verbosity);
    16491062
  • branches/dev/Algorithm/IpBacktrackingLineSearch.hpp

    r541 r542  
    66//
    77// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
    8 
    9 #ifndef __IPFILTERLINESEARCH_HPP__
    10 #define __IPFILTERLINESEARCH_HPP__
    11 
    12 #include "IpFilter.hpp"
     8//           Andreas Waechter                 IBM    2005-10-13
     9//               derived file from IpFilterLineSearch.hpp
     10
     11#ifndef __IPBACKTRACKINGLINESEARCH_HPP__
     12#define __IPBACKTRACKINGLINESEARCH_HPP__
     13
    1314#include "IpLineSearch.hpp"
     15#include "IpBacktrackingLSAcceptor.hpp"
    1416#include "IpRestoPhase.hpp"
    15 #include "IpPDSystemSolver.hpp"
    1617#include "IpConvCheck.hpp"
    1718
     
    1920{
    2021
    21   /** Filter line search.  This class implements the filter line
    22    *  search procedure.
     22  /** General implementation of a backtracking line search.  This
     23   *  class can be used to perform the filter line search procedure or
     24   *  other procedures.  The BacktrackingLSAcceptor is used to
     25   *  determine whether trial points are acceptable (e.g., based on a
     26   *  filter or other methods).
     27   *
     28   *  This backtracking line search knows of a restoration phase
     29   *  (which is called when the trial step size becomes too small or
     30   *  no search direction could be computed).  It also has the notion
     31   *  of a "soft restoration phase," which uses the regular steps but
     32   *  decides on the acceptability based on other measures than the
     33   *  regular ones (e.g., reduction of the PD error instead of
     34   *  acceptability to a filter mechanism).
    2335   */
    24   class FilterLineSearch : public LineSearch
     36  class BacktrackingLineSearch : public LineSearch
    2537  {
    2638  public:
    2739    /**@name Constructors/Destructors */
    2840    //@{
    29     /** Constructor.  The PDSystemSolver object only needs to be
     41    /** Constructor.  The acceptor implements the acceptance test for
     42     *  the line search. The PDSystemSolver object only needs to be
    3043     *  provided (i.e. not NULL) if second order correction is to be
    3144     *  used.  The ConvergenceCheck object is used to determine
     
    3346     *  restoration phase is not started if the acceptability level
    3447     *  has been reached).  If conv_check is NULL, we assume that the
    35      *  current iterate is not acceptable. */
    36     FilterLineSearch(const SmartPtr<RestorationPhase>& resto_phase,
    37                      const SmartPtr<PDSystemSolver>& pd_solver,
    38                      const SmartPtr<ConvergenceCheck>& conv_check
    39                     );
     48     *  current iterate is not acceptable (in the sense of the
     49     *  acceptable_tol option). */
     50    BacktrackingLineSearch(const SmartPtr<BacktrackingLSAcceptor>& acceptor,
     51                           const SmartPtr<RestorationPhase>& resto_phase,
     52                           const SmartPtr<ConvergenceCheck>& conv_check
     53                          );
    4054
    4155    /** Default destructor */
    42     virtual ~FilterLineSearch();
     56    virtual ~BacktrackingLineSearch();
    4357    //@}
    4458
     
    8599    }
    86100
    87     /**@name Trial Point Accepting Methods. Used internally to check certain
    88      * acceptability criteria and used externally (by the restoration phase
    89      * convergence check object, for instance)
    90      */
    91     //@{
    92     /** Checks if a trial point is acceptable to the current iterate */
    93     bool IsAcceptableToCurrentIterate(Number trial_barr, Number trial_theta,
    94                                       bool called_from_restoration=false) const;
    95 
    96     /** Checks if a trial point is acceptable to the current filter */
    97     bool IsAcceptableToCurrentFilter(Number trial_barr, Number trial_theta) const;
    98     //@}
    99 
    100101    /** Methods for OptionsList */
    101102    //@{
     
    113114    //@{
    114115    /** Copy Constructor */
    115     FilterLineSearch(const FilterLineSearch&);
     116    BacktrackingLineSearch(const BacktrackingLineSearch&);
    116117
    117118    /** Overloaded Equals Operator */
    118     void operator=(const FilterLineSearch&);
    119     //@}
    120 
    121     /** @name Filter information */
    122     //@{
    123     /** Upper bound on infeasibility */
    124     Number theta_max_;
    125     Number theta_max_fact_;
    126 
    127     /** Infeasibility switching bound */
    128     Number theta_min_;
    129     Number theta_min_fact_;
    130     //@}
    131 
    132     /** Method for checking if the current step size satisfies the
    133      *  f-type switching condition.  Here, we use the search direction
    134      *  stored in ip_data
    135      */
    136     bool IsFtype(Number alpha_primal_test);
    137 
    138     /** Method for checking the Armijo condition, given a trial step
    139      *  size.  The test uses the search direction stored in ip_data,
    140      *  and the values of the functions at the trial point in ip_data.
    141      */
    142     bool ArmijoHolds(Number alpha_primal_test);
    143 
    144     /** Method to calculate alpha_min (minimum alpha before going to
    145      *  restoration
    146      */
    147     Number CalculateAlphaMin();
    148 
    149     /** Augment the filter used on the current values of the barrier
    150      *  objective function and the contraint violation */
    151     void AugmentFilter();
     119    void operator=(const BacktrackingLineSearch&);
     120    //@}
    152121
    153122    /** Method performing the backtracking line search.  The return
     
    194163    /** Try a step for the soft restoration phase and check if it is
    195164     *  acceptable.  The step size is identical for all variables.  A
    196      *  point is accepted if it is acceptable for the original filter
    197      *  (in which case satisfies_original_filter = true on return), or
    198      *  if the primal-dual system error was decrease by at least the
    199      *  factor soft_resto_pderror_reduction_factor_.  The return value is
    200      *  true, if the trial point was acceptable. */
     165     *  point is accepted if it is acceptable for the original
     166     *  acceptability criterion (in which case
     167     *  satisfies_original_criterion = true on return), or if the
     168     *  primal-dual system error was decrease by at least the factor
     169     *  soft_resto_pderror_reduction_factor_.  The return value is
     170     *  true, if the trial point was acceptable for the soft
     171     *  restoration phase. */
    201172    bool TrySoftRestoStep(SmartPtr<IteratesVector>& actual_delta,
    202                           bool &satisfies_original_filter);
     173                          bool &satisfies_original_criterion);
    203174
    204175    /** Try a second order correction for the constraints.  If the
     
    245216    bool RestoreAcceptablePoint();
    246217
    247     /** Method for determining if the current iterate is acceptable.
    248      *  This is a wrapper for same method from ConvergenceCheck, but
    249      *  returns false, if no ConvergenceCheck object is provided. */
     218    /** Method for determining if the current iterate is acceptable
     219     *  (in the sense of the acceptable_tol options).  This is a
     220     *  wrapper for same method from ConvergenceCheck, but returns
     221     *  false, if no ConvergenceCheck object is provided. */
    250222    bool CurrentIsAcceptable();
    251223
    252224    /** @name Parameters for the filter algorithm.  Names as in the paper */
    253225    //@{
    254     /** \f$ \eta_{\varphi} \f$ */
    255     Number eta_phi_;
    256     /** \f$ \delta \f$ */
    257     Number delta_;
    258     /** \f$ s_{\varphi} \f$ */
    259     Number s_phi_;
    260     /** \f$ s_{\Theta} \f$ */
    261     Number s_theta_;
    262     /** \f$ \gamma_{\varphi} \f$ */
    263     Number gamma_phi_;
    264     /** \f$ \gamma_{\Theta} \f$ */
    265     Number gamma_theta_;
    266     /** \f$ \gamma_{\alpha} \f$ */
    267     Number alpha_min_frac_;
    268226    /** factor by which search direction is to be shortened if trial
    269227     *  point is rejected. */
    270228    Number alpha_red_factor_;
    271     /** Maximal number of second order correction steps */
    272     Index max_soc_;
    273     /** Required reduction in constraint violation before trying
    274      *  multiple second order correction steps \f$ \kappa_{soc}\f$.
    275      */
    276     Number kappa_soc_;
    277     /** Maximal increase in objective function in orders of magnitute
    278      *  (log10).  If the log10(barrier objective function) is
    279      *  increased more than this compared to the current point, the
    280      *  trial point is rejected. */
    281     Number obj_max_inc_;
    282229
    283230    /** enumeration for the different alpha_for_y_ settings */
     
    297244    AlphaForYEnum alpha_for_y_;
    298245
     246    /** Reduction factor for the restoration phase that accepts steps
     247     *  reducing the optimality error ("soft restoration phase"). If
     248     *  0., then this restoration phase is not enabled. */
     249    Number soft_resto_pderror_reduction_factor_;
     250
    299251    /** Flag indicating whether magic steps should be used. */
    300252    bool magic_steps_;
    301     /** enumeration for the corrector type */
    302     enum CorrectorTypeEnum {
    303       NO_CORRECTOR=0,
    304       AFFINE_CORRECTOR,
    305       PRIMAL_DUAL_CORRECTOR
    306     };
    307     /** Type of corrector steps that should be tried. */
    308     CorrectorTypeEnum corrector_type_;
    309253    /** Flag indicating whether the line search should always accept
    310254     *  the full (fraction-to-the-boundary) step. */
    311255    bool accept_every_trial_step_;
    312     /** parameter in heurstic that determines whether corrector step
    313     should be tried. */
    314     Number corrector_compl_avrg_red_fact_;
    315     /** Flag indicating whether the corrector should be skipped in an
    316      *  iteration in which negative curvature is detected */
    317     bool skip_corr_if_neg_curv_;
    318     /** Flag indicating whether the corrector should be skipped during
    319      *  the monotone mu mode. */
    320     bool skip_corr_in_monotone_mode_;
    321256    /** Indicates whether problem can be expected to be infeasible.
    322257     *  This will trigger requesting a tighter reduction in
     
    329264     *  disabled for the rest of the optimization run. */
    330265    Number expect_infeasible_problem_ctol_;
    331     /** Reduction factor for the restoration phase that accepts steps
    332      *  reducing the optimality error ("soft restoration phase"). If
    333      *  0., then this restoration phase is not enabled. */
    334     Number soft_resto_pderror_reduction_factor_;
    335266
    336267    /** Tolerance for detecting tiny steps. */
     
    349280    /** @name Information related to watchdog procedure */
    350281    //@{
    351     /** Constraint violation at the point with respect to which
    352      *  progress is to be made */
    353     Number reference_theta_;
    354     /** Barrier objective function at the point with respect to which
    355      *  progress is to be made */
    356     Number reference_barr_;
    357     /** Barrier gradient transpose search direction at the point with
    358      *  respect to which progress is to be made */
    359     Number reference_gradBarrTDelta_;
    360282    /** Flag indicating if the watchdog is active */
    361283    bool in_watchdog_;
     
    366288    /** Step size for Armijo test in watch dog */
    367289    Number watchdog_alpha_primal_test_;
    368     /** Constraint violation at reference point */
    369     Number watchdog_theta_;
    370     /** Barrier objective function at reference point */
    371     Number watchdog_barr_;
    372     /** Barrier gradient transpose search direction at reference point */
    373     Number watchdog_gradBarrTDelta_;
    374290    /** Watchdog reference iterate */
    375291    SmartPtr<const IteratesVector> watchdog_iterate_;
     
    383299    SmartPtr<const IteratesVector> acceptable_iterate_;
    384300    //@}
    385 
    386     /** Filter with entries */
    387     Filter filter_;
    388301
    389302    /** Flag indicating whether the line search is to be performed
     
    413326    /** @name Strategy objective that are used */
    414327    //@{
     328    SmartPtr<BacktrackingLSAcceptor> acceptor_;
    415329    SmartPtr<RestorationPhase> resto_phase_;
    416     SmartPtr<PDSystemSolver> pd_solver_;
    417330    SmartPtr<ConvergenceCheck> conv_check_;
    418331    //@}
  • branches/dev/Algorithm/IpFilterLSAcceptor.cpp

    r541 r542  
    66//
    77// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
    8 
    9 #include "IpFilterLineSearch.hpp"
     8//           Andreas Waechter                 IBM    2005-10-13
     9//               derived file from IpFilterLineSearch.cpp
     10
     11#include "IpFilterLSAcceptor.hpp"
    1012#include "IpJournalist.hpp"
    1113#include "IpRestoPhase.hpp"
     
    2224#endif
    2325
    24 #ifdef HAVE_CCTYPE
    25 # include <cctype>
    26 #else
    27 # ifdef HAVE_CTYPE_H
    28 #  include <ctype.h>
    29 # else
    30 #  error "don't have header file for ctype"
    31 # endif
    32 #endif
    33 
    3426namespace Ipopt
    3527{
     
    3931#endif
    4032
    41   FilterLineSearch::FilterLineSearch(const SmartPtr<RestorationPhase>& resto_phase,
    42                                      const SmartPtr<PDSystemSolver>& pd_solver,
    43                                      const SmartPtr<ConvergenceCheck>& conv_check)
     33  FilterLSAcceptor::FilterLSAcceptor(const SmartPtr<PDSystemSolver>& pd_solver)
    4434      :
    45       LineSearch(),
    4635      filter_(2),
    47       resto_phase_(resto_phase),
    48       pd_solver_(pd_solver),
    49       conv_check_(conv_check)
    50   {
    51     DBG_START_FUN("FilterLineSearch::FilterLineSearch",
     36      pd_solver_(pd_solver)
     37  {
     38    DBG_START_FUN("FilterLSAcceptor::FilterLSAcceptor",
    5239                  dbg_verbosity);
    5340  }
    5441
    55   FilterLineSearch::~FilterLineSearch()
    56   {
    57     DBG_START_FUN("FilterLineSearch::~FilterLineSearch()",
     42  FilterLSAcceptor::~FilterLSAcceptor()
     43  {
     44    DBG_START_FUN("FilterLSAcceptor::~FilterLSAcceptor()",
    5845                  dbg_verbosity);
    5946  }
    6047
    61   void FilterLineSearch::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
     48  void FilterLSAcceptor::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
    6249  {
    6350    roptions->AddLowerBoundedNumberOption(
     
    11299      0.0, true, 1.0, true, 0.05,
    113100      "(This is gamma_alpha in Eqn. (20) in implementation paper.)");
    114     roptions->AddBoundedNumberOption(
    115       "alpha_red_factor",
    116       "Fractional reduction of the trial step size in the backtracking line search.",
    117       0.0, true, 1.0, true, 0.5,
    118       "At every step of the backtracking line search, the trial step size is "
    119       "reduced by this factor.");
    120101    roptions->AddLowerBoundedIntegerOption(
    121102      "max_soc",
     
    140121      "of magnitude.");
    141122
    142     std::string prev_category = roptions->RegisteringCategory();
    143     roptions->SetRegisteringCategory("Undocumented");
    144     roptions->AddStringOption2(
    145       "magic_steps",
    146       "Enables magic steps.",
    147       "no",
    148       "no", "don't take magic steps",
    149       "yes", "take magic steps",
    150       "DOESN'T REALLY WORK YET!");
    151     roptions->SetRegisteringCategory(prev_category);
    152123    roptions->AddStringOption3(
    153124      "corrector_type",
     
    181152      "This option is only used if \"mu_strategy\" is \"adaptive\".");
    182153
    183     roptions->AddStringOption2(
    184       "accept_every_trial_step",
    185       "Always accept the frist trial step.",
    186       "no",
    187       "no", "don't arbitrarily accept the full step",
    188       "yes", "always accept the full step",
    189       "Setting this option to \"yes\" essentially disables the line search "
    190       "and makes the algorithm take aggressive steps.");
    191 
    192     roptions->AddStringOption7(
    193       "alpha_for_y",
    194       "Method to determine the step size for constraint multipliers.",
    195       "primal",
    196       "primal", "use primal step size",
    197       "bound_mult", "use step size for the bound multipliers",
    198       "min", "use the min of primal and bound multipliers",
    199       "max", "use the max of primal and bound multipliers",
    200       "full", "take a full step of size one",
    201       "min_dual_infeas", "choose step size minimizing new dual infeasibility",
    202       "safe_min_dual_infeas", "like \"min_dual_infeas\", but safeguarded by \"min\" and \"max\"",
    203       "This option determines how the step size (alpha_y) will be calculated when updating the "
    204       "constraint multipliers.");
    205 
    206154    roptions->AddLowerBoundedNumberOption(
    207155      "corrector_compl_avrg_red_fact",
     
    210158      "This option determines the factor by which complementarity is allowed to increase "
    211159      "for a corrector step to be accepted.");
    212 
    213     roptions->AddStringOption2(
    214       "expect_infeasible_problem",
    215       "Enable heuristics to quickly detect an infeasible problem.",
    216       "no",
    217       "no", "the problem probably be feasible",
    218       "yes", "the problem has a good chance to be infeasible",
    219       "This options is meant to activate heuristics that may speed up the "
    220       "infeasibility determination if you expect that there is a good chance for the problem to be "
    221       "infeasible.  In the filter line search procedure, the restoration "
    222       "phase is called more qucikly than usually, and more reduction in "
    223       "the constraint violation is enforced. If the problem is square, this "
    224       "option is enabled automatically.");
    225     roptions->AddLowerBoundedNumberOption(
    226       "expect_infeasible_problem_ctol",
    227       "Threshold for disabling \"expect_infeasible_problem\" option.",
    228       0.0, false, 1e-3,
    229       "If the constraint violation becomes smaller than this threshold, "
    230       "the \"expect_infeasible_problem\" heuristics in the filter line "
    231       "search are disabled. If the problem is square, this options is set to "
    232       "0.");
    233     roptions->AddLowerBoundedNumberOption(
    234       "soft_resto_pderror_reduction_factor",
    235       "Required reduction in primal-dual error in the soft restoration phase.",
    236       0.0, false, (1.0 - 1e-4),
    237       "The soft restoration phase attempts to reduce the "
    238       "primal-dual error with regular steps. If the damped "
    239       "primal-dual step (damped only to satisfy the "
    240       "fraction-to-the-boundary rule) is not decreasing the primal-dual error "
    241       "by at least this factor, then the regular restoration phase is called. "
    242       "Choosing \"0\" here disables the soft "
    243       "restoration phase.");
    244     roptions->AddStringOption2(
    245       "start_with_resto",
    246       "Tells algorithm to switch to restoration phase in first iteration.",
    247       "no",
    248       "no", "don't force start in restoration phase",
    249       "yes", "force start in restoration phase",
    250       "Setting this option to \"yes\" forces the algorithm to switch to the "
    251       "feasibility restoration phase in the first iteration. If the initial "
    252       "point is feasible, the algorithm will abort with a failure.");
    253     roptions->AddLowerBoundedNumberOption(
    254       "tiny_step_tol",
    255       "Tolerance for detecting numerically insignificant steps.",
    256       0.0, false, 10.0*std::numeric_limits<double>::epsilon(),
    257       "If the search direction in the primal variables (x and s) is, in "
    258       "relative terms for each component, less than this value, the "
    259       "algorithm accepts the full step without line search.  If this happens "
    260       "repeatedly, the algorithm will terminate with a corresponding exit "
    261       "message. The default value is 10 times machine precision.");
    262     roptions->AddLowerBoundedIntegerOption(
    263       "watchdog_shortened_iter_trigger",
    264       "Number of shortened iterations that trigger the watchdog.",
    265       0, 10,
    266       "If the number of iterations in which the backtracking line search "
    267       "did not accept the first trial point exceedes this number, the "
    268       "watchdog procedure is activated.  Choosing \"0\" here disables the "
    269       "watchdog procedure.");
    270     roptions->AddLowerBoundedIntegerOption(
    271       "watchdog_trial_iter_max",
    272       "Maximum number of watchdog iterations.",
    273       1, 3,
    274       "This option determines the number of trial iterations "
    275       "allowed before the watchdog "
    276       "procedure is aborted and the algorithm returns to the stored point.");
    277   }
    278 
    279   bool FilterLineSearch::InitializeImpl(const OptionsList& options,
     160  }
     161
     162  bool FilterLSAcceptor::InitializeImpl(const OptionsList& options,
    280163                                        const std::string& prefix)
    281164  {
     
    291174    options.GetNumericValue("gamma_theta", gamma_theta_, prefix);
    292175    options.GetNumericValue("alpha_min_frac", alpha_min_frac_, prefix);
    293     options.GetNumericValue("alpha_red_factor", alpha_red_factor_, prefix);
    294176    options.GetIntegerValue("max_soc", max_soc_, prefix);
    295177    if (max_soc_>0) {
    296178      ASSERT_EXCEPTION(IsValid(pd_solver_), OPTION_INVALID,
    297                        "Option \"max_soc\": This option is non-negative, but no linear solver for computing the SOC given to FilterLineSearch object.");
     179                       "Option \"max_soc\": This option is non-negative, but no linear solver for computing the SOC given to FilterLSAcceptor object.");
    298180    }
    299181    options.GetNumericValue("kappa_soc", kappa_soc_, prefix);
    300182    options.GetNumericValue("obj_max_inc", obj_max_inc_, prefix);
    301     options.GetBoolValue("magic_steps", magic_steps_, prefix);
    302183    Index enum_int;
    303184    options.GetEnumValue("corrector_type", enum_int, prefix);
     
    305186    options.GetBoolValue("skip_corr_if_neg_curv", skip_corr_if_neg_curv_, prefix);
    306187    options.GetBoolValue("skip_corr_in_monotone_mode", skip_corr_in_monotone_mode_, prefix);
    307     options.GetBoolValue("accept_every_trial_step", accept_every_trial_step_, prefix);
    308     options.GetEnumValue("alpha_for_y", enum_int, prefix);
    309     alpha_for_y_ = AlphaForYEnum(enum_int);
    310188    options.GetNumericValue("corrector_compl_avrg_red_fact", corrector_compl_avrg_red_fact_, prefix);
    311     options.GetNumericValue("expect_infeasible_problem_ctol", expect_infeasible_problem_ctol_, prefix);
    312     options.GetBoolValue("expect_infeasible_problem", expect_infeasible_problem_, prefix);
    313 
    314     options.GetBoolValue("start_with_resto", start_with_resto_, prefix);
    315 
    316     bool retvalue = true;
    317     if (IsValid(resto_phase_)) {
    318       retvalue = resto_phase_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
    319                                           options, prefix);
    320     }
    321 
    322     options.GetNumericValue("soft_resto_pderror_reduction_factor",
    323                             soft_resto_pderror_reduction_factor_, prefix);
    324     options.GetNumericValue("tiny_step_tol", tiny_step_tol_, prefix);
    325     options.GetIntegerValue("watchdog_trial_iter_max", watchdog_trial_iter_max_, prefix);
    326     options.GetIntegerValue("watchdog_shortened_iter_trigger", watchdog_shortened_iter_trigger_, prefix);
    327 
    328     rigorous_ = true;
    329     skipped_line_search_ = false;
    330     tiny_step_last_iteration_ = false;
     189
    331190    theta_min_ = -1.;
    332191    theta_max_ = -1.;
     
    334193    Reset();
    335194
    336     count_successive_shortened_steps_ = 0;
    337 
    338     acceptable_iterate_ = NULL;
    339 
    340     return retvalue;
    341   }
    342 
    343   void FilterLineSearch::FindAcceptableTrialPoint()
    344   {
    345     DBG_START_METH("FilterLineSearch::FindAcceptableTrialPoint",
     195    return true;
     196  }
     197
     198  void FilterLSAcceptor::InitThisLineSearch(bool in_watchdog)
     199  {
     200    DBG_START_METH("FilterLSAcceptor::InitThisLineSearch",
    346201                   dbg_verbosity);
    347     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    348                    "--> Starting filter line search in iteration %d <--\n",
    349                    IpData().iter_count());
    350 
    351     // If the problem is square, we want to enable the
    352     // expect_infeasible_problem option automatically so that the
    353     // restoration phase is entered soon
    354     if (IpCq().IsSquareProblem()) {
    355       expect_infeasible_problem_ = true;
    356       expect_infeasible_problem_ctol_ = 0.;
    357     }
    358 
    359     // Store current iterate if the optimality error is on acceptable
    360     // level to restored if things fail later
    361     if (CurrentIsAcceptable()) {
    362       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    363                      "Storing current iterate as backup acceptable point.\n");
    364       StoreAcceptablePoint();
    365     }
    366 
    367     // First assume that line search will find an acceptable trial point
    368     skipped_line_search_ = false;
    369202
    370203    // Set the values for the reference point
    371     if (!in_watchdog_) {
     204    if (!in_watchdog) {
    372205      reference_theta_ = IpCq().curr_constraint_violation();
    373206      reference_barr_ = IpCq().curr_barrier_obj();
     
    379212      reference_gradBarrTDelta_ = watchdog_gradBarrTDelta_;
    380213    }
    381 
    382     // Get the search directions (this will store the actual search
    383     // direction, possibly including higher order corrections)
    384     SmartPtr<IteratesVector> actual_delta = IpData().delta()->MakeNewContainer();
    385 
    386     bool goto_resto = false;
    387     if (actual_delta->Asum() == 0. ) {
    388       // In this case, a search direction could not be computed, and
    389       // we should immediately go to the restoration phase.  ToDo: Cue
    390       // off of a something else than the norm of the search direction
    391       goto_resto = true;
    392     }
    393 
    394     if (start_with_resto_) {
    395       // If the use requested to start with the restoration phase,
    396       // skip the lin e search and do exactly that.  Reset the flag so
    397       // that this happens only once.
    398       goto_resto = true;
    399       start_with_resto_= false;
    400     }
    401 
    402     bool accept = false;
    403     bool corr_taken = false;
    404     bool soc_taken = false;
    405     Index n_steps = 0;
    406     Number alpha_primal = 0.;
    407 
    408     // Check if search direction becomes too small
    409     // ToDo: move this into place independent of this particular line search?
    410     bool tiny_step = (!goto_resto && DetectTinyStep());
    411 
    412     if (in_watchdog_ && (goto_resto || tiny_step)) {
    413       // If the step could not be computed or is too small and the
    414       // watchdog is active, stop the watch dog and resume everything
    415       // from reference point
    416       StopWatchDog(actual_delta);
    417       goto_resto = false;
    418       tiny_step = false;
    419     }
    420 
    421     // Check if we want to wake up the watchdog
    422     if (watchdog_shortened_iter_trigger_ > 0 &&
    423         !in_watchdog_ && !goto_resto && !tiny_step &&
    424         !in_soft_resto_phase_ && !expect_infeasible_problem_ &&
    425         watchdog_shortened_iter_ >= watchdog_shortened_iter_trigger_) {
    426       StartWatchDog();
    427     }
    428 
    429     if (tiny_step) {
    430       alpha_primal =
    431         IpCq().curr_primal_frac_to_the_bound(IpData().curr_tau());
    432       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    433                      "Tiny step detected. Use step size alpha = %e unchecked\n",
    434                      alpha_primal);
    435       IpData().SetTrialPrimalVariablesFromStep(alpha_primal, *IpData().delta()->x(), *IpData().delta()->s());
    436       IpData().Set_info_ls_count(0);
    437 
    438       if (tiny_step_last_iteration_) {
    439         IpData().Set_info_alpha_primal_char('T');
    440         IpData().Set_tiny_step_flag(true);
    441       }
    442       else {
    443         IpData().Set_info_alpha_primal_char('t');
    444       }
    445 
    446       tiny_step_last_iteration_ = true;
    447       accept = true;
    448     }
    449     else {
    450       tiny_step_last_iteration_ = false;
    451     }
    452 
    453     if (!goto_resto && !tiny_step) {
    454 
    455       if (in_soft_resto_phase_) {
    456         bool satisfies_original_filter = false;
    457         // ToDo use tiny_step in TrySoftRestoStep?
    458         accept = TrySoftRestoStep(actual_delta,
    459                                   satisfies_original_filter);
    460         if (accept) {
    461           IpData().Set_info_alpha_primal_char('s');
    462           if (satisfies_original_filter) {
    463             in_soft_resto_phase_ = false;
    464             IpData().Set_info_alpha_primal_char('S');
    465           }
    466         }
    467       }
    468       else {
    469         bool done = false;
    470         bool skip_first_trial_point = false;
    471         bool evaluation_error;
    472         while (!done) {
    473           accept = DoBacktrackingLineSearch(skip_first_trial_point,
    474                                             alpha_primal,
    475                                             corr_taken,
    476                                             soc_taken,
    477                                             n_steps,
    478                                             evaluation_error,
    479                                             actual_delta);
    480           DBG_PRINT((1, "evaluation_error = %d\n", evaluation_error));
    481           if (in_watchdog_) {
    482             if (accept) {
    483               in_watchdog_ = false;
    484               IpData().Append_info_string("W");
    485               Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    486                              "Watch dog procedure successful!\n");
    487               done = true;
    488             }
    489             else {
    490               watchdog_trial_iter_++;
    491               if (evaluation_error ||
    492                   watchdog_trial_iter_ > watchdog_trial_iter_max_) {
    493                 StopWatchDog(actual_delta);
    494                 skip_first_trial_point = true;
    495               }
    496               else {
    497                 done = true;
    498                 accept = true;
    499               }
    500             }
    501           }
    502           else {
    503             done = true;
    504           }
    505         }
    506       } /* else: if (in_soft_resto_phase_) { */
    507     } /* if (!goto_resto && !tiny_step) { */
    508 
    509     // If line search has been aborted because the step size becomes too small,
    510     // go to the restoration phase
    511     if (!accept) {
    512       // If we are not asked to do a rigorous line search, do no call
    513       // the restoration phase.
    514       if (!rigorous_) {
    515         Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Skipping call of restoration phase...\n");
    516         skipped_line_search_=true;
    517       }
    518       else {
    519         // Check if we should start the soft restoration phase
    520         if (!in_soft_resto_phase_ && soft_resto_pderror_reduction_factor_>0.
    521             && !goto_resto && !expect_infeasible_problem_) {
    522           Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    523                          "--> Starting soft restoration phase <--\n");
    524           // Augment the filter with the current point
    525           AugmentFilter();
    526 
    527           // Try the current search direction for the soft restoration phase
    528           bool satisfies_original_filter;
    529           accept = TrySoftRestoStep(actual_delta,
    530                                     satisfies_original_filter);
    531           // If it has been accepted: If the original filter is also
    532           // satisfied, we can just take that step and continue with
    533           // the regular algorithm, otherwise we stay in the soft
    534           // restoration phase
    535           if (accept) {
    536             if (satisfies_original_filter) {
    537               IpData().Set_info_alpha_primal_char('S');
    538             }
    539             else {
    540               in_soft_resto_phase_ = true;
    541               IpData().Set_info_alpha_primal_char('s');
    542             }
    543           }
    544         }
    545 
    546         if (!accept) {
    547           if (!in_soft_resto_phase_) {
    548             // Augment the filter with the current point if we are
    549             // already in the soft restoration phase, this has been
    550             // done earlier
    551             AugmentFilter();
    552           }
    553           if (CurrentIsAcceptable()) {
    554             THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
    555                             "Restoration phase called at acceptable point.");
    556           }
    557 
    558           if (!IsValid(resto_phase_)) {
    559             THROW_EXCEPTION(IpoptException, "No Restoration Phase given to this Filter Line Search Object!");
    560           }
    561           // ToDo make the 1e-2 below a parameter?
    562           if (IpCq().curr_constraint_violation()<=
    563               1e-2*IpData().tol()) {
    564             bool found_acceptable = RestoreAcceptablePoint();
    565             if (found_acceptable) {
    566               THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
    567                               "Restoration phase called at almost feasible point, but acceptable point could be restored.\n");
    568             }
    569             else {
    570               // ToDo does that happen too often?
    571               THROW_EXCEPTION(RESTORATION_FAILED,
    572                               "Restoration phase called, but point is almost feasible.");
    573             }
    574           }
    575 
    576           // Set the info fields for the first output line in the
    577           // restoration phase which reflects why the restoration phase
    578           // was called
    579           IpData().Set_info_alpha_primal(alpha_primal);
    580           IpData().Set_info_alpha_dual(0.);
    581           IpData().Set_info_alpha_primal_char('R');
    582           IpData().Set_info_ls_count(n_steps+1);
    583 
    584           accept = resto_phase_->PerformRestoration();
    585           if (!accept) {
    586             bool found_acceptable = RestoreAcceptablePoint();
    587             if (found_acceptable) {
    588               THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
    589                               "Restoration phase failed, but acceptable point could be restore.\n");
    590             }
    591             else {
    592               THROW_EXCEPTION(RESTORATION_FAILED,
    593                               "Failed restoration phase!!!");
    594             }
    595           }
    596           count_successive_shortened_steps_ = 0;
    597           if (expect_infeasible_problem_) {
    598             expect_infeasible_problem_ = false;
    599           }
    600           in_soft_resto_phase_ = false;
    601           watchdog_shortened_iter_ = 0;
    602         }
    603       }
    604     }
    605     else if (!in_soft_resto_phase_ || tiny_step) {
    606       // we didn't do the restoration phase and are now updating the
    607       // dual variables of the trial point
    608       Number alpha_dual_max =
    609         IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
    610                                       *actual_delta->z_L(), *actual_delta->z_U(),
    611                                       *actual_delta->v_L(), *actual_delta->v_U());
    612 
    613       PerformDualStep(alpha_primal, alpha_dual_max, actual_delta);
    614 
    615       if (n_steps==0) {
    616         // accepted this if a full step was
    617         // taken
    618         count_successive_shortened_steps_ = 0;
    619         watchdog_shortened_iter_ = 0;
    620       }
    621       else {
    622         count_successive_shortened_steps_++;
    623         watchdog_shortened_iter_++;
    624       }
    625 
    626       if (expect_infeasible_problem_ &&
    627           IpCq().curr_constraint_violation() <= expect_infeasible_problem_ctol_) {
    628         Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    629                        "Constraint violation is with %e less than expect_infeasible_problem_ctol.\nDisable expect_infeasible_problem_heuristic.\n", IpCq().curr_constraint_violation());
    630         expect_infeasible_problem_ = false;
    631       }
    632     }
    633   }
    634 
    635   bool FilterLineSearch::DoBacktrackingLineSearch(bool skip_first_trial_point,
    636       Number& alpha_primal,
    637       bool& corr_taken,
    638       bool& soc_taken,
    639       Index& n_steps,
    640       bool& evaluation_error,
    641       SmartPtr<IteratesVector>& actual_delta)
    642   {
    643     evaluation_error = false;
    644     bool accept = false;
    645 
    646     DBG_START_METH("FilterLineSearch::DoBacktrackingLineSearch",
    647                    dbg_verbosity);
    648 
    649     // Compute primal fraction-to-the-boundary value
    650     Number alpha_primal_max =
    651       IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
    652                                       *actual_delta->x(),
    653                                       *actual_delta->s());
    654 
    655     // Compute smallest step size allowed
    656     Number alpha_min = alpha_primal_max;
    657     if (!in_watchdog_) {
    658       alpha_min = CalculateAlphaMin();
    659     }
    660     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    661                    "minimal step size ALPHA_MIN = %E\n", alpha_min);
    662 
    663     // Start line search from maximal step size
    664     alpha_primal = alpha_primal_max;
    665 
    666     // Step size used in ftype and armijo tests
    667     Number alpha_primal_test = alpha_primal;
    668     if (in_watchdog_) {
    669       alpha_primal_test = watchdog_alpha_primal_test_;
    670     }
    671 
    672     if (skip_first_trial_point) {
    673       alpha_primal *= alpha_red_factor_;
    674     }
    675 
    676214    filter_.Print(Jnlst());
    677 
    678     if (corrector_type_!=NO_CORRECTOR && !skip_first_trial_point &&
    679         (!skip_corr_if_neg_curv_ || IpData().info_regu_x()==0.) &&
    680         (!skip_corr_in_monotone_mode_ || IpData().FreeMuMode()) ) {
    681       // Before we do the actual backtracking line search for the
    682       // regular primal-dual search direction, let's see if a step
    683       // including a higher-order correctior is already acceptable
    684       accept = TryCorrector(alpha_primal_test,
    685                             alpha_primal,
    686                             actual_delta);
    687     }
    688     if (accept) {
    689       corr_taken = true;
    690     }
    691 
    692     if (!accept) {
    693       // Loop over decreaseing step sizes until acceptable point is
    694       // found or until step size becomes too small
    695 
    696       while (alpha_primal>alpha_min ||
    697              n_steps == 0) { // always allow the "full" step if it is
    698         // acceptable (even if alpha_primal<=alpha_min)
    699         Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    700                        "Starting checks for alpha (primal) = %8.2e\n",
    701                        alpha_primal);
    702 
    703         try {
    704           // Compute the primal trial point
    705           IpData().SetTrialPrimalVariablesFromStep(alpha_primal, *actual_delta->x(), *actual_delta->s());
    706 
    707           if (magic_steps_) {
    708             PerformMagicStep();
    709           }
    710 
    711           // If it is acceptable, stop the search
    712           alpha_primal_test = alpha_primal;
    713           accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test);
    714         }
    715         catch(IpoptNLP::Eval_Error& e) {
    716           e.ReportException(Jnlst(), J_DETAILED);
    717           Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
    718                          "Warning: Cutting back alpha due to evaluation error\n");
    719           IpData().Append_info_string("e");
    720           accept = false;
    721           evaluation_error = true;
    722         }
    723 
    724         if (accept) {
    725           break;
    726         }
    727 
    728         if (in_watchdog_) {
    729           break;
    730         }
    731 
    732         // Decide if we want to go to the restoration phase in a
    733         // short cut to check if the problem is infeasible
    734         if (expect_infeasible_problem_) {
    735           if (count_successive_shortened_steps_>=5) {
    736             break;
    737           }
    738         }
    739 
    740         // try second order correction step if the function could
    741         // be evaluated
    742         // DoTo: check if we want to do SOC when watchdog is active
    743         if (!evaluation_error) {
    744           Number theta_curr = IpCq().curr_constraint_violation();
    745           Number theta_trial = IpCq().trial_constraint_violation();
    746           if (alpha_primal==alpha_primal_max &&       // i.e. first trial point
    747               theta_curr<=theta_trial && max_soc_>0) {
    748             // Try second order correction
    749             accept = TrySecondOrderCorrection(alpha_primal_test,
    750                                               alpha_primal,
    751                                               actual_delta);
    752           }
    753           if (accept) {
    754             soc_taken = true;
    755             break;
    756           }
    757         }
    758 
    759         // Point is not yet acceptable, try a shorter one
    760         alpha_primal *= alpha_red_factor_;
    761         n_steps++;
    762       }
    763     } /* if (!accept) */
    764 
    765     char info_alpha_primal_char;
    766     if (!accept && in_watchdog_) {
    767       info_alpha_primal_char = 'w';
    768     }
    769     else {
    770       // Augment the filter if required
    771       if (!IsFtype(alpha_primal_test) ||
    772           !ArmijoHolds(alpha_primal_test)) {
    773         AugmentFilter();
    774         info_alpha_primal_char = 'h';
    775       }
    776       else {
    777         info_alpha_primal_char = 'f';
    778       }
    779     }
    780     if (soc_taken) {
    781       info_alpha_primal_char = toupper(info_alpha_primal_char);
    782     }
    783     IpData().Set_info_alpha_primal_char(info_alpha_primal_char);
    784     IpData().Set_info_ls_count(n_steps+1);
    785     if (corr_taken) {
    786       IpData().Append_info_string("C");
    787     }
    788 
    789     return accept;
    790   }
    791 
    792   bool FilterLineSearch::IsFtype(Number alpha_primal_test)
    793   {
    794     DBG_START_METH("FilterLineSearch::IsFtype",
     215  }
     216
     217  bool FilterLSAcceptor::IsFtype(Number alpha_primal_test)
     218  {
     219    DBG_START_METH("FilterLSAcceptor::IsFtype",
    795220                   dbg_verbosity);
    796221    DBG_ASSERT(reference_theta_>0. || reference_gradBarrTDelta_ < 0.0);
     
    800225  }
    801226
    802   void FilterLineSearch::AugmentFilter()
    803   {
    804     DBG_START_METH("FilterLineSearch::AugmentFilter",
     227  void FilterLSAcceptor::AugmentFilter()
     228  {
     229    DBG_START_METH("FilterLSAcceptor::AugmentFilter",
    805230                   dbg_verbosity);
    806231
     
    812237
    813238  bool
    814   FilterLineSearch::CheckAcceptabilityOfTrialPoint(Number alpha_primal_test)
    815   {
    816     DBG_START_METH("FilterLineSearch::CheckAcceptabilityOfTrialPoint",
     239  FilterLSAcceptor::CheckAcceptabilityOfTrialPoint(Number alpha_primal_test)
     240  {
     241    DBG_START_METH("FilterLSAcceptor::CheckAcceptabilityOfTrialPoint",
    817242                   dbg_verbosity);
    818 
    819     if (accept_every_trial_step_) {
    820 
    821       // We call the evaluation at the trial point here, so that an
    822       // exception will the thrown if there are problem during the
    823       // evaluation of the functions (in that case, we want to further
    824       // reduce the step size
    825       /* Number trial_barr = */ IpCq().trial_barrier_obj();
    826       /* Number trial_theta = */
    827       IpCq().trial_constraint_violation();
    828       return true;
    829     }
    830243
    831244    bool accept;
     
    862275
    863276    // Check if point is acceptable w.r.t current iterate
    864     if (IsFtype(alpha_primal_test) && reference_theta_ <= theta_min_) {
     277    if (alpha_primal_test>0. && IsFtype(alpha_primal_test) &&
     278        reference_theta_ <= theta_min_) {
    865279      // Armijo condition for the barrier function has to be satisfied
    866280      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Checking Armijo Condition...\n");
     
    894308  }
    895309
    896   bool FilterLineSearch::ArmijoHolds(Number alpha_primal_test)
     310  bool FilterLSAcceptor::ArmijoHolds(Number alpha_primal_test)
    897311  {
    898312    /*
     
    905319  }
    906320
    907   Number FilterLineSearch::CalculateAlphaMin()
     321  Number FilterLSAcceptor::CalculateAlphaMin()
    908322  {
    909323    Number gBD = IpCq().curr_gradBarrTDelta();
     
    924338  }
    925339
    926   bool FilterLineSearch::IsAcceptableToCurrentIterate(Number trial_barr,
     340  bool FilterLSAcceptor::IsAcceptableToCurrentIterate(Number trial_barr,
    927341      Number trial_theta,
    928342      bool called_from_restoration /*=false*/) const
    929343  {
    930     DBG_START_METH("FilterLineSearch::IsAcceptableToCurrentIterate",
     344    DBG_START_METH("FilterLSAcceptor::IsAcceptableToCurrentIterate",
    931345                   dbg_verbosity);
    932346
     
    950364  }
    951365
    952   bool FilterLineSearch::IsAcceptableToCurrentFilter(Number trial_barr, Number trial_theta) const
     366  bool FilterLSAcceptor::IsAcceptableToCurrentFilter(Number trial_barr, Number trial_theta) const
    953367  {
    954368    return filter_.Acceptable(trial_barr, trial_theta);
    955369  }
    956370
    957   bool FilterLineSearch::Compare_le(Number lhs, Number rhs, Number BasVal)
    958   {
    959     DBG_START_FUN("FilterLineSearch::Compare_le",
     371  bool FilterLSAcceptor::Compare_le(Number lhs, Number rhs, Number BasVal)
     372  {
     373    DBG_START_FUN("FilterLSAcceptor::Compare_le",
    960374                  dbg_verbosity);
    961375    DBG_PRINT((1,"lhs = %27.16e rhs = %27.16e  BasVal = %27.16e\n",lhs,rhs,BasVal));
    962     // ToDo: Comparison based on machine precision
    963     return (lhs - rhs <= 1e-15*fabs(BasVal));
    964   }
    965 
    966   bool FilterLineSearch::TrySoftRestoStep(SmartPtr<IteratesVector>& actual_delta,
    967                                           bool &satisfies_original_filter)
    968   {
    969     DBG_START_FUN("FilterLineSearch::TrySoftRestoStep", dbg_verbosity);
    970 
    971     satisfies_original_filter = false;
    972 
    973     // ToDo: Need to decide if we want to try a corrector step first
    974 
    975     // Compute the maximal step sizes (we use identical step sizes for
    976     // primal and dual variables
    977     Number alpha_primal_max =
    978       IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
    979                                       *actual_delta->x(),
    980                                       *actual_delta->s());
    981     Number alpha_dual_max =
    982       IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
    983                                     *actual_delta->z_L(),
    984                                     *actual_delta->z_U(),
    985                                     *actual_delta->v_L(),
    986                                     *actual_delta->v_U());
    987     Number alpha_max =  Min(alpha_primal_max, alpha_dual_max);
    988 
    989     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    990                    "Trying soft restoration phase step with step length %13.6e\n",
    991                    alpha_max);
    992 
    993     // Set the trial point
    994     IpData().SetTrialPrimalVariablesFromStep(alpha_max, *actual_delta->x(), *actual_delta->s());
    995     PerformDualStep(alpha_max, alpha_max,
    996                     actual_delta);
    997 
    998     // Check if that point is acceptable with respect to the current
    999     // original filter
    1000 
    1001     Number trial_barr;
    1002     Number trial_theta;
    1003     try {
    1004       trial_barr = IpCq().trial_barrier_obj();
    1005       trial_theta = IpCq().trial_constraint_violation();
    1006     }
    1007     catch(IpoptNLP::Eval_Error& e) {
    1008       e.ReportException(Jnlst(), J_DETAILED);
    1009       Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
    1010                      "Warning: Evaluation error during soft restoration phase step.\n");
    1011       IpData().Append_info_string("e");
    1012       return false;
    1013     }
    1014     if (theta_max_<=0 || trial_theta<=theta_max_) {
    1015       if (IsAcceptableToCurrentIterate(trial_barr, trial_theta)) {
    1016         Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1017                        "  Trial step acceptable with respect to original filter.\n");
    1018         satisfies_original_filter = true;
    1019         return true;
    1020       }
    1021     }
    1022 
    1023     // Evaluate the optimality error at the new point
    1024     Number mu = .0;
    1025     if (!IpData().FreeMuMode()) {
    1026       mu = IpData().curr_mu();
    1027     }
    1028     Number trial_pderror;
    1029     Number curr_pderror;
    1030     try {
    1031       trial_pderror = IpCq().trial_primal_dual_system_error(mu);
    1032       curr_pderror = IpCq().curr_primal_dual_system_error(mu);
    1033     }
    1034     catch(IpoptNLP::Eval_Error& e) {
    1035       e.ReportException(Jnlst(), J_DETAILED);
    1036       Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
    1037                      "Warning: Evaluation error during soft restoration phase step.\n");
    1038       IpData().Append_info_string("e");
    1039       return false;
    1040     }
    1041 
    1042     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1043                    "  Primal-dual error at current point:  %23.16e\n", curr_pderror);
    1044     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1045                    "  Primal-dual error at trial point  :  %23.16e\n", trial_pderror);
    1046     // Check if there is sufficient reduction in the optimality error
    1047     if (trial_pderror <= soft_resto_pderror_reduction_factor_*curr_pderror) {
    1048       Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1049                      "  Trial step accepted.\n");
    1050       return true;
    1051     }
    1052 
    1053     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1054                    "  Trial step rejected.\n");
    1055     return false;
    1056   }
    1057 
    1058   void FilterLineSearch::StartWatchDog()
    1059   {
    1060     DBG_START_FUN("FilterLineSearch::StartWatchDog", dbg_verbosity);
    1061 
    1062     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1063                    "Starting Watch Dog\n");
    1064 
    1065     in_watchdog_ = true;
    1066     watchdog_iterate_ = IpData().curr();
    1067     watchdog_delta_ = IpData().delta();
    1068     watchdog_trial_iter_ = 0;
    1069     watchdog_alpha_primal_test_ =
    1070       IpCq().curr_primal_frac_to_the_bound(IpData().curr_tau());
     376
     377    Number mach_eps = std::numeric_limits<Number>::epsilon();
     378    return (lhs - rhs <= 1e10*mach_eps*fabs(BasVal));
     379  }
     380
     381  void FilterLSAcceptor::StartWatchDog()
     382  {
     383    DBG_START_FUN("FilterLSAcceptor::StartWatchDog", dbg_verbosity);
     384
    1071385    watchdog_theta_ = IpCq().curr_constraint_violation();
    1072386    watchdog_barr_ = IpCq().curr_barrier_obj();
     
    1074388  }
    1075389
    1076   void FilterLineSearch::StopWatchDog(SmartPtr<IteratesVector>& actual_delta)
    1077   {
    1078     DBG_START_FUN("FilterLineSearch::StopWatchDog", dbg_verbosity);
    1079 
    1080     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1081                    "Stopping Watch Dog\n");
    1082 
    1083     IpData().Append_info_string("w");
    1084 
    1085     in_watchdog_ = false;
    1086 
    1087     // Reset all fields in IpData to reference point
    1088     SmartPtr<IteratesVector> old_trial = watchdog_iterate_->MakeNewContainer();
    1089     IpData().set_trial(old_trial);
    1090     IpData().AcceptTrialPoint();
    1091     actual_delta = watchdog_delta_->MakeNewContainer();
    1092     IpData().SetHaveAffineDeltas(false);
    1093 
    1094     // reset the stored watchdog iterates
    1095     watchdog_iterate_ = NULL;
    1096     watchdog_delta_ = NULL;
    1097 
    1098     watchdog_shortened_iter_ = 0;
     390  void FilterLSAcceptor::StopWatchDog()
     391  {
     392    DBG_START_FUN("FilterLSAcceptor::StopWatchDog", dbg_verbosity);
    1099393
    1100394    reference_theta_ = watchdog_theta_;
     
    1103397  }
    1104398
    1105   void FilterLineSearch::Reset()
    1106   {
    1107     DBG_START_FUN("FilterLineSearch::Reset", dbg_verbosity);
    1108     in_soft_resto_phase_ = false;
    1109 
    1110     // Inactivate the watchdog and release all stored data
    1111     in_watchdog_ = false;
    1112     watchdog_iterate_ = NULL;
    1113     watchdog_delta_ = NULL;
    1114     watchdog_shortened_iter_ = 0;
     399  void FilterLSAcceptor::Reset()
     400  {
     401    DBG_START_FUN("FilterLSAcceptor::Reset", dbg_verbosity);
    1115402
    1116403    filter_.Clear();
    1117404  }
    1118405
    1119   void FilterLineSearch::PerformDualStep(Number alpha_primal,
    1120                                          Number alpha_dual,
    1121                                          SmartPtr<IteratesVector>& delta)
    1122   {
    1123     DBG_START_FUN("FilterLineSearch::PerformDualStep", dbg_verbosity);
    1124 
    1125     // set the bound multipliers from the step
    1126     IpData().SetTrialBoundMultipliersFromStep(alpha_dual, *delta->z_L(), *delta->z_U(), *delta->v_L(), *delta->v_U());
    1127 
    1128     Number alpha_y=-1.;
    1129     switch (alpha_for_y_) {
    1130       case PRIMAL_ALPHA_FOR_Y:
    1131       alpha_y = alpha_primal;
    1132       break;
    1133       case DUAL_ALPHA_FOR_Y:
    1134       alpha_y = alpha_dual;
    1135       break;
    1136       case MIN_ALPHA_FOR_Y:
    1137       alpha_y = Min(alpha_dual, alpha_primal);
    1138       break;
    1139       case MAX_ALPHA_FOR_Y:
    1140       alpha_y = Max(alpha_dual, alpha_primal);
    1141       break;
    1142       case FULL_STEP_FOR_Y:
    1143       alpha_y = 1;
    1144       break;
    1145       case MIN_DUAL_INFEAS_ALPHA_FOR_Y:
    1146       case SAFE_MIN_DUAL_INFEAS_ALPHA_FOR_Y:
    1147       // Here we compute the step size for y so that the dual
    1148       // infeasibility is minimized along delta_y
    1149 
    1150       // compute the dual infeasibility at new point with old y
    1151       SmartPtr<IteratesVector> temp_trial
    1152       = IpData().trial()->MakeNewContainer();
    1153       temp_trial->Set_y_c(*IpData().curr()->y_c());
    1154       temp_trial->Set_y_d(*IpData().curr()->y_d());
    1155       IpData().set_trial(temp_trial);
    1156       SmartPtr<const Vector> dual_inf_x = IpCq().trial_grad_lag_x();
    1157       SmartPtr<const Vector> dual_inf_s = IpCq().trial_grad_lag_s();
    1158 
    1159       SmartPtr<Vector> new_jac_times_delta_y =
    1160         IpData().curr()->x()->MakeNew();
    1161       new_jac_times_delta_y->AddTwoVectors(1., *IpCq().trial_jac_cT_times_vec(*delta->y_c()),
    1162                                            1., *IpCq().trial_jac_dT_times_vec(*delta->y_d()),
    1163                                            0.);
    1164 
    1165       Number a = pow(new_jac_times_delta_y->Nrm2(), 2.) +
    1166                  pow(delta->y_d()->Nrm2(), 2.);
    1167       Number b = dual_inf_x->Dot(*new_jac_times_delta_y) -
    1168                  dual_inf_s->Dot(*delta->y_d());
    1169 
    1170       Number alpha = - b/a;
    1171 
    1172       if (alpha_for_y_==SAFE_MIN_DUAL_INFEAS_ALPHA_FOR_Y) {
    1173         alpha_y = Min(Max(alpha_primal, alpha_dual), Max(alpha, Min(alpha_primal, alpha_dual)));
    1174       }
    1175       else {
    1176         alpha_y = Min(1., Max(0., alpha));
    1177       }
    1178       break;
    1179     }
    1180 
    1181     // Set the eq multipliers from the step now that alpha_y
    1182     // has been calculated.
    1183     DBG_PRINT((1, "alpha_y = %e\n", alpha_y));
    1184     DBG_PRINT_VECTOR(2, "delta_y_c", *delta->y_c());
    1185     DBG_PRINT_VECTOR(2, "delta_y_d", *delta->y_d());
    1186     IpData().SetTrialEqMultipliersFromStep(alpha_y, *delta->y_c(), *delta->y_d());
    1187 
    1188     // Set some information for iteration summary output
    1189     IpData().Set_info_alpha_primal(alpha_primal);
    1190     IpData().Set_info_alpha_dual(alpha_dual);
    1191   }
    1192 
    1193406  bool
    1194   FilterLineSearch::TrySecondOrderCorrection(
     407  FilterLSAcceptor::TrySecondOrderCorrection(
    1195408    Number alpha_primal_test,
    1196409    Number& alpha_primal,
    1197410    SmartPtr<IteratesVector>& actual_delta)
    1198411  {
    1199     DBG_START_METH("FilterLineSearch::TrySecondOrderCorrection",
     412    DBG_START_METH("FilterLSAcceptor::TrySecondOrderCorrection",
    1200413                   dbg_verbosity);
     414
     415    if (max_soc_==0) {
     416      return false;
     417    }
    1201418
    1202419    bool accept = false;
     
    1274491
    1275492  bool
    1276   FilterLineSearch::TryCorrector(
     493  FilterLSAcceptor::TryCorrector(
    1277494    Number alpha_primal_test,
    1278495    Number& alpha_primal,
    1279496    SmartPtr<IteratesVector>& actual_delta)
    1280497  {
    1281     DBG_START_METH("FilterLineSearch::TryCorrector",
     498    if (corrector_type_==NO_CORRECTOR ||
     499        (skip_corr_if_neg_curv_ && IpData().info_regu_x()!=0.) ||
     500        (skip_corr_in_monotone_mode_ && !IpData().FreeMuMode())) {
     501      return false;
     502    }
     503
     504    DBG_START_METH("FilterLSAcceptor::TryCorrector",
    1282505                   dbg_verbosity);
    1283506
     
    1478701  }
    1479702
    1480   void
    1481   FilterLineSearch::PerformMagicStep()
    1482   {
    1483     DBG_START_METH("FilterLineSearch::PerformMagicStep",
    1484                    2);//dbg_verbosity);
    1485 
    1486     DBG_PRINT((1,"Incoming barr = %e and constrviol %e\n",
    1487                IpCq().trial_barrier_obj(),
    1488                IpCq().trial_constraint_violation()));
    1489     DBG_PRINT_VECTOR(2, "s in", *IpData().trial()->s());
    1490     DBG_PRINT_VECTOR(2, "d minus s in", *IpCq().trial_d_minus_s());
    1491     DBG_PRINT_VECTOR(2, "slack_s_L in", *IpCq().trial_slack_s_L());
    1492     DBG_PRINT_VECTOR(2, "slack_s_U in", *IpCq().trial_slack_s_U());
    1493 
    1494     SmartPtr<const Vector> d_L = IpNLP().d_L();
    1495     SmartPtr<const Matrix> Pd_L = IpNLP().Pd_L();
    1496     SmartPtr<Vector> delta_s_magic_L = d_L->MakeNew();
    1497     delta_s_magic_L->Set(0.);
    1498     SmartPtr<Vector> tmp = d_L->MakeNew();
    1499     Pd_L->TransMultVector(1., *IpCq().trial_d_minus_s(), 0., *tmp);
    1500     delta_s_magic_L->ElementWiseMax(*tmp);
    1501 
    1502     SmartPtr<const Vector> d_U = IpNLP().d_U();
    1503     SmartPtr<const Matrix> Pd_U = IpNLP().Pd_U();
    1504     SmartPtr<Vector> delta_s_magic_U = d_U->MakeNew();
    1505     delta_s_magic_U->Set(0.);
    1506     tmp = d_U->MakeNew();
    1507     Pd_U->TransMultVector(1., *IpCq().trial_d_minus_s(), 0., *tmp);
    1508     delta_s_magic_U->ElementWiseMin(*tmp);
    1509 
    1510     SmartPtr<Vector> delta_s_magic = IpData().trial()->s()->MakeNew();
    1511     Pd_L->MultVector(1., *delta_s_magic_L, 0., *delta_s_magic);
    1512     Pd_U->MultVector(1., *delta_s_magic_U, 1., *delta_s_magic);
    1513     delta_s_magic_L = NULL; // free memory
    1514     delta_s_magic_U = NULL; // free memory
    1515 
    1516     // Now find those entries with both lower and upper bounds, there
    1517     // the step is too large
    1518     // ToDo this should only be done if there are inequality
    1519     // constraints with two bounds
    1520     // also this can be done in a smaller space (d_L or d_U whichever
    1521     // is smaller)
    1522     tmp = delta_s_magic->MakeNew();
    1523     tmp->Copy(*IpData().trial()->s());
    1524     Pd_L->MultVector(1., *d_L, -2., *tmp);
    1525     Pd_U->MultVector(1., *d_U, 1., *tmp);
    1526     SmartPtr<Vector> tmp2 = tmp->MakeNew();
    1527     tmp2->Copy(*tmp);
    1528     tmp2->ElementWiseAbs();
    1529     tmp->Axpy(-2., *delta_s_magic);
    1530     tmp->ElementWiseAbs();
    1531     // now, tmp2 = |d_L + d_u - 2*s| and tmp = |d_L + d_u - 2*(s+Delta s)|
    1532     // we want to throw out those for which tmp2 > tmp
    1533     tmp->Axpy(-1., *tmp2);
    1534     tmp->ElementWiseSgn();
    1535     tmp2->Set(0.);
    1536     tmp2->ElementWiseMax(*tmp);
    1537     tmp = d_L->MakeNew();
    1538     Pd_L->TransMultVector(1., *tmp2, 0., *tmp);
    1539     Pd_L->MultVector(1., *tmp, 0., *tmp2);
    1540     tmp = d_U->MakeNew();
    1541     Pd_U->TransMultVector(1., *tmp2, 0., *tmp);
    1542     Pd_U->MultVector(1., *tmp, 0., *tmp2);
    1543     DBG_PRINT_VECTOR(2, "tmp indicator", *tmp2)
    1544     // tmp2 now is one for those entries with both bounds, for which
    1545     // no step should be taken
    1546 
    1547     tmp = delta_s_magic->MakeNew();
    1548     tmp->Copy(*delta_s_magic);
    1549     tmp->ElementWiseMultiply(*tmp2);
    1550     delta_s_magic->Axpy(-1., *tmp);
    1551 
    1552     Number delta_s_magic_max = delta_s_magic->Amax();
    1553     Number mach_eps = std::numeric_limits<Number>::epsilon();
    1554     if (delta_s_magic_max>0.) {
    1555       if (delta_s_magic_max > 10*mach_eps*IpData().trial()->s()->Amax()) {
    1556         IpData().Append_info_string("M");
    1557         Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Magic step with max-norm %.6e taken.\n", delta_s_magic->Amax());
    1558         delta_s_magic->Print(Jnlst(), J_MOREVECTOR, J_LINE_SEARCH,
    1559                              "delta_s_magic");
    1560       }
    1561 
    1562       // now finally compute the new overall slacks
    1563       delta_s_magic->Axpy(1., *IpData().trial()->s());
    1564       SmartPtr<IteratesVector> trial = IpData().trial()->MakeNewContainer();
    1565       trial->Set_s(*delta_s_magic);
    1566 
    1567       // also update the set in the dual variables
    1568 
    1569 
    1570       IpData().set_trial(trial);
    1571     }
    1572 
    1573     DBG_PRINT((1,"Outgoing barr = %e and constrviol %e\n", IpCq().trial_barrier_obj(), IpCq().trial_constraint_violation()));
    1574     DBG_PRINT_VECTOR(2, "s out", *IpData().trial()->s());
    1575     DBG_PRINT_VECTOR(2, "d minus s out", *IpCq().trial_d_minus_s());
    1576     DBG_PRINT_VECTOR(2, "slack_s_L out", *IpCq().trial_slack_s_L());
    1577     DBG_PRINT_VECTOR(2, "slack_s_U out", *IpCq().trial_slack_s_U());
    1578   }
    1579 
    1580   bool
    1581   FilterLineSearch::DetectTinyStep()
    1582   {
    1583     DBG_START_METH("FilterLineSearch::DetectTinyStep",
    1584                    dbg_verbosity);
    1585 
    1586     Number max_step_x;
    1587     Number max_step_s;
    1588 
    1589     if (tiny_step_tol_==0.)
    1590       return false;
    1591 
    1592     // ToDo try to find more efficient implementation
    1593 
    1594     SmartPtr<Vector> tmp = IpData().curr()->x()->MakeNew();
    1595     tmp->Copy(*IpData().curr()->x());
    1596     tmp->ElementWiseAbs();
    1597     tmp->AddScalar(1.);
    1598 
    1599     SmartPtr<Vector> tmp2 = IpData().curr()->x()->MakeNew();
    1600     tmp2->Copy(*IpData().delta()->x());
    1601     tmp2->ElementWiseDivide(*tmp);
    1602     max_step_x = tmp2->Amax();
    1603     Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
    1604                    "Relative step size for delta_x = %e\n",
    1605                    max_step_x);
    1606     if (max_step_x > tiny_step_tol_)
    1607       return false;
    1608 
    1609     tmp = IpData().curr()->s()->MakeNew();
    1610     tmp->Copy(*IpData().curr()->s());
    1611     tmp->ElementWiseAbs();
    1612     tmp->AddScalar(1.);
    1613 
    1614     tmp2 = IpData().curr()->s()->MakeNew();
    1615     tmp2->Copy(*IpData().delta()->s());
    1616     tmp2->ElementWiseDivide(*tmp);
    1617     max_step_s = tmp2->Amax();
    1618     Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
    1619                    "Relative step size for delta_s = %e\n",
    1620                    max_step_s);
    1621     if (max_step_s > tiny_step_tol_)
    1622       return false;
    1623 
    1624     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    1625                    "Tiny step of relative size %e detected.\n",
    1626                    Max(max_step_x, max_step_s));
    1627 
    1628     return true;
    1629   }
    1630 
    1631   bool FilterLineSearch::CurrentIsAcceptable()
    1632   {
    1633     return (IsValid(conv_check_) &&
    1634             conv_check_->CurrentIsAcceptable());
    1635   }
    1636 
    1637   void FilterLineSearch::StoreAcceptablePoint()
    1638   {
    1639     DBG_START_METH("FilterLineSearch::StoreAcceptablePoint",
    1640                    dbg_verbosity);
    1641 
    1642     acceptable_iterate_ = IpData().curr();
    1643   }
    1644 
    1645   bool FilterLineSearch::RestoreAcceptablePoint()
    1646   {
    1647     DBG_START_METH("FilterLineSearch::RestoreAcceptablePoint",
    1648                    dbg_verbosity);
    1649 
    1650     if (!IsValid(acceptable_iterate_)) {
    1651       return false;
    1652     }
    1653 
    1654     SmartPtr<IteratesVector> prev_iterate = acceptable_iterate_->MakeNewContainer();
    1655     IpData().set_trial(prev_iterate);
    1656     IpData().AcceptTrialPoint();
    1657 
    1658     return true;
     703  char FilterLSAcceptor::UpdateForNextIteration(Number alpha_primal_test)
     704  {
     705    char info_alpha_primal_char;
     706    // Augment the filter if required
     707    if (!IsFtype(alpha_primal_test) ||
     708        !ArmijoHolds(alpha_primal_test)) {
     709      AugmentFilter();
     710      info_alpha_primal_char = 'h';
     711    }
     712    else {
     713      info_alpha_primal_char = 'f';
     714    }
     715    return info_alpha_primal_char;
     716  }
     717
     718  void FilterLSAcceptor::PrepareRestoPhaseStart()
     719  {
     720    AugmentFilter();
    1659721  }
    1660722
  • branches/dev/Algorithm/IpFilterLSAcceptor.hpp

    r541 r542  
    1 // Copyright (C) 2004,2005 International Business Machines and others.
     1// Copyright (C) 2005 International Business Machines and others.
    22// All Rights Reserved.
    33// This code is published under the Common Public License.
     
    55// $Id$
    66//
    7 // Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
    8 
    9 #ifndef __IPFILTERLINESEARCH_HPP__
    10 #define __IPFILTERLINESEARCH_HPP__
     7// Authors:  Andreas Waechter                 IBM    2005-10-13
     8//               derived file from IpFilterLineSearch.hpp
     9
     10#ifndef __IPFILTERLSACCEPTOR_HPP__
     11#define __IPFILTERLSACCEPTOR_HPP__
    1112
    1213#include "IpFilter.hpp"
    13 #include "IpLineSearch.hpp"
    14 #include "IpRestoPhase.hpp"
     14#include "IpBacktrackingLSAcceptor.hpp"
    1515#include "IpPDSystemSolver.hpp"
    16 #include "IpConvCheck.hpp"
    1716
    1817namespace Ipopt
     
    2221   *  search procedure.
    2322   */
    24   class FilterLineSearch : public LineSearch
     23  class FilterLSAcceptor : public BacktrackingLSAcceptor
    2524  {
    2625  public:
     
    2827    //@{
    2928    /** Constructor.  The PDSystemSolver object only needs to be
    30      *  provided (i.e. not NULL) if second order correction is to be
    31      *  used.  The ConvergenceCheck object is used to determine
    32      *  whether the current iterate is acceptable (for example, the
    33      *  restoration phase is not started if the acceptability level
    34      *  has been reached).  If conv_check is NULL, we assume that the
    35      *  current iterate is not acceptable. */
    36     FilterLineSearch(const SmartPtr<RestorationPhase>& resto_phase,
    37                      const SmartPtr<PDSystemSolver>& pd_solver,
    38                      const SmartPtr<ConvergenceCheck>& conv_check
    39                     );
     29     *  provided (i.e. not NULL) if second order correction or
     30     *  corrector steps are to be used. */
     31    FilterLSAcceptor(const SmartPtr<PDSystemSolver>& pd_solver);
    4032
    4133    /** Default destructor */
    42     virtual ~FilterLineSearch();
     34    virtual ~FilterLSAcceptor();
    4335    //@}
    4436
     
    4739                                const std::string& prefix);
    4840
    49     /** Perform the line search.  It is assumed that the search
    50      *  direction is computed in the data object.
    51      */
    52     virtual void FindAcceptableTrialPoint();
    53 
    54     /** Reset the line search.
     41    /** Reset the acceptor.
    5542     *  This function should be called if all previous information
    5643     *  should be discarded when the line search is performed the
     
    6047    virtual void Reset();
    6148
    62     /** Set flag indicating whether a very rigorous line search should
    63      *  be performed.  If this flag is set to true, the line search
    64      *  algorithm might decide to abort the line search and not to
    65      *  accept a new iterate.  If the line search decided not to
    66      *  accept a new iterate, the return value of
    67      *  CheckSkippedLineSearch() is true at the next call.  For
    68      *  example, in the non-monotone barrier parameter update
    69      *  procedure, the filter algorithm should not switch to the
    70      *  restoration phase in the free mode; instead, the algorithm
    71      *  should swtich to the fixed mode.
    72      */
    73     virtual void SetRigorousLineSearch(bool rigorous)
    74     {
    75       rigorous_ = rigorous;
    76     }
    77 
    78     /** Check if the line search procedure didn't accept a new iterate
    79      *  during the last call of FindAcceptableTrialPoint().
    80      * 
    81      */
    82     virtual bool CheckSkippedLineSearch()
    83     {
    84       return skipped_line_search_;
    85     }
     49    /** Initialization for the next line search.  The flag in_watchdog
     50     *  indicates if we are currently in an active watchdog
     51     *  procedure. */
     52    virtual void InitThisLineSearch(bool in_watchdog);
     53
     54    /** Method that is called before the restoration phase is called.
     55     *  Here, we can set up things that are required in the
     56     *  termination test for the restoration phase, such as augmenting
     57     *  a filter. */
     58    virtual void PrepareRestoPhaseStart();
     59
     60    /** Method returning the lower bound on the trial step sizes.  If
     61     *  the backtracking procedure encounters a trial step size below
     62     *  this value after the first trial set, it swtiches to the
     63     *  (soft) restoration phase. */
     64    virtual Number CalculateAlphaMin();
     65
     66    /** Method for checking if current trial point is acceptable.
     67     *  It is assumed that the delta information in ip_data is the
     68     *  search direction used in criteria.  The primal trial point has
     69     *  to be set before the call.
     70     */
     71    virtual bool CheckAcceptabilityOfTrialPoint(Number alpha_primal);
     72
     73    /** Try a second order correction for the constraints.  If the
     74     *  first trial step (with incoming alpha_primal) has been reject,
     75     *  this tries up to max_soc_ second order corrections for the
     76     *  constraints.  Here, alpha_primal_test is the step size that
     77     *  has to be used in the filter acceptance tests.  On output
     78     *  actual_delta_ has been set to the step including the
     79     *  second order correction if it has been accepted, otherwise it
     80     *  is unchanged.  If the SOC step has been accepted, alpha_primal
     81     *  has the fraction-to-the-boundary value for the SOC step on output.
     82     *  The return value is true, if a SOC step has been accepted.
     83     */
     84    virtual bool TrySecondOrderCorrection(Number alpha_primal_test,
     85                                          Number& alpha_primal,
     86                                          SmartPtr<IteratesVector>& actual_delta);
     87
     88    /** Try higher order corrector (for fast local convergence).  In
     89     *  contrast to a second order correction step, which tries to
     90     *  make an unacceptable point acceptable by improving constraint
     91     *  violation, this corrector step is tried even if the regular
     92     *  primal-dual step is acceptable.
     93     */
     94    virtual bool TryCorrector(Number alpha_primal_test,
     95                              Number& alpha_primal,
     96                              SmartPtr<IteratesVector>& actual_delta);
     97
     98    /** Method for ending the current line search.  When it is called,
     99     *  the internal data should be updates, e.g., the filter might be
     100     *  augmented.  alpha_primal_test is the value of alpha that has
     101     *  been used for in the acceptence test ealier. */
     102    virtual char UpdateForNextIteration(Number alpha_primal_test);
     103
     104    /** Method for setting internal data if the watchdog procedure is
     105     *  started. */
     106    virtual void StartWatchDog();
     107
     108    /** Method for setting internal data if the watchdog procedure is
     109     *  stopped. */
     110    virtual void StopWatchDog();
    86111
    87112    /**@name Trial Point Accepting Methods. Used internally to check certain
     
    113138    //@{
    114139    /** Copy Constructor */
    115     FilterLineSearch(const FilterLineSearch&);
     140    FilterLSAcceptor(const FilterLSAcceptor&);
    116141
    117142    /** Overloaded Equals Operator */
    118     void operator=(const FilterLineSearch&);
     143    void operator=(const FilterLSAcceptor&);
    119144    //@}
    120145
     
    142167    bool ArmijoHolds(Number alpha_primal_test);
    143168
    144     /** Method to calculate alpha_min (minimum alpha before going to
    145      *  restoration
    146      */
    147     Number CalculateAlphaMin();
    148 
    149169    /** Augment the filter used on the current values of the barrier
    150170     *  objective function and the contraint violation */
    151171    void AugmentFilter();
    152 
    153     /** Method performing the backtracking line search.  The return
    154      *  value indicates if the step acceptance criteria are met.  If
    155      *  the watchdog is active, only one trial step is performed (and
    156      *  the trial values are set accordingly). */
    157     bool DoBacktrackingLineSearch(bool skip_first_trial_point,
    158                                   Number& alpha_primal,
    159                                   bool& corr_taken,
    160                                   bool& soc_taken,
    161                                   Index& n_steps,
    162                                   bool& evaluation_error,
    163                                   SmartPtr<IteratesVector>& actual_delta);
    164 
    165     /** Method for starting the watch dog.  Set all appropriate fields
    166      *  accordingly */
    167     void StartWatchDog();
    168 
    169     /** Method for stopping the watch dog.  Set all appropriate fields
    170      *  accordingly. */
    171     void StopWatchDog(SmartPtr<IteratesVector>& actual_delta);
    172 
    173     /** Method for checking if current trial point is acceptable.
    174      *  It is assumed that the delta information in ip_data is the
    175      *  search direction used in criteria.  The primal trial point has
    176      *  to be set before the call.
    177      */
    178     bool CheckAcceptabilityOfTrialPoint(Number alpha_primal);
    179172
    180173    /** Check comparison "lhs <= rhs", using machine precision based on BasVal */
     
    182175    //     allow for different relaxation parameters values
    183176    static bool Compare_le(Number lhs, Number rhs, Number BasVal);
    184 
    185     /** Method for setting the dual variables in the trial fields in
    186      *  IpData, given the search direction.  The step size for the
    187      *  bound multipliers is alpha_dual (the fraction-to-the-boundary
    188      *  step size), and the step size for the equality constraint
    189      *  multipliers depends on the choice of alpha_for_y. */
    190     void PerformDualStep(Number alpha_primal,
    191                          Number alpha_dual,
    192                          SmartPtr<IteratesVector>& delta);
    193 
    194     /** Try a step for the soft restoration phase and check if it is
    195      *  acceptable.  The step size is identical for all variables.  A
    196      *  point is accepted if it is acceptable for the original filter
    197      *  (in which case satisfies_original_filter = true on return), or
    198      *  if the primal-dual system error was decrease by at least the
    199      *  factor soft_resto_pderror_reduction_factor_.  The return value is
    200      *  true, if the trial point was acceptable. */
    201     bool TrySoftRestoStep(SmartPtr<IteratesVector>& actual_delta,
    202                           bool &satisfies_original_filter);
    203 
    204     /** Try a second order correction for the constraints.  If the
    205      *  first trial step (with incoming alpha_primal) has been reject,
    206      *  this tries up to max_soc_ second order corrections for the
    207      *  constraints.  Here, alpha_primal_test is the step size that
    208      *  has to be used in the filter acceptance tests.  On output
    209      *  actual_delta_... has been set to the steps including the
    210      *  second order correction if it has been accepted, otherwise it
    211      *  is unchanged.  If the SOC step has been accepted, alpha_primal
    212      *  has the fraction-to-the-boundary value for the SOC step on output.
    213      *  The return value is true, if an SOC step has been accepted.
    214      */
    215     bool TrySecondOrderCorrection(Number alpha_primal_test,
    216                                   Number& alpha_primal,
    217                                   SmartPtr<IteratesVector>& actual_delta);
    218 
    219     /** Try higher order corrector (for fast local convergence).  In
    220      *  contrast to a second order correction step, which tries to
    221      *  make an unacceptable point acceptable by improving constraint
    222      *  violation, this corrector step is tried even if the regular
    223      *  primal-dual step is acceptable.
    224      */
    225     bool TryCorrector(Number alpha_primal_test,
    226                       Number& alpha_primal,
    227                       SmartPtr<IteratesVector>& actual_delta);
    228 
    229     /** Perform magic steps.  Take the current values of the slacks in
    230      *  trial and replace them by better ones that lead to smaller
    231      *  values of the barrier function and less constraint
    232      *  violation. */
    233     void PerformMagicStep();
    234 
    235     /** Detect if the search direction is too small.  This should be
    236      *  true if the search direction is so small that if makes
    237      *  numerically no difference. */
    238     bool DetectTinyStep();
    239 
    240     /** Store current iterate as acceptable point */
    241     void StoreAcceptablePoint();
    242 
    243     /** Restore acceptable point into the current fields of IpData if
    244      *  found. Returns true if such as point is available. */
    245     bool RestoreAcceptablePoint();
    246 
    247     /** Method for determining if the current iterate is acceptable.
    248      *  This is a wrapper for same method from ConvergenceCheck, but
    249      *  returns false, if no ConvergenceCheck object is provided. */
    250     bool CurrentIsAcceptable();
    251177
    252178    /** @name Parameters for the filter algorithm.  Names as in the paper */
     
    266192    /** \f$ \gamma_{\alpha} \f$ */
    267193    Number alpha_min_frac_;
    268     /** factor by which search direction is to be shortened if trial
    269      *  point is rejected. */
    270     Number alpha_red_factor_;
    271194    /** Maximal number of second order correction steps */
    272195    Index max_soc_;
     
    281204    Number obj_max_inc_;
    282205
    283     /** enumeration for the different alpha_for_y_ settings */
    284     enum AlphaForYEnum {
    285       PRIMAL_ALPHA_FOR_Y=0,
    286       DUAL_ALPHA_FOR_Y,
    287       MIN_ALPHA_FOR_Y,
    288       MAX_ALPHA_FOR_Y,
    289       FULL_STEP_FOR_Y,
    290       MIN_DUAL_INFEAS_ALPHA_FOR_Y,
    291       SAFE_MIN_DUAL_INFEAS_ALPHA_FOR_Y
    292     };
    293     /** Flag indicating whether the dual step size is to be used for
    294      *  the equality constraint multipliers. If 0, the primal step
    295      *  size is used, if 1 the dual step size, and if 2, the minimum
    296      *  of both. */
    297     AlphaForYEnum alpha_for_y_;
    298 
    299     /** Flag indicating whether magic steps should be used. */
    300     bool magic_steps_;
    301206    /** enumeration for the corrector type */
    302207    enum CorrectorTypeEnum {
     
    307212    /** Type of corrector steps that should be tried. */
    308213    CorrectorTypeEnum corrector_type_;
    309     /** Flag indicating whether the line search should always accept
    310      *  the full (fraction-to-the-boundary) step. */
    311     bool accept_every_trial_step_;
    312214    /** parameter in heurstic that determines whether corrector step
    313215    should be tried. */
     
    319221     *  the monotone mu mode. */
    320222    bool skip_corr_in_monotone_mode_;
    321     /** Indicates whether problem can be expected to be infeasible.
    322      *  This will trigger requesting a tighter reduction in
    323      *  infeasibility the first time the restoration phase is
    324      *  called. */
    325     bool expect_infeasible_problem_;
    326     /** Tolerance on constraint violation for
    327      *  expect_infeasible_problem heuristic.  If the constraint
    328      *  violation becomes that than this value, the heuristic is
    329      *  disabled for the rest of the optimization run. */
    330     Number expect_infeasible_problem_ctol_;
    331     /** Reduction factor for the restoration phase that accepts steps
    332      *  reducing the optimality error ("soft restoration phase"). If
    333      *  0., then this restoration phase is not enabled. */
    334     Number soft_resto_pderror_reduction_factor_;
    335 
    336     /** Tolerance for detecting tiny steps. */
    337     Number tiny_step_tol_;
    338 
    339     /** Number of watch dog trial steps. */
    340     Index watchdog_trial_iter_max_;
    341     /** Number of shortened iterations that trigger the watchdog. */
    342     Index watchdog_shortened_iter_trigger_;
    343 
    344     /** Indicates whether the algorithm should start directly with the
    345      *  restoratin phase */
    346     bool start_with_resto_;
    347223    //@}
    348224
     
    358234     *  respect to which progress is to be made */
    359235    Number reference_gradBarrTDelta_;
    360     /** Flag indicating if the watchdog is active */
    361     bool in_watchdog_;
    362     /** Counter for shortened iterations. */
    363     Index watchdog_shortened_iter_;
    364     /** Counter for watch dog iterations */
    365     Index watchdog_trial_iter_;
    366     /** Step size for Armijo test in watch dog */
    367     Number watchdog_alpha_primal_test_;
    368236    /** Constraint violation at reference point */
    369237    Number watchdog_theta_;
     
    372240    /** Barrier gradient transpose search direction at reference point */
    373241    Number watchdog_gradBarrTDelta_;
    374     /** Watchdog reference iterate */
    375     SmartPtr<const IteratesVector> watchdog_iterate_;
    376     /** Watchdog search direction at reference point */
    377     SmartPtr<const IteratesVector> watchdog_delta_;
    378     //@}
    379 
    380     /** @name Storage for last iterate that satisfies the acceptable
    381      *  level of optimality error. */
    382     //@{
    383     SmartPtr<const IteratesVector> acceptable_iterate_;
    384242    //@}
    385243
     
    387245    Filter filter_;
    388246
    389     /** Flag indicating whether the line search is to be performed
    390      * robust (usually this is true, unless SetRigorousLineSearch is
    391      * called with false).
    392      */
    393     bool rigorous_;
    394 
    395     /** Flag indicating whether no acceptable trial point was found
    396      *  during last line search. */
    397     bool skipped_line_search_;
    398 
    399     /** Flag indicating whether we are currently in the "soft"
    400      *  restoration phase mode, in which steps are accepted if they
    401      *  reduce the optimality error (see
    402      *  soft_resto_pderror_reduction_factor) */
    403     bool in_soft_resto_phase_;
    404 
    405     /** Counter for the number of successive iterations in which the
    406      *  full step was not accepted. */
    407     Index count_successive_shortened_steps_;
    408 
    409     /** Flag indicating if a tiny step was detected in previous
    410      *  iteration */
    411     bool tiny_step_last_iteration_;
    412 
    413247    /** @name Strategy objective that are used */
    414248    //@{
    415     SmartPtr<RestorationPhase> resto_phase_;
    416249    SmartPtr<PDSystemSolver> pd_solver_;
    417     SmartPtr<ConvergenceCheck> conv_check_;
    418250    //@}
    419251  };
  • branches/dev/Algorithm/IpLineSearch.hpp

    r501 r542  
    2626    LineSearch()
    2727    {}
    28     ;
    2928
    3029    /** Default destructor */
    3130    virtual ~LineSearch()
    3231    {}
    33     ;
    3432    //@}
    3533
  • branches/dev/Algorithm/IpRestoFilterConvCheck.cpp

    r490 r542  
    1010#include "IpCompoundVector.hpp"
    1111#include "IpRestoIpoptNLP.hpp"
     12#include "IpRestoPhase.hpp"
    1213
    1314namespace Ipopt
     
    1920  RestoFilterConvergenceCheck::RestoFilterConvergenceCheck()
    2021      :
    21       orig_filter_line_search_(NULL)
     22      orig_filter_ls_acceptor_(NULL)
    2223  {
    2324    DBG_START_FUN("RestoFilterConvergenceCheck::RestoFilterConvergenceCheck()",
     
    3334
    3435  void
    35   RestoFilterConvergenceCheck::SetOrigFilterLineSearch
    36   (const FilterLineSearch& orig_filter_line_search)
     36  RestoFilterConvergenceCheck::SetOrigFilterLSAcceptor
     37  (const FilterLSAcceptor& orig_filter_ls_acceptor)
    3738  {
    38     orig_filter_line_search_ = &orig_filter_line_search;
     39    orig_filter_ls_acceptor_ = &orig_filter_ls_acceptor;
    3940  }
    4041
     
    5455      const std::string& prefix)
    5556  {
    56     DBG_ASSERT(orig_filter_line_search_ && "Need to call RestoFilterConvergenceCheck::SetOrigFilterLineSearch before Initialize");
     57    DBG_ASSERT(orig_filter_ls_acceptor_ && "Need to call RestoFilterConvergenceCheck::SetOrigFilterLineSearch before Initialize");
    5758    options.GetNumericValue("required_infeasibility_reduction", kappa_resto_, prefix);
    5859    options.GetIntegerValue("max_iter", maximum_iters_, prefix);
     
    129130                     "orig_trial_barr = %8.2e\n", orig_trial_barr);
    130131
    131       if (!orig_filter_line_search_->IsAcceptableToCurrentFilter(orig_trial_barr, orig_trial_theta)) {
     132      if (!orig_filter_ls_acceptor_->IsAcceptableToCurrentFilter(orig_trial_barr, orig_trial_theta)) {
    132133        Jnlst().Printf(J_DETAILED, J_MAIN,
    133134                       "Point is not acceptable to the original filter.\n");
    134135        status = CONTINUE;
    135136      }
    136       else if (!orig_filter_line_search_->IsAcceptableToCurrentIterate(orig_trial_barr, orig_trial_theta, true) ) {
     137      else if (!orig_filter_ls_acceptor_->IsAcceptableToCurrentIterate(orig_trial_barr, orig_trial_theta, true) ) {
    137138        Jnlst().Printf(J_DETAILED, J_MAIN,
    138139                       "Point is not acceptable to the original current point.\n");
  • branches/dev/Algorithm/IpRestoFilterConvCheck.hpp

    r501 r542  
    1111
    1212#include "IpOptErrorConvCheck.hpp"
    13 #include "IpFilterLineSearch.hpp"
     13#include "IpFilterLSAcceptor.hpp"
    1414
    1515namespace Ipopt
     
    4747     *  method must be called to finish the definition of the
    4848     *  algorithm, before Initialize is called. */
    49     void SetOrigFilterLineSearch(const FilterLineSearch& orig_filter_line_search);
     49    void SetOrigFilterLSAcceptor(const FilterLSAcceptor& orig_filter_ls_acceptor);
    5050
    5151    /** overloaded from AlgorithmStrategyObject */
     
    9393     *  that prevent the destructor of the line search object to be
    9494     *  called! */
    95     const FilterLineSearch* orig_filter_line_search_;
     95    const FilterLSAcceptor* orig_filter_ls_acceptor_;
    9696  };
    9797
  • branches/dev/Algorithm/Makefile.am

    r523 r542  
    2020        IpAugRestoSystemSolver.cpp IpAugRestoSystemSolver.hpp \
    2121        IpAugSystemSolver.hpp \
     22        IpBacktrackingLSAcceptor.hpp \
     23        IpBacktrackingLineSearch.cpp IpBacktrackingLineSearch.hpp \
    2224        IpConvCheck.hpp \
    2325        IpDefaultIterateInitializer.cpp IpDefaultIterateInitializer.hpp \
    2426        IpEqMultCalculator.hpp \
    2527        IpFilter.cpp IpFilter.hpp \
    26         IpFilterLineSearch.cpp IpFilterLineSearch.hpp \
     28        IpFilterLSAcceptor.cpp IpFilterLSAcceptor.hpp \
    2729        IpGradientScaling.cpp IpGradientScaling.hpp \
    2830        IpIpoptAlg.cpp IpIpoptAlg.hpp \
     
    7577        IpAugRestoSystemSolver.cppbak IpAugRestoSystemSolver.hppbak \
    7678        IpAugSystemSolver.hppbak \
     79        IpBacktrackingLSAcceptor.hppbak \
     80        IpBacktrackingLineSearch.cppbak IpBacktrackingLineSearch.hppbak \
    7781        IpConvCheck.hppbak \
    7882        IpDefaultIterateInitializer.cppbak IpDefaultIterateInitializer.hppbak \
    7983        IpEqMultCalculator.hppbak \
    8084        IpFilter.cppbak IpFilter.hppbak \
    81         IpFilterLineSearch.cppbak IpFilterLineSearch.hppbak \
     85        IpFilterLSAcceptor.cppbak IpFilterLSAcceptor.hppbak \
    8286        IpGradientScaling.cppbak IpGradientScaling.hppbak \
    8387        IpIpoptAlg.cppbak IpIpoptAlg.hppbak \
  • branches/dev/Algorithm/Makefile.in

    r523 r542  
    5959        IpAlgBuilder.$(OBJEXT) IpAlgorithmRegOp.$(OBJEXT) \
    6060        IpAugRestoSystemSolver.$(OBJEXT) \
     61        IpBacktrackingLineSearch.$(OBJEXT) \
    6162        IpDefaultIterateInitializer.$(OBJEXT) IpFilter.$(OBJEXT) \
    62         IpFilterLineSearch.$(OBJEXT) IpGradientScaling.$(OBJEXT) \
     63        IpFilterLSAcceptor.$(OBJEXT) IpGradientScaling.$(OBJEXT) \
    6364        IpIpoptAlg.$(OBJEXT) IpIpoptCalculatedQuantities.$(OBJEXT) \
    6465        IpIpoptData.$(OBJEXT) IpIteratesVector.$(OBJEXT) \
     
    241242        IpAugRestoSystemSolver.cpp IpAugRestoSystemSolver.hpp \
    242243        IpAugSystemSolver.hpp \
     244        IpBacktrackingLSAcceptor.hpp \
     245        IpBacktrackingLineSearch.cpp IpBacktrackingLineSearch.hpp \
    243246        IpConvCheck.hpp \
    244247        IpDefaultIterateInitializer.cpp IpDefaultIterateInitializer.hpp \
    245248        IpEqMultCalculator.hpp \
    246249        IpFilter.cpp IpFilter.hpp \
    247         IpFilterLineSearch.cpp IpFilterLineSearch.hpp \
     250        IpFilterLSAcceptor.cpp IpFilterLSAcceptor.hpp \
    248251        IpGradientScaling.cpp IpGradientScaling.hpp \
    249252        IpIpoptAlg.cpp IpIpoptAlg.hpp \
     
    290293        IpAugRestoSystemSolver.cppbak IpAugRestoSystemSolver.hppbak \
    291294        IpAugSystemSolver.hppbak \
     295        IpBacktrackingLSAcceptor.hppbak \
     296        IpBacktrackingLineSearch.cppbak IpBacktrackingLineSearch.hppbak \
    292297        IpConvCheck.hppbak \
    293298        IpDefaultIterateInitializer.cppbak IpDefaultIterateInitializer.hppbak \
    294299        IpEqMultCalculator.hppbak \
    295300        IpFilter.cppbak IpFilter.hppbak \
    296         IpFilterLineSearch.cppbak IpFilterLineSearch.hppbak \
     301        IpFilterLSAcceptor.cppbak IpFilterLSAcceptor.hppbak \
    297302        IpGradientScaling.cppbak IpGradientScaling.hppbak \
    298303        IpIpoptAlg.cppbak IpIpoptAlg.hppbak \
     
    383388@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/IpAlgorithmRegOp.Po@am__quote@
    384389@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/IpAugRestoSystemSolver.Po@am__quote@
     390@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/IpBacktrackingLineSearch.Po@am__quote@
    385391@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/IpDefaultIterateInitializer.Po@am__quote@
    386392@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/IpFilter.Po@am__quote@
    387 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/IpFilterLineSearch.Po@am__quote@
     393@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/IpFilterLSAcceptor.Po@am__quote@
    388394@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/IpGradientScaling.Po@am__quote@
    389395@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/IpIpoptAlg.Po@am__quote@
Note: See TracChangeset for help on using the changeset viewer.