MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MoochoPack_LineSearch2ndOrderCorrect_Step.cpp
Go to the documentation of this file.
1 #if 0
2 
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
7 // Copyright (2003) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
40 //
41 // ***********************************************************************
42 // @HEADER
43 
44 #include <ostream>
45 #include <typeinfo>
46 
55 #include "ConstrainedOptPack/src/VectorWithNorms.h"
56 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixOp.hpp"
57 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
58 #include "AbstractLinAlgPack/src/max_near_feas_step.hpp"
63 
64 namespace LinAlgOpPack {
66 }
67 
68 namespace MoochoPack {
69 
71  const direct_ls_sqp_ptr_t& direct_ls_sqp
72  ,const merit_func_ptr_t& merit_func
73  ,const feasibility_step_ptr_t& feasibility_step
74  ,const direct_ls_newton_ptr_t& direct_ls_newton
75  ,value_type eta
76  ,ENewtonOutputLevel newton_olevel
77  ,value_type constr_norm_threshold
78  ,value_type constr_incr_ratio
79  ,int after_k_iter
80  ,EForcedConstrReduction forced_constr_reduction
81  ,value_type forced_reduct_ratio
82  ,value_type max_step_ratio
83  ,int max_newton_iter
84  )
85  :direct_ls_sqp_(direct_ls_sqp)
86  ,merit_func_(merit_func)
87  ,feasibility_step_(feasibility_step)
88  ,direct_ls_newton_(direct_ls_newton)
89  ,eta_(eta)
90  ,newton_olevel_(newton_olevel)
91  ,constr_norm_threshold_(constr_norm_threshold)
92  ,constr_incr_ratio_(constr_incr_ratio)
93  ,after_k_iter_(after_k_iter)
94  ,forced_constr_reduction_(forced_constr_reduction)
95  ,forced_reduct_ratio_(forced_reduct_ratio)
96  ,max_step_ratio_(max_step_ratio)
97  ,max_newton_iter_(max_newton_iter)
98 {}
99 
101  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
102  ,poss_type assoc_step_poss
103  )
104 {
105 
106 /* ToDo: Upate the below code!
107 
108  using std::setw;
109 
110  using DenseLinAlgPack::dot;
111  using DenseLinAlgPack::norm_inf;
112  using DenseLinAlgPack::V_VpV;
113  using DenseLinAlgPack::V_VmV;
114  using DenseLinAlgPack::Vp_StV;
115  using DenseLinAlgPack::Vt_S;
116 
117  using LinAlgOpPack::Vp_V;
118  using LinAlgOpPack::V_MtV;
119 
120  using AbstractLinAlgPack::max_near_feas_step;
121 
122  using ConstrainedOptPack::print_vector_change_stats;
123 
124  typedef LineSearch2ndOrderCorrect_Step this_t;
125 
126  NLPAlgo &algo = rsqp_algo(_algo);
127  NLPAlgoState &s = algo.rsqp_state();
128  NLP &nlp = algo.nlp();
129 
130  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
131  std::ostream& out = algo.track().journal_out();
132  out << std::boolalpha;
133 
134  // print step header.
135  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
136  using IterationPack::print_algorithm_step;
137  print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
138  }
139 
140  // /////////////////////////////////////////
141  // Set references to iteration quantities
142  //
143  // Set k+1 first then go back to get k to ensure
144  // we have backward storage.
145 
146  DVector
147  &x_kp1 = s.x().set_k(+1).v();
148  value_type
149  &f_kp1 = s.f().set_k(+1);
150  DVector
151  &c_kp1 = s.c().set_k(+1).v();
152 
153  const value_type
154  &f_k = s.f().get_k(0);
155  const DVector
156  &c_k = s.c().get_k(0).v();
157  const DVector
158  &x_k = s.x().get_k(0).v();
159  const DVector
160  &d_k = s.d().get_k(0).v();
161  value_type
162  &alpha_k = s.alpha().get_k(0);
163 
164  // /////////////////////////////////////
165  // Compute Dphi_k, phi_kp1 and phi_k
166 
167  // Dphi_k
168  const value_type
169  Dphi_k = merit_func().deriv();
170  if( Dphi_k >= 0 ) {
171  throw LineSearchFailure( "LineSearch2ndOrderCorrect_Step::do_step(...) : "
172  "Error, d_k is not a descent direction for the merit function " );
173  }
174 
175  // ph_kp1
176  value_type
177  &phi_kp1 = s.phi().set_k(+1) = merit_func().value( f_kp1, c_kp1 );
178 
179  // Must compute phi(x) at the base point x_k since the penalty parameter may have changed.
180  const value_type
181  &phi_k = s.phi().set_k(0) = merit_func().value( f_k, c_k );
182 
183  // //////////////////////////////////////
184  // Setup the calculation merit function
185 
186  // Here f_kp1, and c_kp1 are updated at the same time the
187  // line search is being performed.
188  nlp.set_f( &f_kp1 );
189  nlp.set_c( &c_kp1 );
190  MeritFuncCalcNLP
191  phi_calc( &merit_func(), &nlp );
192 
193  // //////////////////////////////////////////////////
194  // Concider 2nd order correction if near solution?
195 
196  bool considering_correction = false;
197  {
198  const value_type
199  small_num = std::numeric_limits<value_type>::min(),
200  nrm_c_k = s.c().get_k(0).norm_inf(),
201  nrm_c_kp1 = s.c().get_k(+1).norm_inf();
202  const bool
203  test_1 = nrm_c_k <= constr_norm_threshold(),
204  test_2 = (nrm_c_kp1/(1.0 + nrm_c_k)) < constr_incr_ratio(),
205  test_3 = s.k() >= after_k_iter();
206  considering_correction = test_1 && test_2 && test_3;
207  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
208  out << "\n||c_k||inf = " << nrm_c_k << (test_1 ? " <= " : " > " )
209  << "constr_norm_threshold = " << constr_norm_threshold()
210  << "\n||c_kp1||inf/(1.0+||c_k||inf) = "
211  << "(" << nrm_c_kp1 << ")/(" << 1.0 << " + " << nrm_c_k << ") = "
212  << ( nrm_c_kp1 / (1.0 + nrm_c_k ) ) << (test_2 ? " <= " : " > " )
213  << "constr_incr_ratio = " << constr_incr_ratio()
214  << "\nk = " << s.k() << (test_3 ? " >= " : " < ")
215  << "after_k_iter = " << after_k_iter()
216  << (considering_correction
217  ? ("\nThe computation of a 2nd order correction for x_kp1 = x_k + alpha_k*d_k + alpha_k^2*w"
218  " will be considered ...\n")
219  : "\nThe critera for considering a 2nd order correction has not been met ...\n" );
220  }
221  }
222 
223  // //////////////////////////////
224  // See if we can take a full step
225 
226  bool chose_point = false;
227 
228  const value_type frac_phi = phi_k + eta() * Dphi_k;
229  const bool armijo_test = phi_kp1 <= frac_phi;
230  if( armijo_test ) {
231  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
232  out << "\nAccepting full step x_kp1 = x_k + d_k\n";
233  }
234  chose_point = true; // The point meets the Armijo test.
235  }
236 
237  // This is storage for function and gradient evaluations for
238  // the trial newton points and must be remembered for latter
239  value_type f_xdww;
240  DVector c_xdww;
241  DVector w(x_kp1.size()), // Full correction after completed computation.
242  xdww(x_kp1.size()); // Will be set to xdw + sum( w(newton_i), newton_i = 1... )
243  // where w(itr) is the local corrections for the current
244  // newton iteration.
245  bool use_correction = false;
246  bool newton_failed = true;
247 
248  bool considered_correction = ( considering_correction && !chose_point );
249  if( considered_correction ) {
250 
251  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
252  out << "\nConsidering whether to compute a 2nd order correction for\n"
253  << "x_kp1 = x_k + alpha_k * d_k + alpha_k^2 * w ...\n";
254  }
255 
256  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
257  const value_type obj_descent = dot( s.Gf().get_k(0)(), d_k() );
258  out << "\nGf_k' * d_k = " << obj_descent << std::endl;
259  if( obj_descent >= 0.0 ) {
260  out << "\nWarning, this may not work well with Gf_k'*d_k >= 0.0\n";
261  }
262  }
263 
264  // Merit function for newton line searches
265  ConstrainedOptPack::MeritFuncNLESqrResid
266  phi_c;
267 
268  DVector
269  xdw = x_kp1; // Will be set to x + d + sum(w(i),i=1..itr-1)
270  // where w(i) are previous local corrections
271  value_type
272  phi_c_x = phi_c.value( c_k() ),
273  phi_c_xd = phi_c.value( c_kp1() ),
274  phi_c_xdw = phi_c_xd, // No correction is computed yet so w = 0
275  phi_c_xdww = phi_c_xdw,
276  nrm_d = norm_inf( d_k() );
277 
278  // Merit function for newton line searches
279  nlp.set_f( &(f_xdww = f_kp1) );
280  nlp.set_c( &(c_xdww = c_kp1) );
281  ConstrainedOptPack::MeritFuncCalcNLE
282  phi_c_calc( &phi_c, &nlp );
283 
284  DVector wy(s.con_decomp().size()); // Range space wy (see latter).
285 
286  const bool sufficient_reduction =
287  phi_c_xd < forced_reduct_ratio() * phi_c_x;
288  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
289  out << "\nphi_c(c(x_k+d_k)) = " << phi_c_xd << (sufficient_reduction ? " <= " : " > ")
290  << "forced_reduct_ratio* phi_c(c(x_k)) = " << forced_reduct_ratio() << " * " << phi_c_x
291  << " = " << (forced_reduct_ratio()*phi_c_x)
292  << (sufficient_reduction
293  ? "\nNo need for a 2nd order correciton, perform regular line search ... \n"
294  : "\nWe need to try to compute a correction w ...\n" );
295  }
296  if(sufficient_reduction) {
297  use_correction = false;
298  }
299  else {
300  // Try to compute a second order correction term.
301 
302  // Set print level newton iterations
303  ENewtonOutputLevel use_newton_olevel;
304  if( newton_olevel() == PRINT_USE_DEFAULT ) {
305  switch(olevel) {
306  case PRINT_NOTHING:
307  case PRINT_BASIC_ALGORITHM_INFO:
308  use_newton_olevel = PRINT_NEWTON_NOTHING;
309  break;
310  case PRINT_ALGORITHM_STEPS:
311  case PRINT_ACTIVE_SET:
312  use_newton_olevel = PRINT_NEWTON_SUMMARY_INFO;
313  break;
314  case PRINT_VECTORS:
315  use_newton_olevel = PRINT_NEWTON_STEPS;
316  break;
317  case PRINT_ITERATION_QUANTITIES:
318  use_newton_olevel = PRINT_NEWTON_VECTORS;
319  break;
320  default:
321  TEUCHOS_TEST_FOR_EXCEPT(true);
322  }
323  }
324  else {
325  use_newton_olevel = newton_olevel();
326  }
327  EJournalOutputLevel inner_olevel;
328  switch(use_newton_olevel) {
329  case PRINT_NEWTON_NOTHING:
330  case PRINT_NEWTON_SUMMARY_INFO:
331  inner_olevel = PRINT_NOTHING;
332  break;
333  case PRINT_NEWTON_STEPS:
334  inner_olevel = PRINT_ALGORITHM_STEPS;
335  break;
336  case PRINT_NEWTON_VECTORS:
337  if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES )
338  inner_olevel = PRINT_ITERATION_QUANTITIES;
339  else if( (int)olevel >= (int)PRINT_ACTIVE_SET )
340  inner_olevel = PRINT_ACTIVE_SET;
341  else
342  inner_olevel = PRINT_VECTORS;
343  break;
344  default:
345  TEUCHOS_TEST_FOR_EXCEPT(true);
346  }
347 
348  // Print header for summary information
349  const int dbl_min_w = 21;
350  const int dbl_w = std::_MAX(dbl_min_w,int(out.precision()+8));
351  if( use_newton_olevel == PRINT_NEWTON_SUMMARY_INFO ) {
352  out << "\nStarting Newton iterations ...\n"
353  << "\nphi_c_x = " << phi_c_x
354  << "\nphi_c_xd = " << phi_c_xd
355  << "\n||d_k||nf = " << nrm_d << "\n\n"
356  << setw(5) << "it"
357  << setw(dbl_w) << "||w||inf"
358  << setw(dbl_w) << "u"
359  << setw(dbl_w) << "step_ratio"
360  << setw(5) << "lsit"
361  << setw(dbl_w) << "a"
362  << setw(dbl_w) << "phi_c_xdww"
363  << setw(dbl_w) << "phi_c_xdww-phi_c_x"
364  << setw(dbl_w) << "phi_c_xdww-phi_c_xd\n"
365  << setw(5) << "----"
366  << setw(dbl_w) << "--------------------"
367  << setw(dbl_w) << "-------------------"
368  << setw(dbl_w) << "-------------------"
369  << setw(5) << "----"
370  << setw(dbl_w) << "-------------------"
371  << setw(dbl_w) << "-------------------"
372  << setw(dbl_w) << "-------------------"
373  << setw(dbl_w) << "-------------------\n";
374  }
375 
376  // Perform newton feasibility iterations
377  int newton_i;
378  for( newton_i = 1; newton_i <= max_newton_iter(); ++newton_i ) {
379 
380  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS ) {
381  out << "\n**** newton_i = " << newton_i << std::endl;
382  }
383 
384  // Compute a feasibility step
385  if(!feasibility_step().compute_feasibility_step(
386  out,inner_olevel,&algo,&s,xdw,nlp.c()(),&w() ))
387  {
388  if( (int)use_newton_olevel == (int)PRINT_NEWTON_SUMMARY_INFO ) {
389  out << "\nCould not compute feasible direction!\n";
390  }
391  break; // exit the newton iterations
392  }
393 
394  value_type
395  nrm_w = norm_inf(w());
396 
397  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS ) {
398  out << "\n||w||inf = " << nrm_w << std::endl;
399  }
400 
401  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_VECTORS ) {
402  out << "\nw = " << w();
403  }
404 
405  // ////////////////////////////////
406  // Cutting back w
407 
408  value_type a = 1.0; // This is the alpha for your linesearch
409 
410  // Cut back w to be in the relaxed bounds.
411  std::pair<value_type,value_type>
412  u_steps = max_near_feas_step( s.x().get_k(0)(), w()
413  , algo.nlp().xl(), algo.nlp().xu()
414  , algo.algo_cntr().max_var_bounds_viol() );
415  const value_type u = u_steps.first;
416 
417  if( u < a ) {
418  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS ) {
419  out << "\nCutting back w = (a=u) * w to be within relaxed bounds:\n"
420  << "u = " << u << std::endl;
421  }
422  a = u;
423  }
424 
425  // Cut back step so x+d+sum(w(i)) is not too far from x+d
426  value_type
427  step_ratio = nrm_w / nrm_d;
428  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS ) {
429  out << "\nstep_ratio = ||w||inf/||d||inf = " << step_ratio
430  << std::endl;
431  }
432  if( a * step_ratio > max_step_ratio() ) {
433  const value_type aa = a*(max_step_ratio()/step_ratio);
434  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS ) {
435  out << "\na*step_ratio = " << a*step_ratio
436  << " > max_step_ratio = " << max_step_ratio() << std::endl
437  << "Cutting back a = a*max_step_ratio/step_ratio = "
438  << aa << std::endl;
439  }
440  a = aa;
441  }
442 
443  // /////////////////////////////////////////////////
444  // Perform a line search along xdww = xdw + a * w
445 
446  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS ) {
447  out << "\nPerform linesearch along xdww = xdw + a*w\n"
448  << "starting from a = " << a << " ...\n";
449  }
450 
451  xdww = xdw(); // xdww = xdw + a*w
452  Vp_StV( &xdww(), a, w() );
453  phi_c.calc_deriv(nlp.c()); // Set the directional derivative at c(xdw)
454  phi_c_xdww = phi_c_calc( xdww() ); // phi_c_xdww = phi(xdww)
455  const DVectorSlice xdw_w[2] = { xdw(), w() };
456  MeritFuncCalc1DQuadratic
457  phi_c_calc_1d( phi_c_calc, 1 , xdw_w, &xdww() );
458  bool ls_okay = false;
459  try {
460  ls_okay = direct_ls_newton().do_line_search(
461  phi_c_calc_1d,phi_c_xdw
462  ,&a,&phi_c_xdww
463  , (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_STEPS
464  ? &out : 0
465  );
466  }
467  catch(const DirectLineSearch_Strategy::NotDescentDirection& excpt ) {
468  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO ) {
469  out << "\nThe line search object throw the exception:" << typeName(excpt) << ":\n"
470  << excpt.what() << std::endl;
471  }
472  ls_okay = false;
473  }
474  // Note that the last value c(x) computed by the line search is for
475  // xdw + a*w.
476 
477  // Print line for summary output
478  if( use_newton_olevel == PRINT_NEWTON_SUMMARY_INFO ) {
479  out << setw(5) << newton_i
480  << setw(dbl_w) << nrm_w
481  << setw(dbl_w) << u
482  << setw(dbl_w) << step_ratio
483  << setw(5) << direct_ls_newton().num_iterations()
484  << setw(dbl_w) << a
485  << setw(dbl_w) << phi_c_xdww
486  << setw(dbl_w) << (phi_c_xdww-phi_c_x)
487  << setw(dbl_w) << (phi_c_xdww-phi_c_xd) << std::endl;
488  }
489 
490  if(!ls_okay) {
491  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO ) {
492  out << "\nThe line search failed so forget about computing a correction ...\n";
493  }
494  use_correction = false;
495  newton_failed = true;
496  break;
497  }
498 
499  // See if this point is okay
500  bool good_correction = false;
501  switch( forced_constr_reduction() ) {
502  case CONSTR_LESS_X_D: {
503  good_correction = ( phi_c_xdww < forced_reduct_ratio()*phi_c_xd );
504  if( good_correction
505  && (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO )
506  {
507  out << "\nphi_c(c(x_k+d_k+w)) = " << phi_c_xdww
508  << " < forced_reduct_ratio * phi_c(c(x_k+d_k)) = "
509  << forced_reduct_ratio() << " * " << phi_c_xd
510  << " = " << (forced_reduct_ratio()*phi_c_xd) << std::endl;
511  }
512  break;
513  }
514  case CONSTR_LESS_X: {
515  good_correction = ( phi_c_xdww < forced_reduct_ratio()*phi_c_x );
516  if( good_correction
517  && (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO )
518  {
519  out << "\nphi_c(c(x_k+d_k+w)) = " << phi_c_xdww
520  << " < forced_reduct_ratio * phi_c(c(x_k)) = "
521  << forced_reduct_ratio() << " * " << phi_c_x
522  << " = " << (forced_reduct_ratio()*phi_c_x) << std::endl;
523  }
524  break;
525  }
526  default:
527  TEUCHOS_TEST_FOR_EXCEPT(true);
528  }
529 
530  if(good_correction) {
531  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO ) {
532  out << "\nAccept this point and compute our full correction w ... \n";
533  }
534  // Compute the full correction and do a curved linesearch
535  // w = xdww - x_kp1
536  V_VmV( &w(), xdww(), x_kp1() );
537  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO ) {
538  out << "\n||w||inf = " << norm_inf(w()) << std::endl;
539  }
540  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_VECTORS ) {
541  out << "\nw = " << w();
542  }
543  use_correction = true;
544  newton_failed = false;
545  break;
546  }
547 
548  // Else perform another newton iteration.
549  xdw = xdww;
550  phi_c_xdw = phi_c_xdww;
551 
552  } // end for
553  if( !use_correction ) {
554  newton_failed = true;
555  if( forced_constr_reduction() == CONSTR_LESS_X_D ) {
556  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
557  out << "\nDam! This is really bad!\n"
558  << "We where only looking for point phi_c(c(x_k+d_k+w)"
559  << " < phi_c(c(x_k+k_k) and we could not find it\n"
560  << " in the aloted number of newton iterations!\n"
561  << "Perhaps the Gc_k did not give us a descent direction?\n"
562  << "Just perform a standard line search from here ...\n";
563  }
564  use_correction = false;
565  }
566  else {
567  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
568  out << "\nWe where looking for point phi_c(c(x_k+d_k+w))"
569  << " < phi_c(c(x_k)) and we could not find it.\n";
570  }
571  if( phi_c_xdww < phi_c_xd ) {
572  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
573  out << "However, we did find a point less than phi_c(c(x_k+d_k))\n"
574  << "so lets use the correction anyway.\n";
575  }
576  // Compute the full correction and do a curved linesearch
577  // w = xdww - x_kp1
578  V_VmV( &w(), xdww(), x_kp1() );
579  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_SUMMARY_INFO ) {
580  out << "\n||w||inf = " << norm_inf(w()) << std::endl;
581  }
582  if( (int)use_newton_olevel >= (int)this_t::PRINT_NEWTON_VECTORS ) {
583  out << "\nw = " << w();
584  }
585  use_correction = true;
586  }
587  else {
588  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
589  out << "Dam! We did not even find a point less than phi_c(c(x_k+d_k))\n"
590  << "just perform a standard line search along x_k + alpha_k*d_k.\n";
591  }
592  use_correction = false;
593  }
594  }
595  }
596  } // end else from if phi_c_xdw > phi_c_x
597  } // end considered_correction
598 
599  // //////////////////////////
600  // Set up for the line search
601 
602  if( considered_correction ) {
603  if( use_correction ) {
604  // We are using the correction so setup the full step for the
605  // NLP linesearch to come.
606  Vp_V( &x_kp1(), w() ); // Set point to x_kp1 = x_k + d_k + w
607  nlp.calc_f(x_kp1(),false); // same as last call to calc_c(x)
608  f_kp1 = nlp.f(); // Here f and c where computed at x_k+d_k+w
609  c_kp1 = nlp.c()();
610  phi_kp1 = merit_func().value( f_kp1, c_kp1 );
611  }
612  else {
613  // Just pretend the second order correction never happened
614  // and we don't need to do anything.
615  }
616  // Set back the base point
617  nlp.set_f( &f_kp1 );
618  nlp.set_c( &c_kp1 );
619  }
620 
621  // //////////////////////
622  // Do the line search
623 
624  if( !chose_point ) {
625  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
626  if( use_correction ) {
627  out << "\nPerform a curved linesearch along:\n"
628  << "x_kp1 = x_k + alpha_k * d_k + alpha_k^2 * w ...\n";
629  }
630  else {
631  out << "\nPerform a standard linesearch along:\n"
632  << "x_kp1 = x_k + alpha_k * d_k ...\n";
633  }
634  }
635  const DVectorSlice xdw[3] = { x_k(), d_k(), w() };
636  MeritFuncCalc1DQuadratic
637  phi_calc_1d( phi_calc, (use_correction?2:1) , xdw, &x_kp1() );
638  if( !direct_ls_sqp().do_line_search( phi_calc_1d, phi_k, &alpha_k, &phi_kp1
639  , static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ?
640  &out : static_cast<std::ostream*>(0) ) )
641  {
642  // If the line search failed but the value of the merit function is less than
643  // the base point value then just accept it and move on. This many be a
644  // bad thing to do.
645 
646  const value_type
647  scaled_ared = (s.phi().get_k(0) - s.phi().get_k(+1))/s.phi().get_k(0),
648  keep_on_frac = 1.0e-10; // Make adjustable?
649  bool keep_on = scaled_ared < keep_on_frac;
650 
651  if( (int)olevel >= (int)PRINT_BASIC_ALGORITHM_INFO )
652  {
653  out
654  << "\nThe maximum number of linesearch iterations has been exceeded "
655  << "(k = " << algo.state().k() << ")\n"
656  << "(phi_k - phi_kp1)/phi_k = " << scaled_ared;
657 // if(keep_on) {
658 // out
659 // << " < " << keep_on_frac
660 // << "\nso we will accept to step and move on.\n";
661 // }
662 // else {
663 // out
664 // << " > " << keep_on_frac
665 // << "\nso we will reject the step and declare a line search failure.\n";
666 // }
667  }
668 //
669 // if( keep_on ) return true;
670 
671  throw LineSearchFailure( "LineSearch2ndOrderCorrect_Step::do_step(): "
672  "Error, Line search failure" );
673  }
674  }
675 
676  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
677  out << "\nalpha_k = " << alpha_k << std::endl
678  << "\n||x_kp1||inf = " << norm_inf( x_kp1 ) << std::endl
679  << "\nf_kp1 = " << f_kp1 << std::endl
680  << "\n||c_kp1||inf = " << norm_inf(c_kp1) << std::endl
681  << "\nphi_kp1 = " << phi_kp1 << std::endl;
682  }
683 
684  if( (int)olevel >= (int)PRINT_VECTORS ) {
685  out << "\nx_kp1 =\n" << x_kp1
686  << "\nc_kp1 =\n" << c_kp1;
687  }
688 
689 */
691 
692  return true;
693 }
694 
696  const Algorithm& algo
697  , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
698  , std::ostream& out, const std::string& L ) const
699 {
700  out << L << "*** Calculate a second order correction when near solution.\n"
701  << L << "*** If we can compute a correction w then perform a curved\n"
702  << L << "*** line search along x_kp1 = x_k + alpha_k * d_k + alpha_k^2 * w.\n"
703  << L << "default: eta = " << eta() << std::endl
704  << L << " constr_norm_threshold = " << constr_norm_threshold() << std::endl
705  << L << " constr_incr_ratio = " << constr_incr_ratio() << std::endl
706  << L << " after_k_iter = " << after_k_iter() << std::endl
707  << L << " forced_constr_reduction = " << (forced_constr_reduction()== CONSTR_LESS_X_D
708  ? "CONSTR_LESS_X_D\n"
709  : "CONSTR_LESS_X\n" )
710  << L << " forced_reduct_ratio = " << forced_reduct_ratio() << std::endl
711  << L << " max_step_ratio = " << max_step_ratio() << std::endl
712  << L << " max_newton_iter = " << max_newton_iter() << std::endl
713  << L << "begin definition of NLP merit function phi.value(f(x),c(x)):\n";
714 
715  merit_func().print_merit_func( out, L + " " );
716 
717  out << L << "end definition\n"
718  << L << "Dphi_k = phi.deriv()\n"
719  << L << "if Dphi_k >= 0 then\n"
720  << L << " throw line_search_failure\n"
721  << L << "end\n"
722  << L << "phi_kp1 = phi_k.value(f_kp1,c_kp1)\n"
723  << L << "phi_k = phi.value(f_k,c_k)\n"
724  << L << "if (norm(c_k,inf) <= constr_norm_threshold)\n"
725  << L << "and (norm(c_kp1,inf)/(norm(c_k,inf)+small_num) <= constr_incr_ratio)\n"
726  << L << "and (k >= after_k_iter) then\n"
727  << L << "considering_correction = ( (norm(c_k,inf) <= constr_norm_threshold)\n"
728  << L << " and (norm(c_kp1,inf)/(1.0 + norm(c_k,inf)) <= constr_incr_ratio)\n"
729  << L << " and (k >= after_k_iter) )\n"
730  << L << "chose_point = false\n"
731  << L << "if phi_kp1 < phi_k + eta * Dphi_k then\n"
732  << L << " chose_point = true\n"
733  << L << "else\n"
734  << L << "if (considering_correction == true) and (chose_point == false) then\n"
735  << L << " considered_correction = true\n"
736  << L << " begin definition of c(x) merit function phi_c.value(c(x)):\n";
737 
739  out, L + " " );
740 
741  out << L << " end definition\n"
742  << L << " xdw = x_kp1;\n"
743  << L << " phi_c_x = phi_c.value(c_k);\n"
744  << L << " phi_c_xd = phi_c.value(c_kp1);\n"
745  << L << " phi_c_xdw = phi_c_xd;\n"
746  << L << " phi_c_xdww = phi_c_xdw;\n"
747  << L << " if phi_c_xd < forced_reduct_ratio * phi_c_x then\n"
748  << L << " *** There is no need to perform a correction.\n"
749  << L << " use_correction = false;\n"
750  << L << " else\n"
751  << L << " *** Lets try to compute a correction by performing\n"
752  << L << " *** a series of newton steps to compute local steps w\n"
753  << L << " for newton_i = 1...max_newton_itr\n"
754  << L << " begin feasibility step calculation: \"" << typeName(feasibility_step()) << "\"\n";
755 
756  feasibility_step().print_step( out, L + " " );
757 
758  out << L << " end feasibility step calculation\n"
759  << L << " Find the largest positive step u where the slightly\n"
760  << L << " relaxed variable bounds:\n"
761  << L << " xl - delta <= xdw + u * w <= xu + delta\n"
762  << L << " are strictly satisfied (where delta = max_var_bounds_viol).\n"
763  << L << " a = min(1,u);\n"
764  << L << " step_ratio = norm(w,inf)/norm(d,inf);\n"
765  << L << " a = min(a,max_step_ratio/step_ratio);\n"
766  << L << " Perform line search for phi_c.value(c(xdww = xdw+a*w))\n"
767  << L << " starting from a and compute:\n"
768  << L << " a,\n"
769  << L << " xdww = xdw + a * w,\n"
770  << L << " phi_c_xdww = phi_c.value(c(xdww))\n"
771  << L << " print summary information;\n"
772  << L << " if line search failed then\n"
773  << L << " use_correction = false;\n"
774  << L << " exit for loop;\n"
775  << L << " end\n"
776  << L << " *** Determine if this is sufficent reduction in c(x) error\n"
777  << L << " if forced_constr_reduction == CONSTR_LESS_X_D then\n"
778  << L << " good_correction = (phi_c.value(c(xdww)) < forced_reduct_ratio*phi_c_xd);\n"
779  << L << " else if forced_constr_reduction == CONSTR_LESS_X then\n"
780  << L << " good_correction = (phi_c.value(c(xdww)) < forced_reduct_ratio*phi_c_x);\n"
781  << L << " end\n"
782  << L << " if good_correction == true then\n"
783  << L << " w = xdww - (x_k+d_k);\n"
784  << L << " use_correction = true;\n"
785  << L << " exit for loop;\n"
786  << L << " end\n"
787  << L << " *** This is not sufficient reduction is c(x) error yet.\n"
788  << L << " xdw = xdww;\n"
789  << L << " phi_c_xdw = phi_c_xdww;\n"
790  << L << " end\n"
791  << L << " if use_correction == false then\n"
792  << L << " if forced_constr_reduction == CONSTR_LESS_X_D then\n"
793  << L << " *** Dam! We could not find a point phi_c_xdww < phi_c_xd.\n"
794  << L << " *** Perhaps Gc_k does not give a descent direction for phi_c!\n"
795  << L << " else if forced_constr_reduction == CONSTR_LESS_X then\n"
796  << L << " if phi_c_dww < phi_c_xd then\n"
797  << L << " *** Accept this correction anyway.\n"
798  << L << " use_correction = true\n"
799  << L << " else\n"
800  << L << " *** Dam! we could not find any reduction in phi_c so\n"
801  << L << " *** Perhaps Gc_k does not give a descent direction for phi_c!\n"
802  << L << " end\n"
803  << L << " end\n"
804  << L << " end\n"
805  << L << "end\n"
806  << L << "if chose_point == false then\n"
807  << L << " if use_correction == true then\n"
808  << L << " Perform line search along x_kp1 = x_k + alpha_k * d_k + alpha_k^2 * w\n"
809  << L << " else\n"
810  << L << " Perform line search along x_kp1 = x_k + alpha_k * d_k\n"
811  << L << " end\n"
812  << L << " begin direct line search : \"" << typeName(direct_ls_sqp()) << "\"\n";
813 
814  direct_ls_sqp().print_algorithm( out, L + " " );
815 
816  out
817  << L << " end direct line search\n"
818  << L << " if maximum number of linesearch iterations are exceeded then\n"
819  << L << " throw line_search_failure\n"
820  << L << " end\n"
821  << L << "end\n";
822 }
823 
824 } // end namespace MoochoPack
825 
826 #endif // 0
void print_merit_func(std::ostream &out, const std::string &leading_str) const
std::string typeName(const T &t)
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
std::ostream * out
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
v_lhs = alpha * op(M_rhs1) * v_rhs2 + beta * v_lhs (BLAS xGEMV)
A merit function for the square of the constriant values.
AbstractLinAlgPack::value_type value_type
void print_step(const Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss, std::ostream &out, const std::string &leading_str) const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
LineSearch2ndOrderCorrect_Step(const direct_ls_sqp_ptr_t &direct_ls_sqp=NULL, const merit_func_ptr_t &merit_func=NULL, const feasibility_step_ptr_t &feasibility_step=NULL, const direct_ls_newton_ptr_t &direct_ls_newton=0, value_type eta=1.0e-4, ENewtonOutputLevel newton_olevel=PRINT_USE_DEFAULT, value_type constr_norm_threshold=1.0, value_type constr_incr_ratio=10.0, int after_k_iter=0, EForcedConstrReduction forced_constr_reduction=CONSTR_LESS_X, value_type forced_reduct_ratio=1.0, value_type max_step_ratio=1.0, int max_newton_iter=3)