1 // Copyright (C) 2004, 2009 International Business Machines and others. 2 // All Rights Reserved. 3 // This code is published under the Eclipse Public License. 4 // 5 // $Id$ 6 // 7 // Authors: Carl Laird, Andreas Waechter IBM 2004-11-12 8 9 #include "IpQualityFunctionMuOracle.hpp" 10 11 #ifdef HAVE_CMATH 12 # include <cmath> 13 #else 14 # ifdef HAVE_MATH_H 15 # include <math.h> 16 # else 17 # error "don't have header file for math" 18 # endif 19 #endif 20 21 #ifdef HAVE_CSTDIO 22 # include <cstdio> 23 #else 24 # ifdef HAVE_STDIO_H 25 # include <stdio.h> 26 # else 27 # error "don't have header file for stdio" 28 # endif 29 #endif 30 31 namespace Ipopt 32 { 33 #if COIN_IPOPT_VERBOSITY > 0 34 static const Index dbg_verbosity = 0; 35 #endif 36 QualityFunctionMuOracle(const SmartPtr<PDSystemSolver> & pd_solver)37 QualityFunctionMuOracle::QualityFunctionMuOracle(const SmartPtr<PDSystemSolver>& pd_solver) 38 : 39 MuOracle(), 40 pd_solver_(pd_solver), 41 42 tmp_step_x_L_(NULL), 43 tmp_step_x_U_(NULL), 44 tmp_step_s_L_(NULL), 45 tmp_step_s_U_(NULL), 46 tmp_step_z_L_(NULL), 47 tmp_step_z_U_(NULL), 48 tmp_step_v_L_(NULL), 49 tmp_step_v_U_(NULL), 50 51 tmp_slack_x_L_(NULL), 52 tmp_slack_x_U_(NULL), 53 tmp_slack_s_L_(NULL), 54 tmp_slack_s_U_(NULL), 55 tmp_z_L_(NULL), 56 tmp_z_U_(NULL), 57 tmp_v_L_(NULL), 58 tmp_v_U_(NULL), 59 60 count_qf_evals_(0) 61 { 62 DBG_ASSERT(IsValid(pd_solver_)); 63 } 64 ~QualityFunctionMuOracle()65 QualityFunctionMuOracle::~QualityFunctionMuOracle() 66 {} 67 RegisterOptions(SmartPtr<RegisteredOptions> roptions)68 void QualityFunctionMuOracle::RegisterOptions(SmartPtr<RegisteredOptions> roptions) 69 { 70 roptions->AddLowerBoundedNumberOption( 71 "sigma_max", 72 "Maximum value of the centering parameter.", 73 0.0, true, 1e2, 74 "This is the upper bound for the centering parameter chosen by the " 75 "quality function based barrier parameter update. (Only used if option " 76 "\"mu_oracle\" is set to \"quality-function\".)"); 77 roptions->AddLowerBoundedNumberOption( 78 "sigma_min", 79 "Minimum value of the centering parameter.", 80 0.0, false, 1e-6, 81 "This is the lower bound for the centering parameter chosen by the " 82 "quality function based barrier parameter update. (Only used if option " 83 "\"mu_oracle\" is set to \"quality-function\".)"); 84 roptions->AddStringOption4( 85 "quality_function_norm_type", 86 "Norm used for components of the quality function.", 87 "2-norm-squared", 88 "1-norm", "use the 1-norm (abs sum)", 89 "2-norm-squared", "use the 2-norm squared (sum of squares)", 90 "max-norm", "use the infinity norm (max)", 91 "2-norm", "use 2-norm", 92 "(Only used if option \"mu_oracle\" is set to \"quality-function\".)"); 93 roptions->AddStringOption4( 94 "quality_function_centrality", 95 "The penalty term for centrality that is included in quality function.", 96 "none", 97 "none", "no penalty term is added", 98 "log", "complementarity * the log of the centrality measure", 99 "reciprocal", "complementarity * the reciprocal of the centrality measure", 100 "cubed-reciprocal", "complementarity * the reciprocal of the centrality measure cubed", 101 "This determines whether a term is added to the quality function to " 102 "penalize deviation from centrality with respect to complementarity. The " 103 "complementarity measure here is the xi in the Loqo update rule. (Only used if option " 104 "\"mu_oracle\" is set to \"quality-function\".)"); 105 roptions->AddStringOption2( 106 "quality_function_balancing_term", 107 "The balancing term included in the quality function for centrality.", 108 "none", 109 "none", "no balancing term is added", 110 "cubic", "Max(0,Max(dual_inf,primal_inf)-compl)^3", 111 "This determines whether a term is added to the quality function that " 112 "penalizes situations where the complementarity is much smaller " 113 "than dual and primal infeasibilities. (Only used if option " 114 "\"mu_oracle\" is set to \"quality-function\".)"); 115 roptions->AddLowerBoundedIntegerOption( 116 "quality_function_max_section_steps", 117 "Maximum number of search steps during direct search procedure " 118 "determining the optimal centering parameter.", 119 0, 8, 120 "The golden section search is performed for the quality function based " 121 "mu oracle. (Only used if option " 122 "\"mu_oracle\" is set to \"quality-function\".)"); 123 roptions->AddBoundedNumberOption( 124 "quality_function_section_sigma_tol", 125 "Tolerance for the section search procedure determining " 126 "the optimal centering parameter (in sigma space).", 127 0.0, false, 1.0, true, 128 1e-2, 129 "The golden section search is performed for the quality function based " 130 "mu oracle. (Only used if option " 131 "\"mu_oracle\" is set to \"quality-function\".)"); 132 roptions->AddBoundedNumberOption( 133 "quality_function_section_qf_tol", 134 "Tolerance for the golden section search procedure determining " 135 "the optimal centering parameter (in the function value space).", 136 0.0, false, 1.0, true, 137 0e-2, 138 "The golden section search is performed for the quality function based mu " 139 "oracle. (Only used if option " 140 "\"mu_oracle\" is set to \"quality-function\".)"); 141 } 142 143 InitializeImpl(const OptionsList & options,const std::string & prefix)144 bool QualityFunctionMuOracle::InitializeImpl(const OptionsList& options, 145 const std::string& prefix) 146 { 147 Index enum_int; 148 149 options.GetNumericValue("sigma_max", sigma_max_, prefix); 150 options.GetNumericValue("sigma_min", sigma_min_, prefix); 151 152 options.GetEnumValue("quality_function_norm_type", enum_int, prefix); 153 quality_function_norm_ = NormEnum(enum_int); 154 options.GetEnumValue("quality_function_centrality", enum_int, prefix); 155 quality_function_centrality_ = CentralityEnum(enum_int); 156 options.GetEnumValue("quality_function_balancing_term", enum_int, prefix); 157 quality_function_balancing_term_ = BalancingTermEnum(enum_int); 158 options.GetIntegerValue("quality_function_max_section_steps", 159 quality_function_max_section_steps_, prefix); 160 options.GetNumericValue("quality_function_section_sigma_tol", 161 quality_function_section_sigma_tol_, prefix); 162 options.GetNumericValue("quality_function_section_qf_tol", 163 quality_function_section_qf_tol_, prefix); 164 165 initialized_ = false; 166 167 return true; 168 } 169 CalculateMu(Number mu_min,Number mu_max,Number & new_mu)170 bool QualityFunctionMuOracle::CalculateMu(Number mu_min, Number mu_max, 171 Number& new_mu) 172 { 173 DBG_START_METH("QualityFunctionMuOracle::CalculateMu", 174 dbg_verbosity); 175 176 /////////////////////////////////////////////////////////////////////////// 177 // Reserve memory for temporary vectors used in CalculateQualityFunction // 178 /////////////////////////////////////////////////////////////////////////// 179 180 tmp_step_x_L_ = IpNLP().x_L()->MakeNew(); 181 tmp_step_x_U_ = IpNLP().x_U()->MakeNew(); 182 tmp_step_s_L_ = IpNLP().d_L()->MakeNew(); 183 tmp_step_s_U_ = IpNLP().d_U()->MakeNew(); 184 tmp_step_z_L_ = IpNLP().x_L()->MakeNew(); 185 tmp_step_z_U_ = IpNLP().x_U()->MakeNew(); 186 tmp_step_v_L_ = IpNLP().d_L()->MakeNew(); 187 tmp_step_v_U_ = IpNLP().d_U()->MakeNew(); 188 189 tmp_slack_x_L_ = IpNLP().x_L()->MakeNew(); 190 tmp_slack_x_U_ = IpNLP().x_U()->MakeNew(); 191 tmp_slack_s_L_ = IpNLP().d_L()->MakeNew(); 192 tmp_slack_s_U_ = IpNLP().d_U()->MakeNew(); 193 tmp_z_L_ = IpNLP().x_L()->MakeNew(); 194 tmp_z_U_ = IpNLP().x_U()->MakeNew(); 195 tmp_v_L_ = IpNLP().d_L()->MakeNew(); 196 tmp_v_U_ = IpNLP().d_U()->MakeNew(); 197 198 ///////////////////////////////////// 199 // Compute the affine scaling step // 200 ///////////////////////////////////// 201 202 Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE, 203 "Solving the Primal Dual System for the affine step\n"); 204 // First get the right hand side 205 SmartPtr<IteratesVector> rhs_aff = IpData().curr()->MakeNewIteratesVector(false); 206 rhs_aff->Set_x(*IpCq().curr_grad_lag_x()); 207 rhs_aff->Set_s(*IpCq().curr_grad_lag_s()); 208 rhs_aff->Set_y_c(*IpCq().curr_c()); 209 rhs_aff->Set_y_d(*IpCq().curr_d_minus_s()); 210 rhs_aff->Set_z_L(*IpCq().curr_compl_x_L()); 211 rhs_aff->Set_z_U(*IpCq().curr_compl_x_U()); 212 rhs_aff->Set_v_L(*IpCq().curr_compl_s_L()); 213 rhs_aff->Set_v_U(*IpCq().curr_compl_s_U()); 214 215 // Get space for the affine scaling step 216 SmartPtr<IteratesVector> step_aff = IpData().curr()->MakeNewIteratesVector(true); 217 218 // Now solve the primal-dual system to get the step. We allow a 219 // somewhat inexact solution, iterative refinement will be done 220 // after mu is known 221 bool allow_inexact = true; 222 bool retval = pd_solver_->Solve(-1.0, 0.0, 223 *rhs_aff, 224 *step_aff, 225 allow_inexact 226 ); 227 if (!retval) { 228 Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE, 229 "The linear system could not be solved for the affine step!\n"); 230 return false; 231 } 232 233 DBG_PRINT_VECTOR(2, "step_aff", *step_aff); 234 235 ///////////////////////////////////// 236 // Compute the pure centering step // 237 ///////////////////////////////////// 238 239 Number avrg_compl = IpCq().curr_avrg_compl(); 240 241 Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE, 242 "Solving the Primal Dual System for the centering step\n"); 243 // First get the right hand side 244 SmartPtr<IteratesVector> rhs_cen = IpData().curr()->MakeNewIteratesVector(true); 245 rhs_cen->x_NonConst()->AddOneVector(-avrg_compl, 246 *IpCq().grad_kappa_times_damping_x(), 247 0.); 248 rhs_cen->s_NonConst()->AddOneVector(-avrg_compl, 249 *IpCq().grad_kappa_times_damping_s(), 250 0.); 251 252 rhs_cen->y_c_NonConst()->Set(0.); 253 rhs_cen->y_d_NonConst()->Set(0.); 254 rhs_cen->z_L_NonConst()->Set(avrg_compl); 255 rhs_cen->z_U_NonConst()->Set(avrg_compl); 256 rhs_cen->v_L_NonConst()->Set(avrg_compl); 257 rhs_cen->v_U_NonConst()->Set(avrg_compl); 258 259 // Get space for the centering step 260 SmartPtr<IteratesVector> step_cen = IpData().curr()->MakeNewIteratesVector(true); 261 262 // Now solve the primal-dual system to get the step 263 allow_inexact = true; 264 retval = pd_solver_->Solve(1.0, 0.0, 265 *rhs_cen, 266 *step_cen, 267 allow_inexact 268 ); 269 if (!retval) { 270 Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE, 271 "The linear system could not be solved for the centering step!\n"); 272 return false; 273 } 274 275 DBG_PRINT_VECTOR(2, "step_cen", *step_cen); 276 277 // Start the timing for the quality function search here 278 IpData().TimingStats().QualityFunctionSearch().Start(); 279 280 // Some initializations 281 if (!initialized_) { 282 n_dual_ = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim(); 283 n_pri_ = IpData().curr()->y_c()->Dim() + IpData().curr()->y_d()->Dim(); 284 n_comp_ = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim() + 285 IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim(); 286 initialized_ = true; 287 } 288 289 count_qf_evals_ = 0; 290 291 // Compute some quantities used for the quality function evaluations 292 // (This way we try to avoid retrieving numbers from cache... 293 294 curr_slack_x_L_ = IpCq().curr_slack_x_L(); 295 curr_slack_x_U_ = IpCq().curr_slack_x_U(); 296 curr_slack_s_L_ = IpCq().curr_slack_s_L(); 297 curr_slack_s_U_ = IpCq().curr_slack_s_U(); 298 299 curr_z_L_ = IpData().curr()->z_L(); 300 curr_z_U_ = IpData().curr()->z_U(); 301 curr_v_L_ = IpData().curr()->v_L(); 302 curr_v_U_ = IpData().curr()->v_U(); 303 304 IpData().TimingStats().Task5().Start(); 305 switch (quality_function_norm_) { 306 case NM_NORM_1: 307 curr_grad_lag_x_asum_ = IpCq().curr_grad_lag_x()->Asum(); 308 curr_grad_lag_s_asum_ = IpCq().curr_grad_lag_s()->Asum(); 309 curr_c_asum_ = IpCq().curr_c()->Asum(); 310 curr_d_minus_s_asum_ = IpCq().curr_d_minus_s()->Asum(); 311 break; 312 case NM_NORM_2_SQUARED: 313 case NM_NORM_2: 314 curr_grad_lag_x_nrm2_ = IpCq().curr_grad_lag_x()->Nrm2(); 315 curr_grad_lag_s_nrm2_ = IpCq().curr_grad_lag_s()->Nrm2(); 316 curr_c_nrm2_ = IpCq().curr_c()->Nrm2(); 317 curr_d_minus_s_nrm2_ = IpCq().curr_d_minus_s()->Nrm2(); 318 break; 319 case NM_NORM_MAX: 320 curr_grad_lag_x_amax_ = IpCq().curr_grad_lag_x()->Amax(); 321 curr_grad_lag_s_amax_ = IpCq().curr_grad_lag_s()->Amax(); 322 curr_c_amax_ = IpCq().curr_c()->Amax(); 323 curr_d_minus_s_amax_ = IpCq().curr_d_minus_s()->Amax(); 324 break; 325 default: 326 DBG_ASSERT(false && "Unknown value for quality_function_norm_"); 327 } 328 IpData().TimingStats().Task5().End(); 329 330 // We now compute the step for the slack variables. This safes 331 // time, because we then don't have to do this any more for each 332 // evaluation of the quality function 333 SmartPtr<Vector> step_aff_x_L = step_aff->z_L()->MakeNew(); 334 SmartPtr<Vector> step_aff_x_U = step_aff->z_U()->MakeNew(); 335 SmartPtr<Vector> step_aff_s_L = step_aff->v_L()->MakeNew(); 336 SmartPtr<Vector> step_aff_s_U = step_aff->v_U()->MakeNew(); 337 IpNLP().Px_L()->TransMultVector(1., *step_aff->x(), 0., *step_aff_x_L); 338 IpNLP().Px_U()->TransMultVector(-1., *step_aff->x(), 0., *step_aff_x_U); 339 IpNLP().Pd_L()->TransMultVector(1., *step_aff->s(), 0., *step_aff_s_L); 340 IpNLP().Pd_U()->TransMultVector(-1., *step_aff->s(), 0., *step_aff_s_U); 341 SmartPtr<Vector> step_cen_x_L = step_cen->z_L()->MakeNew(); 342 SmartPtr<Vector> step_cen_x_U = step_cen->z_U()->MakeNew(); 343 SmartPtr<Vector> step_cen_s_L = step_cen->v_L()->MakeNew(); 344 SmartPtr<Vector> step_cen_s_U = step_cen->v_U()->MakeNew(); 345 IpNLP().Px_L()->TransMultVector(1., *step_cen->x(), 0., *step_cen_x_L); 346 IpNLP().Px_U()->TransMultVector(-1., *step_cen->x(), 0., *step_cen_x_U); 347 IpNLP().Pd_L()->TransMultVector(1., *step_cen->s(), 0., *step_cen_s_L); 348 IpNLP().Pd_U()->TransMultVector(-1., *step_cen->s(), 0., *step_cen_s_U); 349 350 Number sigma; 351 352 // First we determine whether we want to search for a value of 353 // sigma larger or smaller than 1. For this, we estimate the 354 // slope of the quality function at sigma=1. 355 Number qf_1 = CalculateQualityFunction(1., 356 *step_aff_x_L, 357 *step_aff_x_U, 358 *step_aff_s_L, 359 *step_aff_s_U, 360 *step_aff->y_c(), 361 *step_aff->y_d(), 362 *step_aff->z_L(), 363 *step_aff->z_U(), 364 *step_aff->v_L(), 365 *step_aff->v_U(), 366 *step_cen_x_L, 367 *step_cen_x_U, 368 *step_cen_s_L, 369 *step_cen_s_U, 370 *step_cen->y_c(), 371 *step_cen->y_d(), 372 *step_cen->z_L(), 373 *step_cen->z_U(), 374 *step_cen->v_L(), 375 *step_cen->v_U()); 376 377 Number sigma_1minus = 1.-Max(1e-4, quality_function_section_sigma_tol_); 378 Number qf_1minus = CalculateQualityFunction(sigma_1minus, 379 *step_aff_x_L, 380 *step_aff_x_U, 381 *step_aff_s_L, 382 *step_aff_s_U, 383 *step_aff->y_c(), 384 *step_aff->y_d(), 385 *step_aff->z_L(), 386 *step_aff->z_U(), 387 *step_aff->v_L(), 388 *step_aff->v_U(), 389 *step_cen_x_L, 390 *step_cen_x_U, 391 *step_cen_s_L, 392 *step_cen_s_U, 393 *step_cen->y_c(), 394 *step_cen->y_d(), 395 *step_cen->z_L(), 396 *step_cen->z_U(), 397 *step_cen->v_L(), 398 *step_cen->v_U()); 399 400 if (qf_1minus > qf_1) { 401 // It seems that the quality function decreases for values 402 // larger than sigma, so perform golden section search for sigma 403 // > 1. 404 Number sigma_up = Min(sigma_max_, mu_max/avrg_compl); 405 Number sigma_lo = 1.; 406 if (sigma_lo >= sigma_up) { 407 sigma = sigma_up; 408 } 409 else { 410 // ToDo maybe we should use different tolerances for sigma>1 411 sigma = PerformGoldenSection(sigma_up, -100., sigma_lo, qf_1, 412 quality_function_section_sigma_tol_, 413 quality_function_section_qf_tol_, 414 *step_aff_x_L, 415 *step_aff_x_U, 416 *step_aff_s_L, 417 *step_aff_s_U, 418 *step_aff->y_c(), 419 *step_aff->y_d(), 420 *step_aff->z_L(), 421 *step_aff->z_U(), 422 *step_aff->v_L(), 423 *step_aff->v_U(), 424 *step_cen_x_L, 425 *step_cen_x_U, 426 *step_cen_s_L, 427 *step_cen_s_U, 428 *step_cen->y_c(), 429 *step_cen->y_d(), 430 *step_cen->z_L(), 431 *step_cen->z_U(), 432 *step_cen->v_L(), 433 *step_cen->v_U()); 434 } 435 } 436 else { 437 // Search for sigma less than 1 438 439 Number sigma_lo = Max(sigma_min_, mu_min/avrg_compl); 440 Number sigma_up = Min(Max(sigma_lo, sigma_1minus), mu_max/avrg_compl); 441 if (sigma_lo >= sigma_up) { 442 // Skip the search, we are already at the minimum 443 sigma = sigma_lo; 444 } 445 else { 446 sigma = PerformGoldenSection(sigma_up, qf_1minus, sigma_lo, -100., 447 quality_function_section_sigma_tol_, 448 quality_function_section_qf_tol_, 449 *step_aff_x_L, 450 *step_aff_x_U, 451 *step_aff_s_L, 452 *step_aff_s_U, 453 *step_aff->y_c(), 454 *step_aff->y_d(), 455 *step_aff->z_L(), 456 *step_aff->z_U(), 457 *step_aff->v_L(), 458 *step_aff->v_U(), 459 *step_cen_x_L, 460 *step_cen_x_U, 461 *step_cen_s_L, 462 *step_cen_s_U, 463 *step_cen->y_c(), 464 *step_cen->y_d(), 465 *step_cen->z_L(), 466 *step_cen->z_U(), 467 *step_cen->v_L(), 468 *step_cen->v_U()); 469 } 470 } 471 472 //#define tracequalityfunction 473 #ifdef tracequalityfunction 474 char fname[100]; 475 Snprintf(fname, 100, "qf_values_%d.dat", IpData().iter_count()); 476 FILE* fid = fopen(fname, "w"); 477 478 Number sigma_1 = sigma_max_; 479 Number sigma_2 = 1e-9/avrg_compl; 480 Number sigma_trace = sigma_1; 481 while (sigma_trace > sigma_2) { 482 Number qf = CalculateQualityFunction(sigma_trace, 483 *step_aff_x_L, 484 *step_aff_x_U, 485 *step_aff_s_L, 486 *step_aff_s_U, 487 *step_aff->y_c(), 488 *step_aff->y_d(), 489 *step_aff->z_L(), 490 *step_aff->z_U(), 491 *step_aff->v_L(), 492 *step_aff->v_U(), 493 *step_cen_x_L, 494 *step_cen_x_U, 495 *step_cen_s_L, 496 *step_cen_s_U, 497 *step_cen->y_c(), 498 *step_cen->y_d(), 499 *step_cen->z_L(), 500 *step_cen->z_U(), 501 *step_cen->v_L(), 502 *step_cen->v_U()); 503 fprintf(fid, "%9.2e %25.16e\n", sigma_trace, qf); 504 sigma_trace /= 1.1; 505 } 506 fclose(fid); 507 #endif 508 509 // End timing of quality function search 510 IpData().TimingStats().QualityFunctionSearch().End(); 511 512 Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE, 513 "Sigma = %e\n", sigma); 514 Number mu = sigma*avrg_compl; 515 516 // Store the affine search direction (in case it is needed in the 517 // line search for a corrector step) 518 IpData().set_delta_aff(step_aff); 519 IpData().SetHaveAffineDeltas(true); 520 521 // Now construct the overall search direction here 522 SmartPtr<IteratesVector> step = IpData().curr()->MakeNewIteratesVector(true); 523 step->AddTwoVectors(sigma, *step_cen, 1.0, *IpData().delta_aff(), 0.0); 524 525 DBG_PRINT_VECTOR(2, "step", *step); 526 IpData().set_delta(step); 527 IpData().SetHaveDeltas(true); 528 529 /////////////////////////////////////////////////////////////////////////// 530 // Release memory for temporary vectors used in CalculateQualityFunction // 531 /////////////////////////////////////////////////////////////////////////// 532 533 tmp_step_x_L_ = NULL; 534 tmp_step_x_U_ = NULL; 535 tmp_step_s_L_ = NULL; 536 tmp_step_s_U_ = NULL; 537 tmp_step_z_L_ = NULL; 538 tmp_step_z_U_ = NULL; 539 tmp_step_v_L_ = NULL; 540 tmp_step_v_U_ = NULL; 541 542 tmp_slack_x_L_ = NULL; 543 tmp_slack_x_U_ = NULL; 544 tmp_slack_s_L_ = NULL; 545 tmp_slack_s_U_ = NULL; 546 tmp_z_L_ = NULL; 547 tmp_z_U_ = NULL; 548 tmp_v_L_ = NULL; 549 tmp_v_U_ = NULL; 550 551 curr_slack_x_L_ = NULL; 552 curr_slack_x_U_ = NULL; 553 curr_slack_s_L_ = NULL; 554 curr_slack_s_U_ = NULL; 555 556 // DELETEME 557 char ssigma[40]; 558 Snprintf(ssigma, 39, " sigma=%8.2e", sigma); 559 IpData().Append_info_string(ssigma); 560 Snprintf(ssigma, 39, " qf=%d", count_qf_evals_); 561 IpData().Append_info_string(ssigma); 562 /* 563 Snprintf(ssigma, 39, " xi=%8.2e ", IpCq().curr_centrality_measure()); 564 IpData().Append_info_string(ssigma); 565 if (sigma>1.) { 566 IpData().Append_info_string("LARGESIGMA"); 567 } 568 */ 569 570 new_mu = mu; 571 return true; 572 } 573 CalculateQualityFunction(Number sigma,const Vector & step_aff_x_L,const Vector & step_aff_x_U,const Vector & step_aff_s_L,const Vector & step_aff_s_U,const Vector & step_aff_y_c,const Vector & step_aff_y_d,const Vector & step_aff_z_L,const Vector & step_aff_z_U,const Vector & step_aff_v_L,const Vector & step_aff_v_U,const Vector & step_cen_x_L,const Vector & step_cen_x_U,const Vector & step_cen_s_L,const Vector & step_cen_s_U,const Vector & step_cen_y_c,const Vector & step_cen_y_d,const Vector & step_cen_z_L,const Vector & step_cen_z_U,const Vector & step_cen_v_L,const Vector & step_cen_v_U)574 Number QualityFunctionMuOracle::CalculateQualityFunction 575 (Number sigma, 576 const Vector& step_aff_x_L, 577 const Vector& step_aff_x_U, 578 const Vector& step_aff_s_L, 579 const Vector& step_aff_s_U, 580 const Vector& step_aff_y_c, 581 const Vector& step_aff_y_d, 582 const Vector& step_aff_z_L, 583 const Vector& step_aff_z_U, 584 const Vector& step_aff_v_L, 585 const Vector& step_aff_v_U, 586 const Vector& step_cen_x_L, 587 const Vector& step_cen_x_U, 588 const Vector& step_cen_s_L, 589 const Vector& step_cen_s_U, 590 const Vector& step_cen_y_c, 591 const Vector& step_cen_y_d, 592 const Vector& step_cen_z_L, 593 const Vector& step_cen_z_U, 594 const Vector& step_cen_v_L, 595 const Vector& step_cen_v_U 596 ) 597 { 598 DBG_START_METH("QualityFunctionMuOracle::CalculateQualityFunction", 599 dbg_verbosity); 600 count_qf_evals_++; 601 602 IpData().TimingStats().Task1().Start(); 603 tmp_step_x_L_->AddTwoVectors(1., step_aff_x_L, sigma, step_cen_x_L, 0.); 604 tmp_step_x_U_->AddTwoVectors(1., step_aff_x_U, sigma, step_cen_x_U, 0.); 605 tmp_step_s_L_->AddTwoVectors(1., step_aff_s_L, sigma, step_cen_s_L, 0.); 606 tmp_step_s_U_->AddTwoVectors(1., step_aff_s_U, sigma, step_cen_s_U, 0.); 607 tmp_step_z_L_->AddTwoVectors(1., step_aff_z_L, sigma, step_cen_z_L, 0.); 608 tmp_step_z_U_->AddTwoVectors(1., step_aff_z_U, sigma, step_cen_z_U, 0.); 609 tmp_step_v_L_->AddTwoVectors(1., step_aff_v_L, sigma, step_cen_v_L, 0.); 610 tmp_step_v_U_->AddTwoVectors(1., step_aff_v_U, sigma, step_cen_v_U, 0.); 611 IpData().TimingStats().Task1().End(); 612 613 // Compute the fraction-to-the-boundary step sizes 614 IpData().TimingStats().Task2().Start(); 615 Number tau = IpData().curr_tau(); 616 Number alpha_primal = IpCq().uncached_slack_frac_to_the_bound(tau, 617 *tmp_step_x_L_, 618 *tmp_step_x_U_, 619 *tmp_step_s_L_, 620 *tmp_step_s_U_); 621 622 Number alpha_dual = IpCq().uncached_dual_frac_to_the_bound(tau, 623 *tmp_step_z_L_, 624 *tmp_step_z_U_, 625 *tmp_step_v_L_, 626 *tmp_step_v_U_); 627 IpData().TimingStats().Task2().End(); 628 629 Number xi = 0.; // centrality measure 630 631 IpData().TimingStats().Task1().Start(); 632 tmp_slack_x_L_->AddTwoVectors(1., *curr_slack_x_L_, 633 alpha_primal, *tmp_step_x_L_, 0.); 634 tmp_slack_x_U_->AddTwoVectors(1., *curr_slack_x_U_, 635 alpha_primal, *tmp_step_x_U_, 0.); 636 tmp_slack_s_L_->AddTwoVectors(1., *curr_slack_s_L_, 637 alpha_primal, *tmp_step_s_L_, 0.); 638 tmp_slack_s_U_->AddTwoVectors(1., *curr_slack_s_U_, 639 alpha_primal, *tmp_step_s_U_, 0.); 640 641 tmp_z_L_->AddTwoVectors(1., *curr_z_L_, 642 alpha_dual, *tmp_step_z_L_, 0.); 643 tmp_z_U_->AddTwoVectors(1., *curr_z_U_, 644 alpha_dual, *tmp_step_z_U_, 0.); 645 tmp_v_L_->AddTwoVectors(1., *curr_v_L_, 646 alpha_dual, *tmp_step_v_L_, 0.); 647 tmp_v_U_->AddTwoVectors(1., *curr_v_U_, 648 alpha_dual, *tmp_step_v_U_, 0.); 649 IpData().TimingStats().Task1().End(); 650 651 IpData().TimingStats().Task3().Start(); 652 tmp_slack_x_L_->ElementWiseMultiply(*tmp_z_L_); 653 tmp_slack_x_U_->ElementWiseMultiply(*tmp_z_U_); 654 tmp_slack_s_L_->ElementWiseMultiply(*tmp_v_L_); 655 tmp_slack_s_U_->ElementWiseMultiply(*tmp_v_U_); 656 IpData().TimingStats().Task3().End(); 657 658 DBG_PRINT_VECTOR(2, "compl_x_L", *tmp_slack_x_L_); 659 DBG_PRINT_VECTOR(2, "compl_x_U", *tmp_slack_x_U_); 660 DBG_PRINT_VECTOR(2, "compl_s_L", *tmp_slack_s_L_); 661 DBG_PRINT_VECTOR(2, "compl_s_U", *tmp_slack_s_U_); 662 663 Number dual_inf=-1.; 664 Number primal_inf=-1.; 665 Number compl_inf=-1.; 666 667 IpData().TimingStats().Task5().Start(); 668 switch (quality_function_norm_) { 669 case NM_NORM_1: 670 dual_inf = (1.-alpha_dual)*(curr_grad_lag_x_asum_ + 671 curr_grad_lag_s_asum_); 672 673 primal_inf = (1.-alpha_primal)*(curr_c_asum_ + 674 curr_d_minus_s_asum_); 675 676 compl_inf = tmp_slack_x_L_->Asum() + tmp_slack_x_U_->Asum() + 677 tmp_slack_s_L_->Asum() + tmp_slack_s_U_->Asum(); 678 679 dual_inf /= n_dual_; 680 if (n_pri_>0) { 681 primal_inf /= n_pri_; 682 } 683 DBG_ASSERT(n_comp_>0); 684 compl_inf /= n_comp_; 685 break; 686 case NM_NORM_2_SQUARED: 687 dual_inf = 688 pow(1.-alpha_dual, 2)*(pow(curr_grad_lag_x_nrm2_, 2) + 689 pow(curr_grad_lag_s_nrm2_, 2)); 690 primal_inf = 691 pow(1.-alpha_primal, 2)*(pow(curr_c_nrm2_, 2) + 692 pow(curr_d_minus_s_nrm2_, 2)); 693 compl_inf = 694 pow(tmp_slack_x_L_->Nrm2(), 2) + pow(tmp_slack_x_U_->Nrm2(), 2) + 695 pow(tmp_slack_s_L_->Nrm2(), 2) + pow(tmp_slack_s_U_->Nrm2(), 2); 696 697 dual_inf /= n_dual_; 698 if (n_pri_>0) { 699 primal_inf /= n_pri_; 700 } 701 DBG_ASSERT(n_comp_>0); 702 compl_inf /= n_comp_; 703 break; 704 case NM_NORM_MAX: 705 dual_inf = 706 (1.-alpha_dual)*Max(curr_grad_lag_x_amax_, 707 curr_grad_lag_s_amax_); 708 primal_inf = 709 (1.-alpha_primal)*Max(curr_c_amax_, 710 curr_d_minus_s_amax_); 711 compl_inf = 712 Max(tmp_slack_x_L_->Amax(), tmp_slack_x_U_->Amax(), 713 tmp_slack_s_L_->Amax(), tmp_slack_s_U_->Amax()); 714 break; 715 case NM_NORM_2: 716 dual_inf = 717 (1.-alpha_dual)*sqrt(pow(curr_grad_lag_x_nrm2_, 2) + 718 pow(curr_grad_lag_s_nrm2_, 2)); 719 primal_inf = 720 (1.-alpha_primal)*sqrt(pow(curr_c_nrm2_, 2) + 721 pow(curr_d_minus_s_nrm2_, 2)); 722 compl_inf = 723 sqrt(pow(tmp_slack_x_L_->Nrm2(), 2) + pow(tmp_slack_x_U_->Nrm2(), 2) + 724 pow(tmp_slack_s_L_->Nrm2(), 2) + pow(tmp_slack_s_U_->Nrm2(), 2)); 725 726 dual_inf /= sqrt((Number)n_dual_); 727 if (n_pri_>0) { 728 primal_inf /= sqrt((Number)n_pri_); 729 } 730 DBG_ASSERT(n_comp_>0); 731 compl_inf /= sqrt((Number)n_comp_); 732 break; 733 default: 734 DBG_ASSERT(false && "Unknown value for quality_function_norm_"); 735 } 736 IpData().TimingStats().Task5().End(); 737 738 Number quality_function = dual_inf + primal_inf + compl_inf; 739 740 if (quality_function_centrality_!=CEN_NONE) { 741 IpData().TimingStats().Task4().Start(); 742 xi = IpCq().CalcCentralityMeasure(*tmp_slack_x_L_, *tmp_slack_x_U_, 743 *tmp_slack_s_L_, *tmp_slack_s_U_); 744 IpData().TimingStats().Task4().End(); 745 } 746 switch (quality_function_centrality_) { 747 case CEN_NONE: 748 //Nothing 749 break; 750 case CEN_LOG: 751 quality_function -= compl_inf*log(xi); 752 break; 753 case CEN_RECIPROCAL: 754 quality_function += compl_inf/xi; 755 case CEN_CUBED_RECIPROCAL: 756 quality_function += compl_inf/pow(xi,3); 757 break; 758 default: 759 DBG_ASSERT(false && "Unknown value for quality_function_centrality_"); 760 } 761 762 switch (quality_function_balancing_term_) { 763 case BT_NONE: 764 //Nothing 765 break; 766 case BT_CUBIC: 767 quality_function += pow(Max(0., Max(dual_inf,primal_inf)-compl_inf),3); 768 break; 769 default: 770 DBG_ASSERT(false && "Unknown value for quality_function_balancing term_"); 771 } 772 773 Jnlst().Printf(J_MOREDETAILED, J_BARRIER_UPDATE, 774 "sigma = %8.2e d_inf = %18.12e p_inf = %18.12e cmpl = %18.12e q = %18.12e a_pri = %8.2e a_dual = %8.2e xi = %8.2e\n", sigma, dual_inf, primal_inf, compl_inf, quality_function, alpha_primal, alpha_dual, xi); 775 776 777 return quality_function; 778 //return compl_inf; 779 } 780 781 Number PerformGoldenSection(Number sigma_up_in,Number q_up,Number sigma_lo_in,Number q_lo,Number sigma_tol,Number qf_tol,const Vector & step_aff_x_L,const Vector & step_aff_x_U,const Vector & step_aff_s_L,const Vector & step_aff_s_U,const Vector & step_aff_y_c,const Vector & step_aff_y_d,const Vector & step_aff_z_L,const Vector & step_aff_z_U,const Vector & step_aff_v_L,const Vector & step_aff_v_U,const Vector & step_cen_x_L,const Vector & step_cen_x_U,const Vector & step_cen_s_L,const Vector & step_cen_s_U,const Vector & step_cen_y_c,const Vector & step_cen_y_d,const Vector & step_cen_z_L,const Vector & step_cen_z_U,const Vector & step_cen_v_L,const Vector & step_cen_v_U)782 QualityFunctionMuOracle::PerformGoldenSection 783 (Number sigma_up_in, 784 Number q_up, 785 Number sigma_lo_in, 786 Number q_lo, 787 Number sigma_tol, 788 Number qf_tol, 789 const Vector& step_aff_x_L, 790 const Vector& step_aff_x_U, 791 const Vector& step_aff_s_L, 792 const Vector& step_aff_s_U, 793 const Vector& step_aff_y_c, 794 const Vector& step_aff_y_d, 795 const Vector& step_aff_z_L, 796 const Vector& step_aff_z_U, 797 const Vector& step_aff_v_L, 798 const Vector& step_aff_v_U, 799 const Vector& step_cen_x_L, 800 const Vector& step_cen_x_U, 801 const Vector& step_cen_s_L, 802 const Vector& step_cen_s_U, 803 const Vector& step_cen_y_c, 804 const Vector& step_cen_y_d, 805 const Vector& step_cen_z_L, 806 const Vector& step_cen_z_U, 807 const Vector& step_cen_v_L, 808 const Vector& step_cen_v_U 809 ) 810 { 811 Number sigma_up = ScaleSigma(sigma_up_in); 812 Number sigma_lo = ScaleSigma(sigma_lo_in); 813 814 Number sigma; 815 Number gfac = (3.-sqrt(5.))/2.; 816 Number sigma_mid1 = sigma_lo + gfac*(sigma_up-sigma_lo); 817 Number sigma_mid2 = sigma_lo + (1.-gfac)*(sigma_up-sigma_lo); 818 819 Number qmid1 = CalculateQualityFunction(UnscaleSigma(sigma_mid1), 820 step_aff_x_L, 821 step_aff_x_U, 822 step_aff_s_L, 823 step_aff_s_U, 824 step_aff_y_c, 825 step_aff_y_d, 826 step_aff_z_L, 827 step_aff_z_U, 828 step_aff_v_L, 829 step_aff_v_U, 830 step_cen_x_L, 831 step_cen_x_U, 832 step_cen_s_L, 833 step_cen_s_U, 834 step_cen_y_c, 835 step_cen_y_d, 836 step_cen_z_L, 837 step_cen_z_U, 838 step_cen_v_L, 839 step_cen_v_U); 840 Number qmid2 = CalculateQualityFunction(UnscaleSigma(sigma_mid2), 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 862 Index nsections = 0; 863 while ((sigma_up-sigma_lo)>=sigma_tol*sigma_up && 864 //while ((sigma_up-sigma_lo)>=sigma_tol && // Note we are using the non-relative criterion here for sigma 865 (1.-Min(q_lo, q_up, qmid1, qmid2)/Max(q_lo, q_up, qmid1, qmid2))>=qf_tol && 866 nsections<quality_function_max_section_steps_) { 867 nsections++; 868 //printf("sigma_lo=%e sigma_mid1=%e sigma_mid2=%e sigma_up=%e\n",sigma_lo, sigma_mid1, sigma_mid2, sigma_up); 869 if (qmid1 > qmid2) { 870 sigma_lo = sigma_mid1; 871 q_lo = qmid1; 872 sigma_mid1 = sigma_mid2; 873 qmid1 = qmid2; 874 sigma_mid2 = sigma_lo + (1.-gfac)*(sigma_up-sigma_lo); 875 qmid2 = CalculateQualityFunction(UnscaleSigma(sigma_mid2), 876 step_aff_x_L, 877 step_aff_x_U, 878 step_aff_s_L, 879 step_aff_s_U, 880 step_aff_y_c, 881 step_aff_y_d, 882 step_aff_z_L, 883 step_aff_z_U, 884 step_aff_v_L, 885 step_aff_v_U, 886 step_cen_x_L, 887 step_cen_x_U, 888 step_cen_s_L, 889 step_cen_s_U, 890 step_cen_y_c, 891 step_cen_y_d, 892 step_cen_z_L, 893 step_cen_z_U, 894 step_cen_v_L, 895 step_cen_v_U); 896 } 897 else { 898 sigma_up = sigma_mid2; 899 q_up = qmid2; 900 sigma_mid2 = sigma_mid1; 901 qmid2 = qmid1; 902 sigma_mid1 = sigma_lo + gfac*(sigma_up-sigma_lo); 903 qmid1 = CalculateQualityFunction(UnscaleSigma(sigma_mid1), 904 step_aff_x_L, 905 step_aff_x_U, 906 step_aff_s_L, 907 step_aff_s_U, 908 step_aff_y_c, 909 step_aff_y_d, 910 step_aff_z_L, 911 step_aff_z_U, 912 step_aff_v_L, 913 step_aff_v_U, 914 step_cen_x_L, 915 step_cen_x_U, 916 step_cen_s_L, 917 step_cen_s_U, 918 step_cen_y_c, 919 step_cen_y_d, 920 step_cen_z_L, 921 step_cen_z_U, 922 step_cen_v_L, 923 step_cen_v_U); 924 } 925 } 926 927 if ((sigma_up-sigma_lo)>=sigma_tol*sigma_up && 928 (1.-Min(q_lo, q_up, qmid1, qmid2)/Max(q_lo, q_up, qmid1, qmid2))<qf_tol) { 929 // The qf tolerance make it stop 930 IpData().Append_info_string("qf_tol "); 931 Number qf_min = Min(q_lo, q_up, qmid1, qmid2); 932 DBG_ASSERT(qf_min>-100.); 933 if (qf_min == q_lo) { 934 sigma = sigma_lo; 935 } 936 else if (qf_min == qmid1) { 937 sigma = sigma_mid1; 938 } 939 else if (qf_min == qmid2) { 940 sigma = sigma_mid2; 941 } 942 else { 943 sigma = sigma_up; 944 } 945 } 946 else { 947 Number q; 948 if (qmid1 < qmid2) { 949 sigma = sigma_mid1; 950 q = qmid1; 951 } 952 else { 953 sigma = sigma_mid2; 954 q = qmid2; 955 } 956 if (sigma_up == ScaleSigma(sigma_up_in)) { 957 Number qtmp; 958 if (q_up<0.) { 959 qtmp = CalculateQualityFunction(UnscaleSigma(sigma_up), 960 step_aff_x_L, 961 step_aff_x_U, 962 step_aff_s_L, 963 step_aff_s_U, 964 step_aff_y_c, 965 step_aff_y_d, 966 step_aff_z_L, 967 step_aff_z_U, 968 step_aff_v_L, 969 step_aff_v_U, 970 step_cen_x_L, 971 step_cen_x_U, 972 step_cen_s_L, 973 step_cen_s_U, 974 step_cen_y_c, 975 step_cen_y_d, 976 step_cen_z_L, 977 step_cen_z_U, 978 step_cen_v_L, 979 step_cen_v_U); 980 } 981 else { 982 qtmp = q_up; 983 } 984 if (qtmp < q) { 985 sigma = sigma_up; 986 q = qtmp; 987 } 988 } 989 else if (sigma_lo == ScaleSigma(sigma_lo_in)) { 990 Number qtmp; 991 if (q_lo<0.) { 992 qtmp = CalculateQualityFunction(UnscaleSigma(sigma_lo), 993 step_aff_x_L, 994 step_aff_x_U, 995 step_aff_s_L, 996 step_aff_s_U, 997 step_aff_y_c, 998 step_aff_y_d, 999 step_aff_z_L, 1000 step_aff_z_U, 1001 step_aff_v_L, 1002 step_aff_v_U, 1003 step_cen_x_L, 1004 step_cen_x_U, 1005 step_cen_s_L, 1006 step_cen_s_U, 1007 step_cen_y_c, 1008 step_cen_y_d, 1009 step_cen_z_L, 1010 step_cen_z_U, 1011 step_cen_v_L, 1012 step_cen_v_U); 1013 } 1014 else { 1015 qtmp = q_lo; 1016 } 1017 if (qtmp < q) { 1018 sigma = sigma_lo; 1019 q = qtmp; 1020 } 1021 } 1022 } 1023 1024 return UnscaleSigma(sigma); 1025 } 1026 1027 /* 1028 Number QualityFunctionMuOracle::ScaleSigma(Number sigma) {return log(sigma);} 1029 Number QualityFunctionMuOracle::UnscaleSigma(Number scaled_sigma) {return exp(scaled_sigma);} 1030 */ ScaleSigma(Number sigma)1031 Number QualityFunctionMuOracle::ScaleSigma(Number sigma) 1032 { 1033 return sigma; 1034 } UnscaleSigma(Number scaled_sigma)1035 Number QualityFunctionMuOracle::UnscaleSigma(Number scaled_sigma) 1036 { 1037 return scaled_sigma; 1038 } 1039 1040 /* AW: Tried search in the log space, but that was even worse than 1041 search in unscaled space */ 1042 /* 1043 Number 1044 QualityFunctionMuOracle::PerformGoldenSectionLog 1045 (Number sigma_up, 1046 Number sigma_lo, 1047 Number tol, 1048 const Vector& step_aff_x_L, 1049 const Vector& step_aff_x_U, 1050 const Vector& step_aff_s_L, 1051 const Vector& step_aff_s_U, 1052 const Vector& step_aff_y_c, 1053 const Vector& step_aff_y_d, 1054 const Vector& step_aff_z_L, 1055 const Vector& step_aff_z_U, 1056 const Vector& step_aff_v_L, 1057 const Vector& step_aff_v_U, 1058 const Vector& step_cen_x_L, 1059 const Vector& step_cen_x_U, 1060 const Vector& step_cen_s_L, 1061 const Vector& step_cen_s_U, 1062 const Vector& step_cen_y_c, 1063 const Vector& step_cen_y_d, 1064 const Vector& step_cen_z_L, 1065 const Vector& step_cen_z_U, 1066 const Vector& step_cen_v_L, 1067 const Vector& step_cen_v_U 1068 ) 1069 { 1070 Number log_sigma; 1071 Number log_sigma_up = log(sigma_up); 1072 Number log_sigma_lo = log(sigma_lo); 1073 1074 Number log_sigma_up_in = log_sigma_up; 1075 Number log_sigma_lo_in = log_sigma_lo; 1076 Number gfac = (3.-sqrt(5.))/2.; 1077 Number log_sigma_mid1 = log_sigma_lo + gfac*(log_sigma_up-log_sigma_lo); 1078 Number log_sigma_mid2 = log_sigma_lo + (1.-gfac)*(log_sigma_up-log_sigma_lo); 1079 1080 Number qmid1 = CalculateQualityFunction(exp(log_sigma_mid1), 1081 step_aff_x_L, 1082 step_aff_x_U, 1083 step_aff_s_L, 1084 step_aff_s_U, 1085 step_aff_y_c, 1086 step_aff_y_d, 1087 step_aff_z_L, 1088 step_aff_z_U, 1089 step_aff_v_L, 1090 step_aff_v_U, 1091 step_cen_x_L, 1092 step_cen_x_U, 1093 step_cen_s_L, 1094 step_cen_s_U, 1095 step_cen_y_c, 1096 step_cen_y_d, 1097 step_cen_z_L, 1098 step_cen_z_U, 1099 step_cen_v_L, 1100 step_cen_v_U); 1101 Number qmid2 = CalculateQualityFunction(exp(log_sigma_mid2), 1102 step_aff_x_L, 1103 step_aff_x_U, 1104 step_aff_s_L, 1105 step_aff_s_U, 1106 step_aff_y_c, 1107 step_aff_y_d, 1108 step_aff_z_L, 1109 step_aff_z_U, 1110 step_aff_v_L, 1111 step_aff_v_U, 1112 step_cen_x_L, 1113 step_cen_x_U, 1114 step_cen_s_L, 1115 step_cen_s_U, 1116 step_cen_y_c, 1117 step_cen_y_d, 1118 step_cen_z_L, 1119 step_cen_z_U, 1120 step_cen_v_L, 1121 step_cen_v_U); 1122 1123 Index nsections = 0; 1124 while ((log_sigma_up-log_sigma_lo)>=tol*log_sigma_up && nsections<quality_function_max_section_steps_) { 1125 nsections++; 1126 if (qmid1 > qmid2) { 1127 log_sigma_lo = log_sigma_mid1; 1128 log_sigma_mid1 = log_sigma_mid2; 1129 qmid1 = qmid2; 1130 log_sigma_mid2 = log_sigma_lo + (1.-gfac)*(log_sigma_up-log_sigma_lo); 1131 qmid2 = CalculateQualityFunction(exp(log_sigma_mid2), 1132 step_aff_x_L, 1133 step_aff_x_U, 1134 step_aff_s_L, 1135 step_aff_s_U, 1136 step_aff_y_c, 1137 step_aff_y_d, 1138 step_aff_z_L, 1139 step_aff_z_U, 1140 step_aff_v_L, 1141 step_aff_v_U, 1142 step_cen_x_L, 1143 step_cen_x_U, 1144 step_cen_s_L, 1145 step_cen_s_U, 1146 step_cen_y_c, 1147 step_cen_y_d, 1148 step_cen_z_L, 1149 step_cen_z_U, 1150 step_cen_v_L, 1151 step_cen_v_U); 1152 } 1153 else { 1154 log_sigma_up = log_sigma_mid2; 1155 log_sigma_mid2 = log_sigma_mid1; 1156 qmid2 = qmid1; 1157 log_sigma_mid1 = log_sigma_lo + gfac*(log_sigma_up-log_sigma_lo); 1158 qmid1 = CalculateQualityFunction(exp(log_sigma_mid1), 1159 step_aff_x_L, 1160 step_aff_x_U, 1161 step_aff_s_L, 1162 step_aff_s_U, 1163 step_aff_y_c, 1164 step_aff_y_d, 1165 step_aff_z_L, 1166 step_aff_z_U, 1167 step_aff_v_L, 1168 step_aff_v_U, 1169 step_cen_x_L, 1170 step_cen_x_U, 1171 step_cen_s_L, 1172 step_cen_s_U, 1173 step_cen_y_c, 1174 step_cen_y_d, 1175 step_cen_z_L, 1176 step_cen_z_U, 1177 step_cen_v_L, 1178 step_cen_v_U); 1179 } 1180 } 1181 1182 Number q; 1183 if (qmid1 < qmid2) { 1184 log_sigma = log_sigma_mid1; 1185 q = qmid1; 1186 } 1187 else { 1188 log_sigma = log_sigma_mid2; 1189 q = qmid2; 1190 } 1191 if (log_sigma_up == log_sigma_up_in) { 1192 Number qtmp = CalculateQualityFunction(exp(log_sigma_up), 1193 step_aff_x_L, 1194 step_aff_x_U, 1195 step_aff_s_L, 1196 step_aff_s_U, 1197 step_aff_y_c, 1198 step_aff_y_d, 1199 step_aff_z_L, 1200 step_aff_z_U, 1201 step_aff_v_L, 1202 step_aff_v_U, 1203 step_cen_x_L, 1204 step_cen_x_U, 1205 step_cen_s_L, 1206 step_cen_s_U, 1207 step_cen_y_c, 1208 step_cen_y_d, 1209 step_cen_z_L, 1210 step_cen_z_U, 1211 step_cen_v_L, 1212 step_cen_v_U); 1213 if (qtmp < q) { 1214 log_sigma = log_sigma_up; 1215 q = qtmp; 1216 } 1217 } 1218 else if (log_sigma_lo == log_sigma_lo_in) { 1219 Number qtmp = CalculateQualityFunction(exp(log_sigma_lo), 1220 step_aff_x_L, 1221 step_aff_x_U, 1222 step_aff_s_L, 1223 step_aff_s_U, 1224 step_aff_y_c, 1225 step_aff_y_d, 1226 step_aff_z_L, 1227 step_aff_z_U, 1228 step_aff_v_L, 1229 step_aff_v_U, 1230 step_cen_x_L, 1231 step_cen_x_U, 1232 step_cen_s_L, 1233 step_cen_s_U, 1234 step_cen_y_c, 1235 step_cen_y_d, 1236 step_cen_z_L, 1237 step_cen_z_U, 1238 step_cen_v_L, 1239 step_cen_v_U); 1240 if (qtmp < q) { 1241 log_sigma = log_sigma_lo; 1242 q = qtmp; 1243 } 1244 } 1245 1246 return exp(log_sigma); 1247 } 1248 */ 1249 1250 1251 } // namespace Ipopt 1252