MoochoPack : Framework for Large-Scale Optimization Algorithms Version of the Day
MoochoPack_TangentialStepWithInequStd_Step.cpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <math.h>
00030 
00031 #include <ostream>
00032 #include <sstream>
00033 
00034 #include "MoochoPack_TangentialStepWithInequStd_Step.hpp"
00035 #include "MoochoPack_moocho_algo_conversion.hpp"
00036 #include "MoochoPack_Exceptions.hpp"
00037 #include "IterationPack_print_algorithm_step.hpp"
00038 #include "ConstrainedOptPack_MatrixIdentConcat.hpp"
00039 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00040 #include "AbstractLinAlgPack_VectorOut.hpp"
00041 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00042 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00043 #include "Teuchos_dyn_cast.hpp"
00044 
00045 namespace {
00046 template< class T >
00047 inline
00048 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00049 template< class T >
00050 inline
00051 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
00052 } // end namespace
00053 
00054 namespace MoochoPack {
00055 
00056 TangentialStepWithInequStd_Step::TangentialStepWithInequStd_Step(
00057   const qp_solver_ptr_t       &qp_solver
00058   ,const qp_tester_ptr_t      &qp_tester
00059   ,value_type                 warm_start_frac
00060   ,EQPTesting                 qp_testing
00061   ,bool                       primal_feasible_point_error
00062   ,bool                       dual_feasible_point_error
00063   )
00064   :qp_solver_(qp_solver)
00065   ,qp_tester_(qp_tester)
00066   ,warm_start_frac_(warm_start_frac)
00067   ,qp_testing_(qp_testing)
00068   ,primal_feasible_point_error_(primal_feasible_point_error)
00069   ,dual_feasible_point_error_(dual_feasible_point_error)
00070   ,dl_iq_(dl_name)
00071   ,du_iq_(du_name)
00072 {}
00073 
00074 bool TangentialStepWithInequStd_Step::do_step(
00075   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00076   ,poss_type assoc_step_poss
00077   )
00078 {
00079 
00080   using Teuchos::RCP;
00081   using Teuchos::dyn_cast;
00082   using ::fabs;
00083   using LinAlgOpPack::Vt_S;
00084   using LinAlgOpPack::V_VpV;
00085   using LinAlgOpPack::V_VmV;
00086   using LinAlgOpPack::Vp_StV;
00087   using LinAlgOpPack::Vp_V;
00088   using LinAlgOpPack::V_StV;
00089   using LinAlgOpPack::V_MtV;
00090 //  using ConstrainedOptPack::min_abs;
00091   using AbstractLinAlgPack::max_near_feas_step;
00092   typedef VectorMutable::vec_mut_ptr_t   vec_mut_ptr_t;
00093 
00094   NLPAlgo &algo = rsqp_algo(_algo);
00095   NLPAlgoState &s = algo.rsqp_state();
00096   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00097   EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_journal_output_level();
00098   std::ostream &out = algo.track().journal_out();
00099   //const bool check_results = algo.algo_cntr().check_results();
00100 
00101   // print step header.
00102   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00103     using IterationPack::print_algorithm_step;
00104     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00105   }
00106 
00107   // problem dimensions
00108   const size_type
00109     //n  = algo.nlp().n(),
00110     m  = algo.nlp().m(),
00111     r  = s.equ_decomp().size();
00112 
00113   // Get the iteration quantity container objects
00114   IterQuantityAccess<value_type>
00115     &alpha_iq = s.alpha(),
00116     &zeta_iq  = s.zeta(),
00117     &eta_iq   = s.eta();
00118   IterQuantityAccess<VectorMutable>
00119     &dl_iq      = dl_iq_(s),
00120     &du_iq      = du_iq_(s),
00121     &nu_iq      = s.nu(),
00122     *c_iq       = m > 0  ? &s.c()       : NULL,
00123     *lambda_iq  = m > 0  ? &s.lambda()  : NULL,
00124     &rGf_iq     = s.rGf(),
00125     &w_iq       = s.w(),
00126     &qp_grad_iq = s.qp_grad(),
00127     &py_iq      = s.py(),
00128     &pz_iq      = s.pz(),
00129     &Ypy_iq     = s.Ypy(),
00130     &Zpz_iq     = s.Zpz();
00131   IterQuantityAccess<MatrixOp>
00132     &Z_iq   = s.Z(),
00133     //*Uz_iq  = (m > r)  ? &s.Uz() : NULL,
00134     *Uy_iq  = (m > r)  ? &s.Uy() : NULL;
00135   IterQuantityAccess<MatrixSymOp>
00136     &rHL_iq = s.rHL();
00137   IterQuantityAccess<ActSetStats>
00138     &act_set_stats_iq = act_set_stats_(s);
00139   
00140   // Accessed/modified/updated (just some)
00141   VectorMutable  *Ypy_k = (m ? &Ypy_iq.get_k(0) : NULL);
00142   const MatrixOp  &Z_k   = Z_iq.get_k(0);
00143   VectorMutable  &pz_k  = pz_iq.set_k(0);
00144   VectorMutable  &Zpz_k = Zpz_iq.set_k(0);
00145 
00146   // Comupte qp_grad which is an approximation to rGf + Z'*HL*Y*py
00147 
00148   // qp_grad = rGf
00149   VectorMutable
00150     &qp_grad_k = ( qp_grad_iq.set_k(0) = rGf_iq.get_k(0) );
00151 
00152   // qp_grad += zeta * w
00153   if( w_iq.updated_k(0) ) {
00154     if(zeta_iq.updated_k(0))
00155       Vp_StV( &qp_grad_k, zeta_iq.get_k(0), w_iq.get_k(0) );
00156     else
00157       Vp_V( &qp_grad_k, w_iq.get_k(0) );
00158   }
00159 
00160   //
00161   // Set the bounds for:
00162   //
00163   //   dl <= Z*pz + Y*py <= du  ->  dl - Ypy <= Z*pz <= du - Ypz
00164 
00165   vec_mut_ptr_t
00166     bl = s.space_x().create_member(),
00167     bu = s.space_x().create_member();
00168 
00169   if(m) {
00170     // bl = dl_k - Ypy_k
00171     V_VmV( bl.get(), dl_iq.get_k(0), *Ypy_k );
00172     // bu = du_k - Ypy_k
00173     V_VmV( bu.get(), du_iq.get_k(0), *Ypy_k );
00174   }
00175   else {
00176     *bl = dl_iq.get_k(0);
00177     *bu = du_iq.get_k(0);
00178   }
00179 
00180   // Print out the QP bounds for the constraints
00181   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00182     out << "\nqp_grad_k = \n" << qp_grad_k;
00183   }
00184   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00185     out << "\nbl = \n" << *bl;
00186     out << "\nbu = \n" << *bu;
00187   }
00188 
00189   //
00190   // Determine if we should perform a warm start or not.
00191   //
00192   bool do_warm_start = false;
00193   if( act_set_stats_iq.updated_k(-1) ) {
00194     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00195       out << "\nDetermining if the QP should use a warm start ...\n";
00196     }
00197     // We need to see if we should preform a warm start for the next iteration
00198     ActSetStats &stats = act_set_stats_iq.get_k(-1);
00199     const size_type
00200       num_active = stats.num_active(),
00201       num_adds   = stats.num_adds(),
00202       num_drops  = stats.num_drops();
00203     const value_type
00204       frac_same
00205       = ( num_adds == ActSetStats::NOT_KNOWN || num_active == 0
00206         ? 0.0
00207         : my_max(((double)(num_active)-num_adds-num_drops) / num_active, 0.0 ) );
00208     do_warm_start = ( num_active > 0 && frac_same >= warm_start_frac() );
00209     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00210       out << "\nnum_active = " << num_active;
00211       if( num_active ) {
00212         out << "\nmax(num_active-num_adds-num_drops,0)/(num_active) = "
00213           << "max("<<num_active<<"-"<<num_adds<<"-"<<num_drops<<",0)/("<<num_active<<") = "
00214           << frac_same;
00215         if( do_warm_start )
00216           out << " >= ";
00217         else
00218           out << " < ";
00219         out << "warm_start_frac = " << warm_start_frac();
00220       }
00221       if( do_warm_start )
00222         out << "\nUse a warm start this time!\n";
00223       else
00224         out << "\nDon't use a warm start this time!\n";
00225     }
00226   }
00227 
00228   // Use active set from last iteration as an estimate for current active set
00229   // if we are to use a warm start.
00230   // 
00231   // ToDo: If the selection of dependent and independent variables changes
00232   // then you will have to adjust this or not perform a warm start at all!
00233   if( do_warm_start ) {
00234     nu_iq.set_k(0,-1);
00235   }
00236   else {
00237     nu_iq.set_k(0) = 0.0; // No guess of the active set
00238   }
00239   VectorMutable  &nu_k  = nu_iq.get_k(0);
00240 
00241   //
00242   // Setup the reduced QP subproblem
00243   //
00244   // The call to the QP is setup for the more flexible call to the QPSolverRelaxed
00245   // interface to deal with the three independent variabilities: using simple
00246   // bounds for pz or not, general inequalities included or not, and extra equality
00247   // constraints included or not.
00248   // If this method of calling the QP solver were not used then 4 separate
00249   // calls to solve_qp(...) would have to be included to handle the four possible
00250   // QP formulations.
00251   //
00252 
00253   // The numeric arguments for the QP solver (in the nomenclatrue of QPSolverRelaxed)
00254 
00255   const value_type  qp_bnd_inf = NLP::infinite_bound();
00256 
00257   const Vector            &qp_g       = qp_grad_k;
00258   const MatrixSymOp       &qp_G       = rHL_iq.get_k(0);
00259   const value_type        qp_etaL     = 0.0;
00260   vec_mut_ptr_t           qp_dL       = Teuchos::null;
00261   vec_mut_ptr_t           qp_dU       = Teuchos::null;
00262   Teuchos::RCP<const MatrixOp>
00263                           qp_E        = Teuchos::null;
00264   BLAS_Cpp::Transp        qp_trans_E  = BLAS_Cpp::no_trans;
00265   vec_mut_ptr_t           qp_b        = Teuchos::null;
00266   vec_mut_ptr_t           qp_eL       = Teuchos::null;
00267   vec_mut_ptr_t           qp_eU       = Teuchos::null;
00268   Teuchos::RCP<const MatrixOp>
00269                           qp_F        = Teuchos::null;
00270   BLAS_Cpp::Transp        qp_trans_F  = BLAS_Cpp::no_trans;
00271   vec_mut_ptr_t           qp_f        = Teuchos::null;
00272   value_type              qp_eta      = 0.0;
00273   VectorMutable           &qp_d       = pz_k;  // pz_k will be updated directly!
00274   vec_mut_ptr_t           qp_nu       = Teuchos::null;
00275   vec_mut_ptr_t           qp_mu       = Teuchos::null;
00276   vec_mut_ptr_t           qp_Ed       = Teuchos::null;
00277   vec_mut_ptr_t           qp_lambda   = Teuchos::null;
00278 
00279   //
00280   // Determine if we can use simple bounds on pz.
00281   // 
00282   // If we have a variable-reduction null-space matrix
00283   // (with any choice for Y) then:
00284   // 
00285   // d = Z*pz + (1-eta) * Y*py
00286   // 
00287   // [ d(var_dep)   ]  = [ D ] * pz  + (1-eta) * [ Ypy(var_dep)   ]
00288   // [ d(var_indep) ]    [ I ]                   [ Ypy(var_indep) ]
00289   // 
00290   // For a cooridinate decomposition (Y = [ I ; 0 ]) then Ypy(var_indep) ==
00291   // 0.0 and in this case the bounds on d(var_indep) become simple bounds on
00292   // pz even with the relaxation.  Also, if dl(var_dep) and du(var_dep) are
00293   // unbounded, then we can also use simple bounds since we don't need the
00294   // relaxation and we can set eta=0. In this case we just have to subtract
00295   // from the upper and lower bounds on pz!
00296   // 
00297   // Otherwise, we can not use simple variable bounds and implement the
00298   // relaxation properly.
00299   // 
00300 
00301   const MatrixIdentConcat
00302     *Zvr = dynamic_cast<const MatrixIdentConcat*>( &Z_k );
00303   const Range1D
00304     var_dep   = Zvr ? Zvr->D_rng() : Range1D::Invalid,
00305     var_indep = Zvr ? Zvr->I_rng() : Range1D();
00306 
00307   RCP<Vector> Ypy_indep;
00308   const value_type
00309     Ypy_indep_norm_inf
00310     = ( m ? (Ypy_indep=Ypy_k->sub_view(var_indep))->norm_inf() : 0.0);
00311 
00312   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00313     out
00314       << "\nDetermine if we can use simple bounds on pz ...\n"
00315       << "    m = " << m << std::endl
00316       << "    dynamic_cast<const MatrixIdentConcat*>(&Z_k) = " << Zvr << std::endl
00317       << "    ||Ypy_k(var_indep)||inf = " << Ypy_indep_norm_inf << std::endl;
00318 
00319   const bool
00320     bounded_var_dep
00321     = (
00322       m > 0
00323       &&
00324       num_bounded( *bl->sub_view(var_dep), *bu->sub_view(var_dep), qp_bnd_inf )
00325       );
00326 
00327   const bool
00328     use_simple_pz_bounds
00329     = (
00330       m == 0
00331       ||
00332       ( Zvr != NULL && ( Ypy_indep_norm_inf == 0.0 || bounded_var_dep == 0 ) )
00333       );
00334 
00335   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00336     out
00337       << (use_simple_pz_bounds
00338           ? "\nUsing simple bounds on pz ...\n"
00339           : "\nUsing bounds on full Z*pz ...\n")
00340       << (bounded_var_dep
00341           ? "\nThere are finite bounds on dependent variables.  Adding extra inequality constrints for D*pz ...\n" 
00342           : "\nThere are no finite bounds on dependent variables.  There will be no extra inequality constraints added on D*pz ...\n" ) ;
00343 
00344   if( use_simple_pz_bounds ) {
00345     // Set simple bound constraints on pz
00346     qp_dL = bl->sub_view(var_indep);
00347     qp_dU = bu->sub_view(var_indep);
00348     qp_nu = nu_k.sub_view(var_indep); // nu_k(var_indep) will be updated directly!
00349     if( m && bounded_var_dep ) {
00350       // Set general inequality constraints for D*pz
00351       qp_E   = Teuchos::rcp(&Zvr->D(),false);
00352       qp_b   = Ypy_k->sub_view(var_dep);
00353       qp_eL  = bl->sub_view(var_dep);
00354       qp_eU  = bu->sub_view(var_dep);
00355       qp_mu  = nu_k.sub_view(var_dep);  // nu_k(var_dep) will be updated directly!
00356       qp_Ed  = Zpz_k.sub_view(var_dep); // Zpz_k(var_dep) will be updated directly!
00357     }
00358     else {
00359       // Leave these as NULL since there is no extra general inequality constraints
00360     }
00361   }
00362   else if( !use_simple_pz_bounds ) { // ToDo: Leave out parts for unbounded dependent variables!
00363     // There are no simple bounds! (leave qp_dL, qp_dU and qp_nu as null)
00364     // Set general inequality constraints for Z*pz
00365     qp_E   = Teuchos::rcp(&Z_k,false);
00366     qp_b   = Teuchos::rcp(Ypy_k,false);
00367     qp_eL  = bl;
00368     qp_eU  = bu;
00369     qp_mu  = Teuchos::rcp(&nu_k,false);
00370     qp_Ed  = Teuchos::rcp(&Zpz_k,false); // Zpz_k will be updated directly!
00371   }
00372   else {
00373     TEST_FOR_EXCEPT(true);
00374   }
00375 
00376   // Set the general equality constriants (if they exist)
00377   Range1D equ_undecomp = s.equ_undecomp();
00378   if( m > r && m > 0 ) {
00379     // qp_f = Uy_k * py_k + c_k(equ_undecomp)
00380     qp_f = s.space_c().sub_space(equ_undecomp)->create_member();
00381     V_MtV( qp_f.get(), Uy_iq->get_k(0), BLAS_Cpp::no_trans, py_iq.get_k(0) );
00382     Vp_V( qp_f.get(), *c_iq->get_k(0).sub_view(equ_undecomp) );
00383     // Must resize for the undecomposed constriants if it has not already been
00384     qp_F       = Teuchos::rcp(&Uy_iq->get_k(0),false);
00385     qp_lambda  = lambda_iq->set_k(0).sub_view(equ_undecomp); // lambda_k(equ_undecomp), will be updated directly!
00386   }
00387 
00388   // Setup the rest of the arguments
00389 
00390   QPSolverRelaxed::EOutputLevel
00391     qp_olevel;
00392   switch( olevel ) {
00393     case PRINT_NOTHING:
00394       qp_olevel = QPSolverRelaxed::PRINT_NONE;
00395       break;
00396     case PRINT_BASIC_ALGORITHM_INFO:
00397       qp_olevel = QPSolverRelaxed::PRINT_NONE;
00398       break;
00399     case PRINT_ALGORITHM_STEPS:
00400       qp_olevel = QPSolverRelaxed::PRINT_BASIC_INFO;
00401       break;
00402     case PRINT_ACTIVE_SET:
00403       qp_olevel = QPSolverRelaxed::PRINT_ITER_SUMMARY;
00404       break;
00405     case PRINT_VECTORS:
00406       qp_olevel = QPSolverRelaxed::PRINT_ITER_VECTORS;
00407       break;
00408     case PRINT_ITERATION_QUANTITIES:
00409       qp_olevel = QPSolverRelaxed::PRINT_EVERY_THING;
00410       break;
00411     default:
00412       TEST_FOR_EXCEPT(true);
00413   }
00414   // ToDo: Set print options so that only vectors matrices etc
00415   // are only printed in the null space
00416 
00417   //
00418   // Solve the QP
00419   // 
00420   qp_solver().infinite_bound(qp_bnd_inf);
00421   const QPSolverStats::ESolutionType
00422     solution_type =
00423     qp_solver().solve_qp(
00424       int(olevel) == int(PRINT_NOTHING) ? NULL : &out
00425       ,qp_olevel
00426       ,( algo.algo_cntr().check_results()
00427          ? QPSolverRelaxed::RUN_TESTS :  QPSolverRelaxed::NO_TESTS )
00428       ,qp_g, qp_G, qp_etaL, qp_dL.get(), qp_dU.get()
00429       ,qp_E.get(), qp_trans_E, qp_E.get() ? qp_b.get() : NULL
00430       ,qp_E.get() ? qp_eL.get() : NULL, qp_E.get() ? qp_eU.get() : NULL
00431       ,qp_F.get(), qp_trans_F, qp_F.get() ? qp_f.get() : NULL
00432       ,NULL // obj_d
00433       ,&qp_eta, &qp_d
00434       ,qp_nu.get()
00435       ,qp_mu.get(), qp_E.get() ? qp_Ed.get() : NULL
00436       ,qp_F.get() ? qp_lambda.get() : NULL
00437       ,NULL // qp_Fd
00438       );
00439 
00440   //
00441   // Check the optimality conditions for the QP
00442   //
00443   std::ostringstream omsg;
00444   bool throw_qp_failure = false;
00445   if(   qp_testing() == QP_TEST
00446     || ( qp_testing() == QP_TEST_DEFAULT && algo.algo_cntr().check_results() )  )
00447   {
00448     if( int(olevel) >= int(PRINT_ALGORITHM_STEPS) ) {
00449       out << "\nChecking the optimality conditions of the reduced QP subproblem ...\n";
00450     }
00451     if(!qp_tester().check_optimality_conditions(
00452       solution_type,qp_solver().infinite_bound()
00453       ,int(olevel) == int(PRINT_NOTHING) ? NULL : &out
00454       ,int(olevel) >= int(PRINT_VECTORS) ? true : false
00455       ,int(olevel) >= int(PRINT_ITERATION_QUANTITIES) ? true : false
00456       ,qp_g, qp_G, qp_etaL, qp_dL.get(), qp_dU.get()
00457       ,qp_E.get(), qp_trans_E, qp_E.get() ? qp_b.get() : NULL
00458       ,qp_E.get() ? qp_eL.get() : NULL, qp_E.get() ? qp_eU.get() : NULL
00459       ,qp_F.get(), qp_trans_F, qp_F.get() ? qp_f.get() : NULL
00460       ,NULL // obj_d
00461       ,&qp_eta, &qp_d
00462       ,qp_nu.get()
00463       ,qp_mu.get(), qp_E.get() ? qp_Ed.get() : NULL
00464       ,qp_F.get() ? qp_lambda.get() : NULL
00465       ,NULL // qp_Fd
00466       ))
00467     {
00468       omsg << "\n*** Alert! at least one of the QP optimality conditions did not check out.\n";
00469       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00470         out << omsg.str();
00471       }
00472       throw_qp_failure = true;
00473     }
00474   }
00475 
00476   //
00477   // Set the solution
00478   //
00479   if( !use_simple_pz_bounds ) {
00480     // Everything is already updated!
00481   }
00482   else if( use_simple_pz_bounds ) {
00483     // Just have to set Zpz_k(var_indep) = pz_k
00484     *Zpz_k.sub_view(var_indep) = pz_k;
00485     if( m && !bounded_var_dep ) {
00486       // Must compute Zpz_k(var_dep) = D*pz
00487       LinAlgOpPack::V_MtV( &*Zpz_k.sub_view(var_dep), Zvr->D(), BLAS_Cpp::no_trans, pz_k );
00488       // ToDo: Remove the compuation of Zpz here unless you must
00489     }
00490   }
00491   else {
00492     TEST_FOR_EXCEPT(true);
00493   }
00494 
00495   // Set the solution statistics
00496   qp_solver_stats_(s).set_k(0) = qp_solver().get_qp_stats();
00497 
00498   // Cut back Ypy_k = (1-eta) * Ypy_k
00499   const value_type eps = std::numeric_limits<value_type>::epsilon();
00500   if( fabs(qp_eta - 0.0) > eps ) {
00501     if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00502       out
00503         << "\n*** Alert! the QP was infeasible (eta = "<<qp_eta<<").  Cutting back Ypy_k = (1.0 - eta)*Ypy_k  ...\n";
00504     }
00505     Vt_S( Ypy_k, 1.0 - qp_eta );
00506   }
00507 
00508   // eta_k
00509   eta_iq.set_k(0) = qp_eta;
00510 
00511   //
00512   // Modify the solution if we have to!
00513   // 
00514   switch(solution_type) {
00515     case QPSolverStats::OPTIMAL_SOLUTION:
00516       break;  // we are good!
00517     case QPSolverStats::PRIMAL_FEASIBLE_POINT:
00518     {
00519       omsg
00520         << "\n*** Alert! the returned QP solution is PRIMAL_FEASIBLE_POINT but not optimal!\n";
00521       if( primal_feasible_point_error() )
00522         omsg
00523           << "\n*** primal_feasible_point_error == true, this is an error!\n";
00524       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00525         out << omsg.str();
00526       }
00527       throw_qp_failure = primal_feasible_point_error();
00528       break;
00529     } 
00530     case QPSolverStats::DUAL_FEASIBLE_POINT:
00531     {
00532       omsg
00533         << "\n*** Alert! the returned QP solution is DUAL_FEASIBLE_POINT"
00534         << "\n*** but not optimal so we cut back the step ...\n";
00535       if( dual_feasible_point_error() )
00536         omsg
00537           << "\n*** dual_feasible_point_error == true, this is an error!\n";
00538       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00539         out << omsg.str();
00540       }
00541       // Cut back the step to fit in the bounds
00542       // 
00543       // dl <= u*(Ypy_k+Zpz_k) <= du
00544       //
00545       vec_mut_ptr_t
00546         zero  = s.space_x().create_member(0.0),
00547         d_tmp = s.space_x().create_member();
00548       V_VpV( d_tmp.get(), *Ypy_k, Zpz_k );
00549       const std::pair<value_type,value_type>
00550         u_steps = max_near_feas_step( *zero, *d_tmp, dl_iq.get_k(0), du_iq.get_k(0), 0.0 );
00551       const value_type
00552         u = my_min( u_steps.first, 1.0 ); // largest positive step size
00553       alpha_iq.set_k(0) = u;
00554       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00555         out << "\nFinding u s.t. dl <= u*(Ypy_k+Zpz_k) <= du\n"
00556           << "max step length u = " << u << std::endl
00557           << "alpha_k = u = " << alpha_iq.get_k(0) << std::endl;
00558       }
00559       throw_qp_failure = dual_feasible_point_error();
00560       break;
00561     } 
00562     case QPSolverStats::SUBOPTIMAL_POINT:
00563     {
00564       omsg
00565         << "\n*** Alert!, the returned QP solution is SUBOPTIMAL_POINT!\n";
00566       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00567         out << omsg.str();
00568       }
00569       throw_qp_failure = true;
00570       break;
00571     }
00572     default:
00573       TEST_FOR_EXCEPT(true);  // should not happen!
00574   }
00575 
00576   //
00577   // Output the final solution!
00578   //
00579   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00580     out << "\n||pz_k||inf    = " << s.pz().get_k(0).norm_inf()
00581       << "\nnu_k.nz()      = " << s.nu().get_k(0).nz()
00582       << "\nmax(|nu_k(i)|) = " << s.nu().get_k(0).norm_inf()
00583 //      << "\nmin(|nu_k(i)|) = " << min_abs( s.nu().get_k(0)() )
00584       ;
00585     if( m > r ) out << "\n||lambda_k(undecomp)||inf = " << s.lambda().get_k(0).norm_inf();
00586     out << "\n||Zpz_k||2     = " << s.Zpz().get_k(0).norm_2()
00587       ;
00588     if(qp_eta > 0.0) out << "\n||Ypy||2 = " << s.Ypy().get_k(0).norm_2();
00589     out << std::endl;
00590   }
00591 
00592   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00593     out << "\npz_k = \n" << s.pz().get_k(0);
00594     if(var_indep.size())
00595       out << "\nnu_k(var_indep) = \n" << *s.nu().get_k(0).sub_view(var_indep);
00596   }
00597 
00598   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00599     if(var_indep.size())
00600       out << "\nZpz(var_indep)_k = \n" << *s.Zpz().get_k(0).sub_view(var_indep);
00601     out << std::endl;
00602   }
00603 
00604   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00605     if(var_dep.size())
00606       out << "\nZpz(var_dep)_k = \n" << *s.Zpz().get_k(0).sub_view(var_dep);
00607     out << "\nZpz_k = \n" << s.Zpz().get_k(0);
00608     out << std::endl;
00609   }
00610 
00611   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00612     out << "\nnu_k = \n" << s.nu().get_k(0);
00613     if(var_dep.size())
00614       out << "\nnu_k(var_dep) = \n" << *s.nu().get_k(0).sub_view(var_dep);
00615     if( m > r )
00616       out << "\nlambda_k(equ_undecomp) = \n" << *s.lambda().get_k(0).sub_view(equ_undecomp);
00617     if(qp_eta > 0.0) out << "\nYpy = \n" << s.Ypy().get_k(0);
00618   }
00619 
00620   if( qp_eta == 1.0 ) {
00621     omsg
00622       << "TangentialStepWithInequStd_Step::do_step(...) : Error, a QP relaxation parameter\n"
00623       << "of eta = " << qp_eta << " was calculated and therefore it must be assumed\n"
00624       << "that the NLP's constraints are infeasible\n"
00625       << "Throwing an InfeasibleConstraints exception!\n";
00626     if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00627       out << omsg.str();
00628     }
00629     throw InfeasibleConstraints(omsg.str());
00630   }
00631 
00632   if( throw_qp_failure )
00633     throw QPFailure( omsg.str(), qp_solver().get_qp_stats() );
00634 
00635   return true;
00636 }
00637 
00638 void TangentialStepWithInequStd_Step::print_step(
00639   const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type
00640   ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
00641   ) const
00642 {
00643   out
00644     << L << "*** Calculate the null-space step by solving a constrained QP\n"
00645     << L << "qp_grad_k = rGf_k\n"
00646     << L << "if w_k is updated then set qp_grad_k = qp_grad_k + zeta_k * w_k\n"
00647     << L << "bl = dl_k - Ypy_k\n"
00648     << L << "bu = du_k - Ypy_k\n"
00649     << L << "etaL = 0.0\n"
00650     << L << "*** Determine if we can use simple bounds on pz or not\n"
00651     << L << "if num_bounded(bl_k(var_dep),bu_k(var_dep)) > 0 then\n"
00652     << L << "  bounded_var_dep = true\n"
00653     << L << "else\n"
00654     << L << "  bounded_var_dep = false\n"
00655     << L << "end\n"
00656     << L << "if( m==0\n"
00657     << L << "    or\n"
00658     << L << "    ( Z_k is a variable reduction null space matrix\n"
00659     << L << "      and\n"
00660     << L << "      ( ||Ypy_k(var_indep)||inf == 0 or bounded_var_dep==false ) )\n"
00661     << L << "  ) then\n"
00662     << L << "  use_simple_pz_bounds = true\n"
00663     << L << "else\n"
00664     << L << "  use_simple_pz_bounds = false\n"
00665     << L << "end\n"
00666     << L << "*** Setup QP arguments\n"
00667     << L << "qp_g = qp_grad_k\n"
00668     << L << "qp_G = rHL_k\n"
00669     << L << "if (use_simple_pz_bounds == true) then\n"
00670     << L << "  qp_dL = bl(var_indep), qp_dU = bu(var_indep))\n"
00671     << L << "  if (m > 0) then\n"
00672     << L << "    qp_E  = Z_k.D,       qp_b  = Ypy_k(var_dep)\n"
00673     << L << "    qp_eL = bl(var_dep), qp_eU = bu(var_dep)\n"
00674     << L << "  end\n"
00675     << L << "elseif (use_simple_pz_bounds == false) then\n"
00676     << L << "  qp_dL = -inf,  qp_dU = +inf\n"
00677     << L << "  qp_E  = Z_k,   qp_b  = Ypy_k\n"
00678     << L << "  qp_eL = bl,    qp_eU = bu\n"
00679     << L << "end\n"
00680     << L << "if (m > r) then\n"
00681     << L << "  qp_F  = V_k,     qp_f  = Uy_k*py_k + c_k(equ_undecomp)\n"
00682     << L << "else\n"
00683     << L << "  qp_F  = empty,   qp_f  = empty\n"
00684     << L << "end\n"
00685     << L << "Given active-set statistics (act_set_stats_km1)\n"
00686     << L << "  frac_same = max(num_active-num_adds-num_drops,0)/(num_active)\n"
00687     << L << "Use a warm start when frac_same >= warm_start_frac\n"
00688     << L << "Solve the following QP to compute qp_d, qp_eta, qp_Ed = qp_E * qp_d\n"
00689     << L << ",qp_nu, qp_mu and qp_lambda (" << typeName(qp_solver()) << "):\n"
00690     << L << "  min  qp_g' * qp_d + 1/2 * qp_d' * qp_G * qp_d + M(eta)\n"
00691     << L << "  qp_d <: R^(n-r)\n"
00692     << L << "  s.t.\n"
00693     << L << "       etaL  <=  eta\n"
00694     << L << "       qp_dL <= qp_d                         <= qp_dU   [qp_nu]\n"
00695     << L << "       qp_eL <= qp_E * qp_d + (1-eta)*qp_b   <= qp_eU   [qp_mu]\n"
00696     << L << "                qp_F * d_qp + (1-eta) * qp_f  = 0       [qp_lambda]\n"
00697     << L << "if (qp_testing==QP_TEST) or (fd_deriv_testing==QP_TEST_DEFAULT\n"
00698     << L << "and check_results==true) then\n"
00699     << L << "  Check the optimality conditions of the above QP\n"
00700     << L << "  if the optimality conditions do not check out then\n"
00701     << L << "    set throw_qp_failure = true\n"
00702     << L << "  end\n"
00703     << L << "end\n"
00704     << L << "*** Set the solution to the QP subproblem\n"
00705     << L << "pz_k  = qp_d\n"
00706     << L << "eta_k = qp_eta\n"
00707     << L << "if (use_simple_pz_bounds == true) then\n"
00708     << L << "  nu_k(var_dep)  = qp_mu,  nu_k(var_indep)  = qp_nu\n"
00709     << L << "  Zpz_k(var_dep) = qp_Ed,  Zpz_k(var_indep) = pz_k\n"
00710     << L << "elseif (use_simple_pz_bounds == false) then\n"
00711     << L << "  nu_k  = qp_mu\n"
00712     << L << "  Zpz_k = qp_Ed\n"
00713     << L << "end\n"
00714     << L << "if m > r then\n"
00715     << L << "  lambda_k(equ_undecomp) = qp_lambda\n"
00716     << L << "end\n"
00717     << L << "if (eta_k > 0) then set Ypy_k = (1-eta_k) * Ypy_k\n"
00718     << L << "if QP solution is suboptimal then\n"
00719     << L << "  throw_qp_failure = true\n"
00720     << L << "elseif QP solution is primal feasible (not optimal) then\n"
00721     << L << "  throw_qp_failure = primal_feasible_point_error\n"
00722     << L << "elseif QP solution is dual feasible (not optimal) then\n"
00723     << L << "  find max u s.t.\n"
00724     << L << "    dl_k <= u*(Ypy_k+Zpz_k) <= du_k\n"
00725     << L << "  alpha_k = u\n"
00726     << L << "  throw_qp_failure = dual_feasible_point_error\n"
00727     << L << "end\n"
00728     << L << "if (eta_k == 1.0) then\n"
00729     << L << "  The constraints are infeasible!\n"
00730     << L << "  throw InfeasibleConstraints(...)\n"
00731     << L << "end\n"
00732     << L << "if (throw_qp_failure == true) then\n"
00733     << L << "  throw QPFailure(...)\n"
00734     << L << "end\n"
00735     ;
00736 }
00737 
00738 } // end namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends