Changeset 526


Ignore:
Timestamp:
Sep 27, 2005 7:03:20 PM (15 years ago)
Author:
andreasw
Message:
  • (minor) changes in the quality function search
  • added more timing statistics
Location:
branches/dev
Files:
17 edited

Legend:

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

    r520 r526  
    348348
    349349      // Compute the new barrier parameter via the oracle
    350       Number mu = free_mu_oracle_->CalculateMu();
     350      Number mu = free_mu_oracle_->CalculateMu(mu_min_, mu_max_);
    351351
    352352      mu = Max(mu, mu_min_);
     
    546546
    547547    if (IsValid(fix_mu_oracle_)) {
    548       new_mu = fix_mu_oracle_->CalculateMu();
     548      new_mu = fix_mu_oracle_->CalculateMu(mu_min_, mu_max_);
    549549    }
    550550    else {
  • branches/dev/Algorithm/IpIpoptAlg.cpp

    r523 r526  
    151151    }
    152152
     153    IpData().TimingStats().InitializeIterates.Start();
    153154    // Initialize the iterates
    154155    InitializeIterates();
     156    IpData().TimingStats().InitializeIterates.End();
    155157
    156158    if (!skip_print_problem_stats_) {
     159      IpData().TimingStats().PrintProblemStatistics.Start();
    157160      PrintProblemStatistics();
    158     }
    159 
     161      IpData().TimingStats().PrintProblemStatistics.End();
     162    }
     163
     164    IpData().TimingStats().CheckConvergence.Start();
    160165    ConvergenceCheck::ConvergenceStatus conv_status
    161166    = conv_check_->CheckConvergence();
     167    IpData().TimingStats().CheckConvergence.End();
    162168
    163169    try {
     
    165171      while (conv_status == ConvergenceCheck::CONTINUE) {
    166172        // Set the Hessian Matrix
     173        IpData().TimingStats().ActualizeHessian.Start();
    167174        ActualizeHessian();
     175        IpData().TimingStats().ActualizeHessian.End();
    168176
    169177        // do all the output for this iteration
     178        IpData().TimingStats().OutputIteration.Start();
    170179        OutputIteration();
    171180        IpData().ResetInfo();
     181        IpData().TimingStats().OutputIteration.End();
    172182
    173183        // update the barrier parameter if necessary
     184        IpData().TimingStats().UpdateBarrierParameter.Start();
    174185        UpdateBarrierParameter();
     186        IpData().TimingStats().UpdateBarrierParameter.End();
    175187
    176188        // solve the primal-dual system to get the full step
     189        IpData().TimingStats().ComputeSearchDirection.Start();
    177190        ComputeSearchDirection();
     191        IpData().TimingStats().ComputeSearchDirection.End();
    178192
    179193        // Compute the new iterate
     194        IpData().TimingStats().ComputeAcceptableTrialPoint.Start();
    180195        ComputeAcceptableTrialPoint();
    181         //ApplyFractionToBoundary();
     196        IpData().TimingStats().ComputeAcceptableTrialPoint.End();
    182197
    183198        // Accept the new iterate
     199        IpData().TimingStats().AcceptTrialPoint.Start();
    184200        AcceptTrialPoint();
     201        IpData().TimingStats().AcceptTrialPoint.End();
    185202
    186203        IpData().Set_iter_count(IpData().iter_count()+1);
    187204
     205        IpData().TimingStats().CheckConvergence.Start();
    188206        conv_status  = conv_check_->CheckConvergence();
    189       }
    190 
     207        IpData().TimingStats().CheckConvergence.End();
     208      }
     209
     210      IpData().TimingStats().OutputIteration.Start();
    191211      OutputIteration();
     212      IpData().TimingStats().OutputIteration.End();
    192213
    193214      IpData().TimingStats().OverallAlgorithm.End();
     
    205226    catch(TINY_STEP_DETECTED& exc) {
    206227      exc.ReportException(Jnlst(), J_DETAILED);
     228      IpData().TimingStats().UpdateBarrierParameter.End();
    207229      IpData().TimingStats().OverallAlgorithm.End();
    208230      return STOP_AT_TINY_STEP;
     
    210232    catch(ACCEPTABLE_POINT_REACHED& exc) {
    211233      exc.ReportException(Jnlst(), J_DETAILED);
     234      IpData().TimingStats().ComputeAcceptableTrialPoint.End();
    212235      IpData().TimingStats().OverallAlgorithm.End();
    213236      return STOP_AT_ACCEPTABLE_POINT;
     
    215238    catch(LOCALLY_INFEASIBLE& exc) {
    216239      exc.ReportException(Jnlst(), J_DETAILED);
     240      IpData().TimingStats().ComputeAcceptableTrialPoint.EndIfStarted();
     241      IpData().TimingStats().CheckConvergence.EndIfStarted();
    217242      IpData().TimingStats().OverallAlgorithm.End();
    218243      return LOCAL_INFEASIBILITY;
     
    220245    catch(RESTORATION_FAILED& exc) {
    221246      exc.ReportException(Jnlst(), J_DETAILED);
     247      IpData().TimingStats().ComputeAcceptableTrialPoint.End();
    222248      IpData().TimingStats().OverallAlgorithm.End();
    223249      return RESTORATION_FAILURE;
     
    225251    catch(FEASIBILITY_PROBLEM_SOLVED& exc) {
    226252      exc.ReportException(Jnlst(), J_DETAILED);
     253      IpData().TimingStats().ComputeAcceptableTrialPoint.End();
    227254      IpData().TimingStats().OverallAlgorithm.End();
    228255      return SUCCESS;
  • branches/dev/Algorithm/IpIpoptCalculatedQuantities.cpp

    r517 r526  
    116116      primal_frac_to_the_bound_cache_(5),
    117117      dual_frac_to_the_bound_cache_(5),
    118       slack_frac_to_the_bound_cache_(5),
    119118
    120119      curr_sigma_x_cache_(1),
     
    29462945
    29472946  Number
    2948   IpoptCalculatedQuantities::dual_frac_to_the_bound(Number tau,
    2949       const Vector& delta_z_L,
    2950       const Vector& delta_z_U,
    2951       const Vector& delta_v_L,
    2952       const Vector& delta_v_U)
     2947  IpoptCalculatedQuantities::uncached_dual_frac_to_the_bound(
     2948    Number tau,
     2949    const Vector& delta_z_L,
     2950    const Vector& delta_z_U,
     2951    const Vector& delta_v_L,
     2952    const Vector& delta_v_U)
     2953  {
     2954    DBG_START_METH("IpoptCalculatedQuantities::uncached_dual_frac_to_the_bound",
     2955                   dbg_verbosity);
     2956    Number result;
     2957
     2958    result = ip_data_->curr()->z_L()->FracToBound(delta_z_L, tau);
     2959    result = Min(result, ip_data_->curr()->z_U()->FracToBound(delta_z_U, tau));
     2960    result = Min(result, ip_data_->curr()->v_L()->FracToBound(delta_v_L, tau));
     2961    result = Min(result, ip_data_->curr()->v_U()->FracToBound(delta_v_U, tau));
     2962
     2963    return result;
     2964  }
     2965
     2966  Number
     2967  IpoptCalculatedQuantities::dual_frac_to_the_bound(
     2968    Number tau,
     2969    const Vector& delta_z_L,
     2970    const Vector& delta_z_U,
     2971    const Vector& delta_v_L,
     2972    const Vector& delta_v_U)
    29532973  {
    29542974    DBG_START_METH("IpoptCalculatedQuantities::dual_frac_to_the_bound",
     
    29973017
    29983018  Number
    2999   IpoptCalculatedQuantities::slack_frac_to_the_bound(Number tau,
    3000       const Vector& delta_x_L,
    3001       const Vector& delta_x_U,
    3002       const Vector& delta_s_L,
    3003       const Vector& delta_s_U)
     3019  IpoptCalculatedQuantities::uncached_slack_frac_to_the_bound(
     3020    Number tau,
     3021    const Vector& delta_x_L,
     3022    const Vector& delta_x_U,
     3023    const Vector& delta_s_L,
     3024    const Vector& delta_s_U)
    30043025  {
    30053026    DBG_START_METH("IpoptCalculatedQuantities::slack_frac_to_the_bound",
     
    30113032    SmartPtr<const Vector> s_L = curr_slack_s_L();
    30123033    SmartPtr<const Vector> s_U = curr_slack_s_U();
    3013     std::vector<const TaggedObject*> tdeps(8);
    3014     tdeps[0] = GetRawPtr(x_L);
    3015     tdeps[1] = GetRawPtr(x_U);
    3016     tdeps[2] = GetRawPtr(s_L);
    3017     tdeps[3] = GetRawPtr(s_U);
    3018     tdeps[4] = &delta_x_L;
    3019     tdeps[5] = &delta_x_U;
    3020     tdeps[6] = &delta_s_L;
    3021     tdeps[7] = &delta_s_U;
    3022 
    3023     std::vector<Number> sdeps(1);
    3024     sdeps[0] = tau;
    3025 
    3026     if (!slack_frac_to_the_bound_cache_.GetCachedResult(result, tdeps, sdeps)) {
    3027       result = x_L->FracToBound(delta_x_L, tau);
    3028       result = Min(result, x_U->FracToBound(delta_x_U, tau));
    3029       result = Min(result, s_L->FracToBound(delta_s_L, tau));
    3030       result = Min(result, s_U->FracToBound(delta_s_U, tau));
    3031 
    3032       slack_frac_to_the_bound_cache_.AddCachedResult(result, tdeps, sdeps);
    3033     }
     3034
     3035    result = x_L->FracToBound(delta_x_L, tau);
     3036    result = Min(result, x_U->FracToBound(delta_x_U, tau));
     3037    result = Min(result, s_L->FracToBound(delta_s_L, tau));
     3038    result = Min(result, s_U->FracToBound(delta_s_U, tau));
    30343039
    30353040    return result;
  • branches/dev/Algorithm/IpIpoptCalculatedQuantities.hpp

    r517 r526  
    315315                                  const Vector& delta_v_U);
    316316    /** Fraction to the boundary from (current) dual variables z and v
     317     *  for a given step, without caching */
     318    Number uncached_dual_frac_to_the_bound(Number tau,
     319                                           const Vector& delta_z_L,
     320                                           const Vector& delta_z_U,
     321                                           const Vector& delta_v_L,
     322                                           const Vector& delta_v_U);
     323    /** Fraction to the boundary from (current) dual variables z and v
    317324     *  for internal (current) step */
    318325    Number curr_dual_frac_to_the_bound(Number tau);
     
    322329     *  to the boundary step size, but if it is cheaper to provide the
    323330     *  steps in the slacks directly (e.g. when the primal step sizes
    324      *  are only temporary), the this method is more efficient. */
    325     Number slack_frac_to_the_bound(Number tau,
    326                                    const Vector& delta_x_L,
    327                                    const Vector& delta_x_U,
    328                                    const Vector& delta_s_L,
    329                                    const Vector& delta_s_U);
     331     *  are only temporary), the this method is more efficient.  This
     332     *  method does not cache computations. */
     333    Number uncached_slack_frac_to_the_bound(Number tau,
     334                                            const Vector& delta_x_L,
     335                                            const Vector& delta_x_U,
     336                                            const Vector& delta_s_L,
     337                                            const Vector& delta_s_U);
    330338    //@}
    331339
     
    516524    CachedResults<Number> primal_frac_to_the_bound_cache_;
    517525    CachedResults<Number> dual_frac_to_the_bound_cache_;
    518     CachedResults<Number> slack_frac_to_the_bound_cache_;
    519526    //@}
    520527
  • branches/dev/Algorithm/IpLoqoMuOracle.cpp

    r465 r526  
    4242  }
    4343
    44   Number LoqoMuOracle::CalculateMu()
     44  Number LoqoMuOracle::CalculateMu(Number mu_min, Number mu_max)
    4545  {
    4646    DBG_START_METH("LoqoMuOracle::CalculateMu",
     
    7070    IpData().Append_info_string(ssigma);
    7171
    72     return mu;
     72    return Max(Min(mu_max, mu), mu_min);
    7373  }
    7474
  • branches/dev/Algorithm/IpLoqoMuOracle.hpp

    r501 r526  
    3636     *  could be used in the current iteration (using the LOQO formula).
    3737     */
    38     virtual Number CalculateMu();
     38    virtual Number CalculateMu(Number mu_min, Number mu_max);
    3939
    4040  private:
  • branches/dev/Algorithm/IpMuOracle.hpp

    r501 r526  
    3939
    4040    /** Method for computing the value of the barrier parameter that
    41      *  could be used in the current iteration.
     41     *  could be used in the current iteration.  Here, mu_min and
     42     *  mu_max are the lower and upper bounds on acceptable values for
     43     *  the barrier parameter
    4244     */
    43     virtual Number CalculateMu() = 0;
     45    virtual Number CalculateMu(Number mu_min, Number mu_max) = 0;
    4446
    4547  private:
  • branches/dev/Algorithm/IpOrigIpoptNLP.cpp

    r517 r526  
    308308      f_evals_++;
    309309      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
     310      f_eval_time_.Start();
    310311      bool success = nlp_->Eval_f(*unscaled_x, ret);
     312      f_eval_time_.End();
    311313      DBG_PRINT((1, "success = %d ret = %e\n", success, ret));
    312314      ASSERT_EXCEPTION(success && IsFiniteNumber(ret), Eval_Error,
     
    334336
    335337      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
     338      grad_f_eval_time_.Start();
    336339      bool success = nlp_->Eval_grad_f(*unscaled_x, *unscaled_grad_f);
     340      grad_f_eval_time_.End();
    337341      ASSERT_EXCEPTION(success && IsFiniteNumber(unscaled_grad_f->Nrm2()),
    338342                       Eval_Error, "Error evaluating the gradient of the objective function");
     
    369373        c_evals_++;
    370374        SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
     375        c_eval_time_.Start();
    371376        bool success = nlp_->Eval_c(*unscaled_x, *unscaled_c);
     377        c_eval_time_.End();
    372378        ASSERT_EXCEPTION(success && IsFiniteNumber(unscaled_c->Nrm2()),
    373379                         Eval_Error, "Error evaluating the equality constraints");
     
    401407        DBG_PRINT_VECTOR(2, "scaled_x", x);
    402408        SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
     409        d_eval_time_.Start();
    403410        bool success = nlp_->Eval_d(*unscaled_x, *unscaled_d);
     411        d_eval_time_.End();
    404412        DBG_PRINT_VECTOR(2, "unscaled_d", *unscaled_d);
    405413        ASSERT_EXCEPTION(success && IsFiniteNumber(unscaled_d->Nrm2()),
     
    432440
    433441        SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
     442        jac_c_eval_time_.Start();
    434443        bool success = nlp_->Eval_jac_c(*unscaled_x, *unscaled_jac_c);
     444        jac_c_eval_time_.End();
    435445        ASSERT_EXCEPTION(success, Eval_Error, "Error evaluating the jacobian of the equality constraints");
    436446        retValue = NLP_scaling()->apply_jac_c_scaling(ConstPtr(unscaled_jac_c));
     
    461471
    462472        SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
     473        jac_d_eval_time_.Start();
    463474        bool success = nlp_->Eval_jac_d(*unscaled_x, *unscaled_jac_d);
     475        jac_d_eval_time_.End();
    464476        ASSERT_EXCEPTION(success, Eval_Error, "Error evaluating the jacobian of the inequality constraints");
    465477        retValue = NLP_scaling()->apply_jac_d_scaling(ConstPtr(unscaled_jac_d));
     
    493505      SmartPtr<const Vector> unscaled_yd = NLP_scaling()->apply_vector_scaling_d(&yd);
    494506      Number scaled_obj_factor = NLP_scaling()->apply_obj_scaling(obj_factor);
     507      h_eval_time_.Start();
    495508      bool success = nlp_->Eval_h(*unscaled_x, scaled_obj_factor, *unscaled_yc, *unscaled_yd, *unscaled_h);
     509      h_eval_time_.End();
    496510      ASSERT_EXCEPTION(success, Eval_Error, "Error evaluating the hessian of the lagrangian");
    497511      retValue = NLP_scaling()->apply_hessian_scaling(ConstPtr(unscaled_h));
     
    625639  }
    626640
     641  void
     642  OrigIpoptNLP::PrintTimingStatistics(
     643    Journalist& jnlst,
     644    EJournalLevel level,
     645    EJournalCategory category) const
     646  {
     647    if (!jnlst.ProduceOutput(level, category))
     648      return;
     649
     650    jnlst.Printf(level, category,
     651                 "Function Evaluations................: %10.2f\n",
     652                 f_eval_time_.TotalTime()+
     653                 c_eval_time_.TotalTime()+
     654                 d_eval_time_.TotalTime()+
     655                 jac_c_eval_time_.TotalTime()+
     656                 jac_d_eval_time_.TotalTime()+
     657                 h_eval_time_.TotalTime());
     658    jnlst.Printf(level, category,
     659                 " Objective function.................: %10.2f\n",
     660                 f_eval_time_.TotalTime());
     661    jnlst.Printf(level, category,
     662                 " Equality constraints...............: %10.2f\n",
     663                 c_eval_time_.TotalTime());
     664    jnlst.Printf(level, category,
     665                 " Inequality constraints.............: %10.2f\n",
     666                 d_eval_time_.TotalTime());
     667    jnlst.Printf(level, category,
     668                 " Equality constraint Jacobian.......: %10.2f\n",
     669                 jac_c_eval_time_.TotalTime());
     670    jnlst.Printf(level, category,
     671                 " Inequality constraint Jacobian.....: %10.2f\n",
     672                 jac_d_eval_time_.TotalTime());
     673    jnlst.Printf(level, category,
     674                 " Lagrangian Hessian.................: %10.2f\n",
     675                 h_eval_time_.TotalTime());
     676  }
     677
    627678} // namespace Ipopt
  • branches/dev/Algorithm/IpOrigIpoptNLP.hpp

    r517 r526  
    1212#include "IpIpoptNLP.hpp"
    1313#include "IpException.hpp"
     14#include "IpTimingStatistics.hpp"
    1415
    1516namespace Ipopt
     
    225226    }
    226227
     228    void PrintTimingStatistics(Journalist& jnlst,
     229                               EJournalLevel level,
     230                               EJournalCategory category) const;
     231
    227232  private:
    228233    /** journalist */
     
    350355    Index jac_d_evals_;
    351356    Index h_evals_;
     357    //@}
    352358
    353359    /** Flag indicating if initialization method has been called */
    354360    bool initialized_;
     361
     362    /**@name Timing statistics for the function evaluations. */
     363    //@{
     364    TimedTask f_eval_time_;
     365    TimedTask grad_f_eval_time_;
     366    TimedTask c_eval_time_;
     367    TimedTask jac_c_eval_time_;
     368    TimedTask d_eval_time_;
     369    TimedTask jac_d_eval_time_;
     370    TimedTask h_eval_time_;
    355371    //@}
    356372  };
  • branches/dev/Algorithm/IpProbingMuOracle.cpp

    r465 r526  
    4949  }
    5050
    51   Number ProbingMuOracle::CalculateMu()
     51  Number ProbingMuOracle::CalculateMu(Number mu_min, Number mu_max)
    5252  {
    5353    DBG_START_METH("ProbingMuOracle::CalculateMu",
     
    107107
    108108    // now compute the average complementarity at the affine step
     109    // ToDo shoot for mu_min instead of 0?
    109110    Number mu_aff = CalculateAffineMu(alpha_primal_aff, alpha_dual_aff, *step);
    110111    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
     
    138139    IpData().Append_info_string(ssigma);
    139140
    140     return mu;
     141    return Max(Min(mu, mu_max), mu_min);
    141142  }
    142143
  • branches/dev/Algorithm/IpProbingMuOracle.hpp

    r501 r526  
    3737     *  could be used in the current iteration (using the LOQO formula).
    3838     */
    39     virtual Number CalculateMu();
     39    virtual Number CalculateMu(Number mu_min, Number mu_max);
    4040
    4141    /** Methods for IpoptType */
  • branches/dev/Algorithm/IpQualityFunctionMuOracle.cpp

    r524 r526  
    4646      tmp_z_U_(NULL),
    4747      tmp_v_L_(NULL),
    48       tmp_v_U_(NULL)
     48      tmp_v_U_(NULL),
     49
     50      count_qf_evals_(0)
    4951  {
    5052    DBG_ASSERT(IsValid(pd_solver_));
     
    8991      "than dual and primal infeasibilities.");
    9092    roptions->AddLowerBoundedIntegerOption(
    91       "max_bisection_steps",
     93      "quality_function_max_section_steps",
    9294      "Maximum number of search steps during direct search procedure "
    9395      "determining the optimal centering parameter.",
    9496      0, 8,
    95       "The bisection search is performed for the quality function based mu "
    96       "oractle.");
     97      "The golden section search is performed for the quality function based "
     98      "mu oractle.");
    9799    roptions->AddBoundedNumberOption(
    98       "bisection_tol",
    99       "Tolerance for the bisection search procedure determining "
    100       "the optimal centering parameter.",
    101       0.0, true, 1.0, true,
    102       1e-3,
    103       "The bisection search is performed for the quality function based mu "
     100      "quality_function_section_sigma_tol",
     101      "Tolerance for the section search procedure determining "
     102      "the optimal centering parameter (in sigma space).",
     103      0.0, false, 1.0, true,
     104      1e-2,
     105      "The golden section search is performed for the quality function based "
     106      "mu oractle.");
     107    roptions->AddBoundedNumberOption(
     108      "quality_function_section_qf_tol",
     109      "Tolerance for the golden section search procedure determining "
     110      "the optimal centering parameter (in the function value space).",
     111      0.0, false, 1.0, true,
     112      0e-2,
     113      "The section search is performed for the quality function based mu "
    104114      "oractle.");
    105115  }
     
    119129    options.GetEnumValue("quality_function_balancing_term", enum_int, prefix);
    120130    quality_function_balancing_term_ = BalancingTermEnum(enum_int);
    121     options.GetIntegerValue("max_bisection_steps", max_bisection_steps_, prefix);
    122     options.GetNumericValue("bisection_tol", bisection_tol_, prefix);
     131    options.GetIntegerValue("quality_function_max_section_steps",
     132                            quality_function_max_section_steps_, prefix);
     133    options.GetNumericValue("quality_function_section_sigma_tol",
     134                            quality_function_section_sigma_tol_, prefix);
     135    options.GetNumericValue("quality_function_section_qf_tol",
     136                            quality_function_section_qf_tol_, prefix);
    123137
    124138    return true;
    125139  }
    126140
    127   Number QualityFunctionMuOracle::CalculateMu()
     141  Number QualityFunctionMuOracle::CalculateMu(Number mu_min, Number mu_max)
    128142  {
    129143    DBG_START_METH("QualityFunctionMuOracle::CalculateMu",
    130144                   dbg_verbosity);
     145
     146    count_qf_evals_ = 0;
    131147
    132148    ///////////////////////////////////////////////////////////////////////////
     
    247263
    248264    Number sigma;
    249     if (max_bisection_steps_>0) {
    250       // Now we do an search for the best centering parameter, that
    251       // gives us the lower value of a quality function, using golden
    252       // bisection
    253 
    254       Number sigma_up = Min(1., sigma_max_);
    255       Number sigma_lo = 1e-9/avrg_compl;
    256       sigma = PerformGoldenBisection(sigma_up, sigma_lo, bisection_tol_,
    257                                      //sigma = PerformGoldenBisectionLog(sigma_up, sigma_lo, bisection_tol_,
    258                                      *step_aff_x_L,
    259                                      *step_aff_x_U,
    260                                      *step_aff_s_L,
    261                                      *step_aff_s_U,
    262                                      *step_aff->y_c(),
    263                                      *step_aff->y_d(),
    264                                      *step_aff->z_L(),
    265                                      *step_aff->z_U(),
    266                                      *step_aff->v_L(),
    267                                      *step_aff->v_U(),
    268                                      *step_cen_x_L,
    269                                      *step_cen_x_U,
    270                                      *step_cen_s_L,
    271                                      *step_cen_s_U,
    272                                      *step_cen->y_c(),
    273                                      *step_cen->y_d(),
    274                                      *step_cen->z_L(),
    275                                      *step_cen->z_U(),
    276                                      *step_cen->v_L(),
    277                                      *step_cen->v_U());
    278 
    279       if (sigma_max_ > 1. && sigma >= 1.-2*bisection_tol_) {
    280         // It seems that the optimal value might be larger than one.
    281         sigma_up = sigma_max_;
    282         sigma_lo = sigma;
    283         sigma = PerformGoldenBisection(sigma_up, sigma_lo, bisection_tol_,
    284                                        //sigma = PerformGoldenBisectionLog(sigma_up, sigma_lo, bisection_tol_,
    285                                        *step_aff_x_L,
    286                                        *step_aff_x_U,
    287                                        *step_aff_s_L,
    288                                        *step_aff_s_U,
    289                                        *step_aff->y_c(),
    290                                        *step_aff->y_d(),
    291                                        *step_aff->z_L(),
    292                                        *step_aff->z_U(),
    293                                        *step_aff->v_L(),
    294                                        *step_aff->v_U(),
    295                                        *step_cen_x_L,
    296                                        *step_cen_x_U,
    297                                        *step_cen_s_L,
    298                                        *step_cen_s_U,
    299                                        *step_cen->y_c(),
    300                                        *step_cen->y_d(),
    301                                        *step_cen->z_L(),
    302                                        *step_cen->z_U(),
    303                                        *step_cen->v_L(),
    304                                        *step_cen->v_U());
    305       }
    306 
    307       //#define tracequalityfunction
     265
     266    // First we determine whether we want to search for a value of
     267    // sigma larger or smaller than 1.  For this, we estimate the
     268    // slope of the quality function at sigma=1.
     269    Number qf_1 = CalculateQualityFunction(1.,
     270                                           *step_aff_x_L,
     271                                           *step_aff_x_U,
     272                                           *step_aff_s_L,
     273                                           *step_aff_s_U,
     274                                           *step_aff->y_c(),
     275                                           *step_aff->y_d(),
     276                                           *step_aff->z_L(),
     277                                           *step_aff->z_U(),
     278                                           *step_aff->v_L(),
     279                                           *step_aff->v_U(),
     280                                           *step_cen_x_L,
     281                                           *step_cen_x_U,
     282                                           *step_cen_s_L,
     283                                           *step_cen_s_U,
     284                                           *step_cen->y_c(),
     285                                           *step_cen->y_d(),
     286                                           *step_cen->z_L(),
     287                                           *step_cen->z_U(),
     288                                           *step_cen->v_L(),
     289                                           *step_cen->v_U());
     290
     291    Number sigma_1minus = 1.-quality_function_section_sigma_tol_;
     292    Number qf_1minus = CalculateQualityFunction(sigma_1minus,
     293                       *step_aff_x_L,
     294                       *step_aff_x_U,
     295                       *step_aff_s_L,
     296                       *step_aff_s_U,
     297                       *step_aff->y_c(),
     298                       *step_aff->y_d(),
     299                       *step_aff->z_L(),
     300                       *step_aff->z_U(),
     301                       *step_aff->v_L(),
     302                       *step_aff->v_U(),
     303                       *step_cen_x_L,
     304                       *step_cen_x_U,
     305                       *step_cen_s_L,
     306                       *step_cen_s_U,
     307                       *step_cen->y_c(),
     308                       *step_cen->y_d(),
     309                       *step_cen->z_L(),
     310                       *step_cen->z_U(),
     311                       *step_cen->v_L(),
     312                       *step_cen->v_U());
     313
     314    if (qf_1minus > qf_1) {
     315      // It seems that the quality function decreases for values
     316      // larger than sigma, so perform golden section search for sigma
     317      // > 1.
     318      Number sigma_up = Min(sigma_max_, mu_max/avrg_compl);
     319      Number sigma_lo = 1.;
     320      // ToDo maybe we should use different tolerances for sigma>1
     321      sigma = PerformGoldenSection(sigma_up, -100., sigma_lo, qf_1,
     322                                   quality_function_section_sigma_tol_,
     323                                   quality_function_section_qf_tol_,
     324                                   *step_aff_x_L,
     325                                   *step_aff_x_U,
     326                                   *step_aff_s_L,
     327                                   *step_aff_s_U,
     328                                   *step_aff->y_c(),
     329                                   *step_aff->y_d(),
     330                                   *step_aff->z_L(),
     331                                   *step_aff->z_U(),
     332                                   *step_aff->v_L(),
     333                                   *step_aff->v_U(),
     334                                   *step_cen_x_L,
     335                                   *step_cen_x_U,
     336                                   *step_cen_s_L,
     337                                   *step_cen_s_U,
     338                                   *step_cen->y_c(),
     339                                   *step_cen->y_d(),
     340                                   *step_cen->z_L(),
     341                                   *step_cen->z_U(),
     342                                   *step_cen->v_L(),
     343                                   *step_cen->v_U());
     344    }
     345    else {
     346      // Search for sigma less than 1
     347
     348      Number sigma_lo = mu_min/avrg_compl;
     349      Number sigma_up = Max(sigma_lo, sigma_1minus);
     350      // Todo: What if lower bound is less than sigma_lo? Skip search?
     351      sigma = PerformGoldenSection(sigma_up, qf_1minus, sigma_lo, -100.,
     352                                   quality_function_section_sigma_tol_,
     353                                   quality_function_section_qf_tol_,
     354                                   *step_aff_x_L,
     355                                   *step_aff_x_U,
     356                                   *step_aff_s_L,
     357                                   *step_aff_s_U,
     358                                   *step_aff->y_c(),
     359                                   *step_aff->y_d(),
     360                                   *step_aff->z_L(),
     361                                   *step_aff->z_U(),
     362                                   *step_aff->v_L(),
     363                                   *step_aff->v_U(),
     364                                   *step_cen_x_L,
     365                                   *step_cen_x_U,
     366                                   *step_cen_s_L,
     367                                   *step_cen_s_U,
     368                                   *step_cen->y_c(),
     369                                   *step_cen->y_d(),
     370                                   *step_cen->z_L(),
     371                                   *step_cen->z_U(),
     372                                   *step_cen->v_L(),
     373                                   *step_cen->v_U());
     374    }
     375
     376    //#define tracequalityfunction
    308377#ifdef tracequalityfunction
    309       //DELETEME
    310       Jnlst().Printf(J_MOREDETAILED, J_BARRIER_UPDATE,
    311                      "\n");
    312       Number base = 1.2;
    313       for (Index l=30;l>=(Index)trunc(-(log(avrg_compl)-log(1e-9))/log(base));l--) {
    314         Number sig = pow(base, l);
    315         CalculateQualityFunction(sig,
    316                                  *step_aff_x_L,
    317                                  *step_aff_x_U,
    318                                  *step_aff_s_L,
    319                                  *step_aff_s_U,
    320                                  *step_aff->y_c(),
    321                                  *step_aff->y_d(),
    322                                  *step_aff->z_L(),
    323                                  *step_aff->z_U(),
    324                                  *step_aff->v_L(),
    325                                  *step_aff->v_U(),
    326                                  *step_cen_x_L,
    327                                  *step_cen_x_U,
    328                                  *step_cen_s_L,
    329                                  *step_cen_s_U,
    330                                  *step_cen->y_c(),
    331                                  *step_cen->y_d(),
    332                                  *step_cen->z_L(),
    333                                  *step_cen->z_U(),
    334                                  *step_cen->v_L(),
    335                                  *step_cen->v_U());
    336       }
     378    char fname[100];
     379    sprintf(fname, "qf_values_%d.dat", IpData().iter_count());
     380    FILE* fid = fopen(fname, "w");
     381
     382    Number sigma_1 = sigma_max_;
     383    Number sigma_2 = 1e-9/avrg_compl;
     384    Number sigma_trace = sigma_1;
     385    while(sigma_trace > sigma_2) {
     386      Number qf = CalculateQualityFunction(sigma_trace,
     387                                           *step_aff_x_L,
     388                                           *step_aff_x_U,
     389                                           *step_aff_s_L,
     390                                           *step_aff_s_U,
     391                                           *step_aff->y_c(),
     392                                           *step_aff->y_d(),
     393                                           *step_aff->z_L(),
     394                                           *step_aff->z_U(),
     395                                           *step_aff->v_L(),
     396                                           *step_aff->v_U(),
     397                                           *step_cen_x_L,
     398                                           *step_cen_x_U,
     399                                           *step_cen_s_L,
     400                                           *step_cen_s_U,
     401                                           *step_cen->y_c(),
     402                                           *step_cen->y_d(),
     403                                           *step_cen->z_L(),
     404                                           *step_cen->z_U(),
     405                                           *step_cen->v_L(),
     406                                           *step_cen->v_U());
     407      fprintf(fid, "%9.2e %25.16e\n", sigma_trace, qf);
     408      sigma_trace /= 1.1;
     409    }
     410    fclose(fid);
    337411#endif
    338 
    339     }
    340     else {
    341       Index l;
    342       Index l_best;
    343       Number q_best;
    344 
    345       Number base = 1.2;
    346       l = 20;
    347       l_best = l;
    348       sigma = pow(base, l);
    349       q_best = CalculateQualityFunction(sigma,
    350                                         *step_aff_x_L,
    351                                         *step_aff_x_U,
    352                                         *step_aff_s_L,
    353                                         *step_aff_s_U,
    354                                         *step_aff->y_c(),
    355                                         *step_aff->y_d(),
    356                                         *step_aff->z_L(),
    357                                         *step_aff->z_U(),
    358                                         *step_aff->v_L(),
    359                                         *step_aff->v_U(),
    360                                         *step_cen_x_L,
    361                                         *step_cen_x_U,
    362                                         *step_cen_s_L,
    363                                         *step_cen_s_U,
    364                                         *step_cen->y_c(),
    365                                         *step_cen->y_d(),
    366                                         *step_cen->z_L(),
    367                                         *step_cen->z_U(),
    368                                         *step_cen->v_L(),
    369                                         *step_cen->v_U());
    370       Index l_min = (Index)(-(log(avrg_compl)-log(1e-9))/log(base))-1;
    371       for (; l>=l_min; l--) {
    372         sigma = pow(base, l);
    373         Number q = CalculateQualityFunction(sigma,
    374                                             *step_aff_x_L,
    375                                             *step_aff_x_U,
    376                                             *step_aff_s_L,
    377                                             *step_aff_s_U,
    378                                             *step_aff->y_c(),
    379                                             *step_aff->y_d(),
    380                                             *step_aff->z_L(),
    381                                             *step_aff->z_U(),
    382                                             *step_aff->v_L(),
    383                                             *step_aff->v_U(),
    384                                             *step_cen_x_L,
    385                                             *step_cen_x_U,
    386                                             *step_cen_s_L,
    387                                             *step_cen_s_U,
    388                                             *step_cen->y_c(),
    389                                             *step_cen->y_d(),
    390                                             *step_cen->z_L(),
    391                                             *step_cen->z_U(),
    392                                             *step_cen->v_L(),
    393                                             *step_cen->v_U());
    394         if (q<=q_best) {
    395           q_best = q;
    396           l_best = l;
    397         }
    398       }
    399 
    400       sigma = pow(base, l_best);
    401     }
    402412
    403413    // End timing of quality function search
     
    446456    sprintf(ssigma, " sigma=%8.2e", sigma);
    447457    IpData().Append_info_string(ssigma);
     458    sprintf(ssigma, " qf=%d", count_qf_evals_);
     459    IpData().Append_info_string(ssigma);
     460    /*
    448461    sprintf(ssigma, " xi=%8.2e ", IpCq().curr_centrality_measure());
    449462    IpData().Append_info_string(ssigma);
     
    451464      IpData().Append_info_string("LARGESIGMA");
    452465    }
     466    */
    453467
    454468    return mu;
     
    481495    DBG_START_METH("QualityFunctionMuOracle::CalculateQualityFunction",
    482496                   dbg_verbosity);
     497    // DELETEME
     498    count_qf_evals_++;
    483499
    484500    Index n_dual = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim();
     
    501517    IpData().TimingStats().Task2.Start();
    502518    Number tau = IpData().curr_tau();
    503     Number alpha_primal = IpCq().slack_frac_to_the_bound(tau,
     519    Number alpha_primal = IpCq().uncached_slack_frac_to_the_bound(tau,
    504520                          *tmp_step_x_L_,
    505521                          *tmp_step_x_U_,
     
    507523                          *tmp_step_s_U_);
    508524
    509     Number alpha_dual = IpCq().dual_frac_to_the_bound(tau,
     525    Number alpha_dual = IpCq().uncached_dual_frac_to_the_bound(tau,
    510526                        *tmp_step_z_L_,
    511527                        *tmp_step_z_U_,
     
    514530    IpData().TimingStats().Task2.End();
    515531
    516     Number xi; // centrality measure
     532    Number xi = 0.; // centrality measure
    517533
    518534    IpData().TimingStats().Task1.Start();
     
    547563    DBG_PRINT_VECTOR(2, "compl_s_L", *tmp_slack_s_L_);
    548564    DBG_PRINT_VECTOR(2, "compl_s_U", *tmp_slack_s_U_);
    549 
    550     IpData().TimingStats().Task4.Start();
    551     xi = IpCq().CalcCentralityMeasure(*tmp_slack_x_L_, *tmp_slack_x_U_,
    552                                       *tmp_slack_s_L_, *tmp_slack_s_U_);
    553     IpData().TimingStats().Task4.End();
    554565
    555566    Number dual_inf=-1.;
     
    630641    Number quality_function = dual_inf + primal_inf + compl_inf;
    631642
     643    if (quality_function_centrality_!=CEN_NONE) {
     644      IpData().TimingStats().Task4.Start();
     645      xi = IpCq().CalcCentralityMeasure(*tmp_slack_x_L_, *tmp_slack_x_U_,
     646                                        *tmp_slack_s_L_, *tmp_slack_s_U_);
     647      IpData().TimingStats().Task4.End();
     648    }
    632649    switch (quality_function_centrality_) {
    633650      case CEN_NONE:
     
    666683
    667684  Number
    668   QualityFunctionMuOracle::PerformGoldenBisection
    669   (Number sigma_up,
    670    Number sigma_lo,
    671    Number tol,
     685  QualityFunctionMuOracle::PerformGoldenSection
     686  (Number sigma_up_in,
     687   Number q_up,
     688   Number sigma_lo_in,
     689   Number q_lo,
     690   Number sigma_tol,
     691   Number qf_tol,
    672692   const Vector& step_aff_x_L,
    673693   const Vector& step_aff_x_U,
     
    692712  )
    693713  {
     714    Number sigma_up = ScaleSigma(sigma_up_in);
     715    Number sigma_lo = ScaleSigma(sigma_lo_in);
     716
    694717    Number sigma;
    695     Number sigma_up_in = sigma_up;
    696     Number sigma_lo_in = sigma_lo;
    697718    Number gfac = (3.-sqrt(5.))/2.;
    698719    Number sigma_mid1 = sigma_lo + gfac*(sigma_up-sigma_lo);
    699720    Number sigma_mid2 = sigma_lo + (1.-gfac)*(sigma_up-sigma_lo);
    700721
    701     Number qmid1 = CalculateQualityFunction(sigma_mid1,
     722    Number qmid1 = CalculateQualityFunction(UnscaleSigma(sigma_mid1),
    702723                                            step_aff_x_L,
    703724                                            step_aff_x_U,
     
    720741                                            step_cen_v_L,
    721742                                            step_cen_v_U);
    722     Number qmid2 = CalculateQualityFunction(sigma_mid2,
     743    Number qmid2 = CalculateQualityFunction(UnscaleSigma(sigma_mid2),
    723744                                            step_aff_x_L,
    724745                                            step_aff_x_U,
     
    742763                                            step_cen_v_U);
    743764
    744     Index nbisections = 0;
    745     while ((sigma_up-sigma_lo)>=tol*sigma_up && nbisections<max_bisection_steps_) {
    746       nbisections++;
     765    Index nsections = 0;
     766    while ((sigma_up-sigma_lo)>=sigma_tol*sigma_up &&
     767           //while ((sigma_up-sigma_lo)>=sigma_tol &&  // Note we are using the non-relative criterion here for sigma
     768           (1.-Min(q_lo, q_up, qmid1, qmid2)/Max(q_lo, q_up, qmid1, qmid2))>=qf_tol &&
     769           nsections<quality_function_max_section_steps_) {
     770      nsections++;
     771      //printf("sigma_lo=%e sigma_mid1=%e sigma_mid2=%e sigma_up=%e\n",sigma_lo, sigma_mid1, sigma_mid2, sigma_up);
    747772      if (qmid1 > qmid2) {
    748773        sigma_lo = sigma_mid1;
     774        q_lo = qmid1;
    749775        sigma_mid1 = sigma_mid2;
    750776        qmid1 = qmid2;
    751777        sigma_mid2 = sigma_lo + (1.-gfac)*(sigma_up-sigma_lo);
    752         qmid2 = CalculateQualityFunction(sigma_mid2,
     778        qmid2 = CalculateQualityFunction(UnscaleSigma(sigma_mid2),
    753779                                         step_aff_x_L,
    754780                                         step_aff_x_U,
     
    774800      else {
    775801        sigma_up = sigma_mid2;
     802        q_up = qmid2;
    776803        sigma_mid2 = sigma_mid1;
    777804        qmid2 = qmid1;
    778805        sigma_mid1 = sigma_lo + gfac*(sigma_up-sigma_lo);
    779         qmid1 = CalculateQualityFunction(sigma_mid1,
     806        qmid1 = CalculateQualityFunction(UnscaleSigma(sigma_mid1),
    780807                                         step_aff_x_L,
    781808                                         step_aff_x_U,
     
    801828    }
    802829
    803     Number q;
    804     if (qmid1 < qmid2) {
    805       sigma = sigma_mid1;
    806       q = qmid1;
     830    if ((sigma_up-sigma_lo)>=sigma_tol*sigma_up &&
     831        (1.-Min(q_lo, q_up, qmid1, qmid2)/Max(q_lo, q_up, qmid1, qmid2))<qf_tol) {
     832      // The qf tolerance make it stop
     833      IpData().Append_info_string("qf_tol ");
     834      Number qf_min = Min(q_lo, q_up, qmid1, qmid2);
     835      DBG_ASSERT(qf_min>-100.);
     836      if (qf_min == q_lo) {
     837        sigma = sigma_lo;
     838      }
     839      else if (qf_min == qmid1) {
     840        sigma = sigma_mid1;
     841      }
     842      else if (qf_min == qmid2) {
     843        sigma = sigma_mid2;
     844      }
     845      else {
     846        sigma = sigma_up;
     847      }
    807848    }
    808849    else {
    809       sigma = sigma_mid2;
    810       q = qmid2;
    811     }
    812     if (sigma_up == sigma_up_in) {
    813       Number qtmp = CalculateQualityFunction(sigma_up,
    814                                              step_aff_x_L,
    815                                              step_aff_x_U,
    816                                              step_aff_s_L,
    817                                              step_aff_s_U,
    818                                              step_aff_y_c,
    819                                              step_aff_y_d,
    820                                              step_aff_z_L,
    821                                              step_aff_z_U,
    822                                              step_aff_v_L,
    823                                              step_aff_v_U,
    824                                              step_cen_x_L,
    825                                              step_cen_x_U,
    826                                              step_cen_s_L,
    827                                              step_cen_s_U,
    828                                              step_cen_y_c,
    829                                              step_cen_y_d,
    830                                              step_cen_z_L,
    831                                              step_cen_z_U,
    832                                              step_cen_v_L,
    833                                              step_cen_v_U);
    834       if (qtmp < q) {
    835         sigma = sigma_up;
    836         q = qtmp;
    837       }
    838     }
    839     else if (sigma_lo == sigma_lo_in) {
    840       Number qtmp = CalculateQualityFunction(sigma_lo,
    841                                              step_aff_x_L,
    842                                              step_aff_x_U,
    843                                              step_aff_s_L,
    844                                              step_aff_s_U,
    845                                              step_aff_y_c,
    846                                              step_aff_y_d,
    847                                              step_aff_z_L,
    848                                              step_aff_z_U,
    849                                              step_aff_v_L,
    850                                              step_aff_v_U,
    851                                              step_cen_x_L,
    852                                              step_cen_x_U,
    853                                              step_cen_s_L,
    854                                              step_cen_s_U,
    855                                              step_cen_y_c,
    856                                              step_cen_y_d,
    857                                              step_cen_z_L,
    858                                              step_cen_z_U,
    859                                              step_cen_v_L,
    860                                              step_cen_v_U);
    861       if (qtmp < q) {
    862         sigma = sigma_lo;
    863         q = qtmp;
    864       }
    865     }
    866 
     850      Number q;
     851      if (qmid1 < qmid2) {
     852        sigma = sigma_mid1;
     853        q = qmid1;
     854      }
     855      else {
     856        sigma = sigma_mid2;
     857        q = qmid2;
     858      }
     859      if (sigma_up == ScaleSigma(sigma_up_in)) {
     860        Number qtmp;
     861        if (q_up<0.) {
     862          qtmp = CalculateQualityFunction(UnscaleSigma(sigma_up),
     863                                          step_aff_x_L,
     864                                          step_aff_x_U,
     865                                          step_aff_s_L,
     866                                          step_aff_s_U,
     867                                          step_aff_y_c,
     868                                          step_aff_y_d,
     869                                          step_aff_z_L,
     870                                          step_aff_z_U,
     871                                          step_aff_v_L,
     872                                          step_aff_v_U,
     873                                          step_cen_x_L,
     874                                          step_cen_x_U,
     875                                          step_cen_s_L,
     876                                          step_cen_s_U,
     877                                          step_cen_y_c,
     878                                          step_cen_y_d,
     879                                          step_cen_z_L,
     880                                          step_cen_z_U,
     881                                          step_cen_v_L,
     882                                          step_cen_v_U);
     883        }
     884        else {
     885          qtmp = q_up;
     886        }
     887        if (qtmp < q) {
     888          sigma = sigma_up;
     889          q = qtmp;
     890        }
     891      }
     892      else if (sigma_lo == ScaleSigma(sigma_lo_in)) {
     893        Number qtmp;
     894        if (q_lo<0.) {
     895          qtmp = CalculateQualityFunction(UnscaleSigma(sigma_lo),
     896                                          step_aff_x_L,
     897                                          step_aff_x_U,
     898                                          step_aff_s_L,
     899                                          step_aff_s_U,
     900                                          step_aff_y_c,
     901                                          step_aff_y_d,
     902                                          step_aff_z_L,
     903                                          step_aff_z_U,
     904                                          step_aff_v_L,
     905                                          step_aff_v_U,
     906                                          step_cen_x_L,
     907                                          step_cen_x_U,
     908                                          step_cen_s_L,
     909                                          step_cen_s_U,
     910                                          step_cen_y_c,
     911                                          step_cen_y_d,
     912                                          step_cen_z_L,
     913                                          step_cen_z_U,
     914                                          step_cen_v_L,
     915                                          step_cen_v_U);
     916        }
     917        else {
     918          qtmp = q_lo;
     919        }
     920        if (qtmp < q) {
     921          sigma = sigma_lo;
     922          q = qtmp;
     923        }
     924      }
     925    }
     926
     927    return UnscaleSigma(sigma);
     928  }
     929
     930  /*
     931  Number QualityFunctionMuOracle::ScaleSigma(Number sigma) {return log(sigma);}
     932  Number QualityFunctionMuOracle::UnscaleSigma(Number scaled_sigma) {return exp(scaled_sigma);}
     933  */
     934  Number QualityFunctionMuOracle::ScaleSigma(Number sigma)
     935  {
    867936    return sigma;
     937  }
     938  Number QualityFunctionMuOracle::UnscaleSigma(Number scaled_sigma)
     939  {
     940    return scaled_sigma;
    868941  }
    869942
     
    872945  /*
    873946  Number
    874   QualityFunctionMuOracle::PerformGoldenBisectionLog
     947  QualityFunctionMuOracle::PerformGoldenSectionLog
    875948  (Number sigma_up,
    876949   Number sigma_lo,
     
    9511024                                            step_cen_v_U);
    9521025
    953     Index nbisections = 0;
    954     while ((log_sigma_up-log_sigma_lo)>=tol*log_sigma_up && nbisections<max_bisection_steps_) {
    955       nbisections++;
     1026    Index nsections = 0;
     1027    while ((log_sigma_up-log_sigma_lo)>=tol*log_sigma_up && nsections<quality_function_max_section_steps_) {
     1028      nsections++;
    9561029      if (qmid1 > qmid2) {
    9571030        log_sigma_lo = log_sigma_mid1;
  • branches/dev/Algorithm/IpQualityFunctionMuOracle.hpp

    r501 r526  
    3838     *  could be used in the current iteration (using the LOQO formula).
    3939     */
    40     virtual Number CalculateMu();
     40    virtual Number CalculateMu(Number mu_min, Number mu_max);
    4141
    4242    /** Methods for IpoptType */
     
    122122                                    const Vector& step_cen_v_U);
    123123
    124     /** Auxilliary function performing the golden bisection */
    125     Number PerformGoldenBisection(Number sigma_up,
    126                                   Number sigma_lo,
    127                                   Number tol,
    128                                   const Vector& step_aff_x_L,
    129                                   const Vector& step_aff_x_U,
    130                                   const Vector& step_aff_s_L,
    131                                   const Vector& step_aff_s_U,
    132                                   const Vector& step_aff_y_c,
    133                                   const Vector& step_aff_y_d,
    134                                   const Vector& step_aff_z_L,
    135                                   const Vector& step_aff_z_U,
    136                                   const Vector& step_aff_v_L,
    137                                   const Vector& step_aff_v_U,
    138                                   const Vector& step_cen_x_L,
    139                                   const Vector& step_cen_x_U,
    140                                   const Vector& step_cen_s_L,
    141                                   const Vector& step_cen_s_U,
    142                                   const Vector& step_cen_y_c,
    143                                   const Vector& step_cen_y_d,
    144                                   const Vector& step_cen_z_L,
    145                                   const Vector& step_cen_z_U,
    146                                   const Vector& step_cen_v_L,
    147                                   const Vector& step_cen_v_U);
    148 
    149     /** Auxilliary function performing the golden bisection in the
     124    /** Auxilliary function performing the golden section */
     125    Number PerformGoldenSection(Number sigma_up,
     126                                Number q_up,
     127                                Number sigma_lo,
     128                                Number q_lo,
     129                                Number sigma_tol,
     130                                Number qf_tol,
     131                                const Vector& step_aff_x_L,
     132                                const Vector& step_aff_x_U,
     133                                const Vector& step_aff_s_L,
     134                                const Vector& step_aff_s_U,
     135                                const Vector& step_aff_y_c,
     136                                const Vector& step_aff_y_d,
     137                                const Vector& step_aff_z_L,
     138                                const Vector& step_aff_z_U,
     139                                const Vector& step_aff_v_L,
     140                                const Vector& step_aff_v_U,
     141                                const Vector& step_cen_x_L,
     142                                const Vector& step_cen_x_U,
     143                                const Vector& step_cen_s_L,
     144                                const Vector& step_cen_s_U,
     145                                const Vector& step_cen_y_c,
     146                                const Vector& step_cen_y_d,
     147                                const Vector& step_cen_z_L,
     148                                const Vector& step_cen_z_U,
     149                                const Vector& step_cen_v_L,
     150                                const Vector& step_cen_v_U);
     151
     152    /** Auxilliary functions for scaling the sigma axis in the golden
     153     *  section procedure */
     154    //@{
     155    Number ScaleSigma(Number sigma);
     156    Number UnscaleSigma(Number scaled_sigma);
     157    //@}
     158
     159    /** Auxilliary function performing the golden section in the
    150160     *  logarithmic scale */
    151161    /* This doesn't seem to work well, so I took it out for now (AW)
    152     Number PerformGoldenBisectionLog(Number sigma_up,
     162    Number PerformGoldenSectionLog(Number sigma_up,
    153163                                     Number sigma_lo,
    154164                                     Number tol,
     
    188198     */
    189199    BalancingTermEnum quality_function_balancing_term_;
    190     /** Relative tolerance for golden bi-section algorithm. */
    191     Number bisection_tol_;
    192     /** Maximal number of bi-section steps in the golden bisection
     200    /** Relative tolerance for golden bi-section algorithm in sigma
     201     *  space. */
     202    Number quality_function_section_sigma_tol_;
     203    /** Relative tolerance for golden bi-section algorithm in function
     204     *  value space. */
     205    Number quality_function_section_qf_tol_;
     206    /** Maximal number of bi-section steps in the golden section
    193207     *  search for sigma. */
    194     Index max_bisection_steps_;
     208    Index quality_function_max_section_steps_;
    195209    //@}
    196210
     
    216230    SmartPtr<Vector> tmp_v_U_;
    217231    //@}
     232
     233    // Deleteme
     234    Index count_qf_evals_;
    218235  };
    219236
  • branches/dev/Algorithm/IpRestoMinC_1Nrm.cpp

    r490 r526  
    109109    DBG_ASSERT(IpCq().curr_constraint_violation()>0.);
    110110
     111    // ToDo set those up during initialize?
    111112    // Create the restoration phase NLP etc objects
    112113    SmartPtr<IpoptData> resto_ip_data = new IpoptData();
  • branches/dev/Algorithm/IpTimingStatistics.cpp

    r524 r526  
    2020      return;
    2121
    22     jnlst.Printf(level, category, "\n\nTiming Statistics:\n\n");
    2322    jnlst.Printf(level, category,
    2423                 "OverallAlgorithm....................: %10.2f\n",
    2524                 OverallAlgorithm.TotalTime());
    2625    jnlst.Printf(level, category,
    27                  "FunctionEvaluations.................: %10.2f\n",
    28                  FunctionEvaluations.TotalTime());
     26                 " PrintProblemStatistics.............: %10.2f\n",
     27                 PrintProblemStatistics.TotalTime());
     28    jnlst.Printf(level, category,
     29                 " InitializeIterates.................: %10.2f\n",
     30                 InitializeIterates.TotalTime());
     31    jnlst.Printf(level, category,
     32                 " ActualizeHessian...................: %10.2f\n",
     33                 ActualizeHessian.TotalTime());
     34    jnlst.Printf(level, category,
     35                 " OutputIteration....................: %10.2f\n",
     36                 OutputIteration.TotalTime());
     37    jnlst.Printf(level, category,
     38                 " UpdateBarrierParameter.............: %10.2f\n",
     39                 UpdateBarrierParameter.TotalTime());
     40    jnlst.Printf(level, category,
     41                 " ComputeSearchDirection.............: %10.2f\n",
     42                 ComputeSearchDirection.TotalTime());
     43    jnlst.Printf(level, category,
     44                 " ComputeAcceptableTrialPoint........: %10.2f\n",
     45                 ComputeAcceptableTrialPoint.TotalTime());
     46    jnlst.Printf(level, category,
     47                 " AcceptTrialPoint...................: %10.2f\n",
     48                 AcceptTrialPoint.TotalTime());
     49    jnlst.Printf(level, category,
     50                 " CheckConvergence...................: %10.2f\n",
     51                 CheckConvergence.TotalTime());
     52
    2953    jnlst.Printf(level, category,
    3054                 "PDSystemSolverTotal.................: %10.2f\n",
    3155                 PDSystemSolverTotal.TotalTime());
    3256    jnlst.Printf(level, category,
    33                  "PDSystemSolverSolveOnce.............: %10.2f\n",
     57                 " PDSystemSolverSolveOnce............: %10.2f\n",
    3458                 PDSystemSolverSolveOnce.TotalTime());
    3559    jnlst.Printf(level, category,
    36                  "LinearSystemScaling.................: %10.2f\n",
     60                 " LinearSystemScaling................: %10.2f\n",
    3761                 LinearSystemScaling.TotalTime());
    3862    jnlst.Printf(level, category,
    39                  "LinearSystemSymbolicFactorization...: %10.2f\n",
     63                 " LinearSystemSymbolicFactorization..: %10.2f\n",
    4064                 LinearSystemSymbolicFactorization.TotalTime());
    4165    jnlst.Printf(level, category,
    42                  "LinearSystemFactorization...........: %10.2f\n",
     66                 " LinearSystemFactorization..........: %10.2f\n",
    4367                 LinearSystemFactorization.TotalTime());
    4468    jnlst.Printf(level, category,
    45                  "LinearSystemBackSolve...............: %10.2f\n",
     69                 " LinearSystemBackSolve..............: %10.2f\n",
    4670                 LinearSystemBackSolve.TotalTime());
    4771    jnlst.Printf(level, category,
  • branches/dev/Algorithm/IpTimingStatistics.hpp

    r524 r526  
    3838namespace Ipopt
    3939{
     40  /** This class is used to collect timing information for a
     41   *  particular task. */
     42  class TimedTask
     43  {
     44  public:
     45    /**@name Constructors/Destructors */
     46    //@{
     47    /** Default constructor. */
     48    TimedTask()
     49        :
     50        total_time_(0.),
     51        start_called_(false),
     52        end_called_(true)
     53    {}
     54
     55    /** Default destructor */
     56    ~TimedTask()
     57    {}
     58    //@}
     59
     60    /** Method that is called before execution of the task. */
     61    void Start()
     62    {
     63      DBG_ASSERT(end_called_);
     64      DBG_ASSERT(!start_called_);
     65      end_called_ = false;
     66      start_called_ = true;
     67      start_time_ = CpuTime();
     68    }
     69
     70    /** Method that is called after execution of the task. */
     71    void End()
     72    {
     73      DBG_ASSERT(!end_called_);
     74      DBG_ASSERT(start_called_);
     75      end_called_ = true;
     76      start_called_ = false;
     77      total_time_ += CpuTime() - start_time_;
     78    }
     79
     80    /** Method that is called after execution of the task for which
     81     *  timing might have been started.  This only updates the timing
     82     *  if the timing has indeed been conducted. This is useful to
     83     *  stop timing after catching exceptions. */
     84    void EndIfStarted()
     85    {
     86      if (start_called_) {
     87        end_called_ = true;
     88        start_called_ = false;
     89        total_time_ += CpuTime() - start_time_;
     90      }
     91      DBG_ASSERT(end_called_);
     92    }
     93
     94    /** Method returning total time spend for task so far. */
     95    Number TotalTime() const
     96    {
     97      DBG_ASSERT(end_called_);
     98      return total_time_;
     99    }
     100
     101  private:
     102    /**@name Default Compiler Generated Methods (Hidden to avoid
     103     * implicit creation/calling).  These methods are not
     104     * implemented and we do not want the compiler to implement them
     105     * for us, so we declare them private and do not define
     106     * them. This ensures that they will not be implicitly
     107     * created/called. */
     108    //@{
     109    /** Copy Constructor */
     110    TimedTask(const TimedTask&);
     111
     112    /** Overloaded Equals Operator */
     113    void operator=(const TimedTask&);
     114    //@}
     115
     116    /** Time at beginning of task. */
     117    Number start_time_;
     118    /** Total time for task measured so far. */
     119    Number total_time_;
     120
     121    /** @name fields for debugging */
     122    //@{
     123    bool start_called_;
     124    bool end_called_;
     125    //@}
     126
     127    // The following lines were taken from CoinTime.hpp in COIN/Coin
     128    /** method determining CPU executed since start of program */
     129    static inline Number CpuTime()
     130    {
     131      double cpu_temp;
     132#if defined(_MSC_VER)
     133
     134      unsigned int ticksnow;        /* clock_t is same as int */
     135
     136      ticksnow = (unsigned int)clock();
     137
     138      cpu_temp = (double)((double)ticksnow/CLOCKS_PER_SEC);
     139#else
     140
     141      struct rusage usage;
     142      getrusage(RUSAGE_SELF,&usage);
     143      cpu_temp = usage.ru_utime.tv_sec;
     144      cpu_temp += 1.0e-6*((double) usage.ru_utime.tv_usec);
     145#endif
     146
     147      return cpu_temp;
     148    }
     149  };
     150
    40151  /** This class collects all timing statistics for Ipopt.
    41152   */
    42153  class TimingStatistics : public ReferencedObject
    43154  {
    44     /** This class is used to collect timing information for a
    45      *  particular task. */
    46     class TimedTask
    47     {
    48     public:
    49       /**@name Constructors/Destructors */
    50       //@{
    51       /** Default constructor. */
    52       TimedTask()
    53           :
    54           total_time_(0.),
    55           start_called_(false),
    56           end_called_(true)
    57       {}
    58 
    59       /** Default destructor */
    60       ~TimedTask()
    61       {}
    62       //@}
    63 
    64       /** Method that is called before execution of the task. */
    65       void Start()
    66       {
    67         DBG_ASSERT(end_called_);
    68         DBG_ASSERT(!start_called_);
    69         DBG_DO(end_called_ = false);
    70         DBG_DO(start_called_ = true);
    71         start_time_ = CpuTime();
    72       }
    73 
    74       /** Method that is called after execution of the task. */
    75       void End()
    76       {
    77         DBG_ASSERT(!end_called_);
    78         DBG_ASSERT(start_called_);
    79         DBG_DO(end_called_ = true);
    80         DBG_DO(start_called_ = false);
    81         total_time_ += CpuTime() - start_time_;
    82       }
    83 
    84       /** Method returning total time spend for task so far. */
    85       Number TotalTime() const
    86       {
    87         DBG_ASSERT(end_called_);
    88         return total_time_;
    89       }
    90 
    91     private:
    92       /**@name Default Compiler Generated Methods (Hidden to avoid
    93        * implicit creation/calling).  These methods are not
    94        * implemented and we do not want the compiler to implement them
    95        * for us, so we declare them private and do not define
    96        * them. This ensures that they will not be implicitly
    97        * created/called. */
    98       //@{
    99       /** Copy Constructor */
    100       TimedTask(const TimedTask&);
    101 
    102       /** Overloaded Equals Operator */
    103       void operator=(const TimedTask&);
    104       //@}
    105 
    106       /** Time at beginning of task. */
    107       Number start_time_;
    108       /** Total time for task measured so far. */
    109       Number total_time_;
    110 
    111       /** @name fields for debugging */
    112       //@{
    113       bool start_called_;
    114       bool end_called_;
    115       //@}
    116 
    117       // The following lines were taken from CoinTime.hpp in COIN/Coin
    118       /** method determining CPU executed since start of program */
    119       static inline Number CpuTime()
    120       {
    121         double cpu_temp;
    122 #if defined(_MSC_VER)
    123 
    124         unsigned int ticksnow;        /* clock_t is same as int */
    125 
    126         ticksnow = (unsigned int)clock();
    127 
    128         cpu_temp = (double)((double)ticksnow/CLOCKS_PER_SEC);
    129 #else
    130 
    131         struct rusage usage;
    132         getrusage(RUSAGE_SELF,&usage);
    133         cpu_temp = usage.ru_utime.tv_sec;
    134         cpu_temp += 1.0e-6*((double) usage.ru_utime.tv_usec);
    135 #endif
    136 
    137         return cpu_temp;
    138       }
    139     };
    140 
    141155  public:
    142156    /**@name Constructors/Destructors */
     
    159173    //@{
    160174    TimedTask OverallAlgorithm;
    161     TimedTask FunctionEvaluations;
     175    TimedTask PrintProblemStatistics;
     176    TimedTask InitializeIterates;
     177    TimedTask ActualizeHessian;
     178    TimedTask OutputIteration;
     179    TimedTask UpdateBarrierParameter;
     180    TimedTask ComputeSearchDirection;
     181    TimedTask ComputeAcceptableTrialPoint;
     182    TimedTask AcceptTrialPoint;
     183    TimedTask CheckConvergence;
     184
    162185    TimedTask PDSystemSolverTotal;
    163186    TimedTask PDSystemSolverSolveOnce;
  • branches/dev/Interfaces/IpIpoptApplication.cpp

    r525 r526  
    406406      IpoptAlgorithm* p2alg = dynamic_cast<IpoptAlgorithm*> (GetRawPtr(alg_));
    407407      IpoptData* p2ip_data = dynamic_cast<IpoptData*> (GetRawPtr(ip_data_));
    408       IpoptNLP* p2ip_nlp = dynamic_cast<IpoptNLP*> (GetRawPtr(ip_nlp_));
     408      OrigIpoptNLP* p2ip_nlp = dynamic_cast<OrigIpoptNLP*> (GetRawPtr(ip_nlp_));
    409409      IpoptCalculatedQuantities* p2ip_cq = dynamic_cast<IpoptCalculatedQuantities*> (GetRawPtr(ip_cq_));
    410410      // Set up the algorithm
     
    490490                             print_timing_statistics, "");
    491491      if (print_timing_statistics) {
     492        jnlst_->Printf(J_SUMMARY, J_TIMING_STATISTICS,
     493                       "\n\nTiming Statistics:\n\n");
    492494        p2ip_data->TimingStats().PrintAllTimingStatistics(*jnlst_, J_SUMMARY,
    493495            J_TIMING_STATISTICS);
     496        p2ip_nlp->PrintTimingStatistics(*jnlst_, J_SUMMARY,
     497                                        J_TIMING_STATISTICS);
    494498      }
    495499
Note: See TracChangeset for help on using the changeset viewer.