MoochoPack_EvalNewPointTailoredApproach_Step.cpp

Go to the documentation of this file.
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 <ostream>
00030 
00031 #include "MoochoPack_EvalNewPointTailoredApproach_Step.hpp"
00032 #include "MoochoPack_Exceptions.hpp"
00033 #include "MoochoPack_moocho_algo_conversion.hpp"
00034 #include "IterationPack_print_algorithm_step.hpp"
00035 #include "ConstrainedOptPack_MatrixIdentConcatStd.hpp"
00036 #include "NLPInterfacePack_NLPDirect.hpp"
00037 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00038 #include "AbstractLinAlgPack_VectorMutable.hpp"
00039 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00040 #include "AbstractLinAlgPack_VectorOut.hpp"
00041 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00042 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00043 #include "Teuchos_dyn_cast.hpp"
00044 #include "Teuchos_TestForException.hpp"
00045 
00046 namespace MoochoPack {
00047 
00048 EvalNewPointTailoredApproach_Step::EvalNewPointTailoredApproach_Step(
00049   const deriv_tester_ptr_t      &deriv_tester
00050   ,const bounds_tester_ptr_t    &bounds_tester
00051   , EFDDerivTesting             fd_deriv_testing
00052   )
00053   :deriv_tester_(deriv_tester)
00054   ,bounds_tester_(bounds_tester)
00055   ,fd_deriv_testing_(fd_deriv_testing)
00056 {}
00057 
00058 bool EvalNewPointTailoredApproach_Step::do_step(
00059   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00060   ,poss_type assoc_step_poss
00061   )
00062 {
00063   using Teuchos::rcp;
00064   using Teuchos::dyn_cast;
00065   using AbstractLinAlgPack::assert_print_nan_inf;
00066   using LinAlgOpPack::V_MtV;
00067   using IterationPack::print_algorithm_step;
00068 
00069   NLPAlgo             &algo   = rsqp_algo(_algo);
00070   NLPAlgoState        &s      = algo.rsqp_state();
00071   NLPDirect           &nlp    = dyn_cast<NLPDirect>(algo.nlp());
00072 
00073   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00074   EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_journal_output_level();
00075   std::ostream& out = algo.track().journal_out();
00076 
00077   // print step header.
00078   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00079     using IterationPack::print_algorithm_step;
00080     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00081   }
00082 
00083   if(!nlp.is_initialized())
00084     nlp.initialize(algo.algo_cntr().check_results());
00085 
00086   Teuchos::VerboseObjectTempState<NLP>
00087     nlpOutputTempState(rcp(&nlp,false),Teuchos::getFancyOStream(rcp(&out,false)),convertToVerbLevel(olevel));
00088 
00089   const Range1D
00090     var_dep = nlp.var_dep(),
00091     var_indep = nlp.var_indep();
00092 
00093   const size_type
00094     n  = nlp.n(),
00095     m  = nlp.m(),
00096     r  = var_dep.size();
00097 
00098   TEST_FOR_EXCEPTION(
00099     m > r, TestFailed
00100     ,"EvalNewPointTailoredApproach_Step::do_step(...) : Error, "
00101     "Undecomposed equalities are supported yet!" );
00102 
00103   IterQuantityAccess<VectorMutable>
00104     &x_iq = s.x();
00105 
00106   if( x_iq.last_updated() == IterQuantity::NONE_UPDATED ) {
00107     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00108       out << "\nx is not updated for any k so set x_k = nlp.xinit() ...\n";
00109     }
00110     x_iq.set_k(0) = nlp.xinit();
00111   }
00112 
00113   // Validate x
00114   if(algo.algo_cntr().check_results()) {
00115     assert_print_nan_inf(
00116       x_iq.get_k(0), "x_k", true
00117       , int(olevel) >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL
00118       );
00119     if( nlp.num_bounded_x() > 0 ) {
00120       if(!bounds_tester().check_in_bounds(
00121            int(olevel)  >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL
00122            ,int(olevel) >= int(PRINT_VECTORS)                // print_all_warnings
00123            ,int(olevel) >= int(PRINT_ITERATION_QUANTITIES)  // print_vectors
00124            ,nlp.xl(),        "xl"
00125            ,nlp.xu(),        "xu"
00126            ,x_iq.get_k(0),   "x_k"
00127            ))
00128       {
00129         TEST_FOR_EXCEPTION(
00130           true, TestFailed
00131           ,"EvalNewPointTailoredApproach_Step::do_step(...) : Error, "
00132           "the variables bounds xl <= x_k <= xu where violated!" );
00133       }
00134     }
00135   }
00136 
00137   Vector &x = x_iq.get_k(0);
00138 
00139   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00140     out << "\n||x_k||inf            = " << x.norm_inf();
00141     if( var_dep.size() )
00142       out << "\n||x(var_dep)_k||inf   = " << x.sub_view(var_dep)->norm_inf();
00143     if( var_indep.size() )
00144       out << "\n||x(var_indep)_k||inf = " << x.sub_view(var_indep)->norm_inf();
00145     out << std::endl;
00146   }
00147   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00148     out << "\nx_k = \n" << x;
00149     if( var_dep.size() )
00150       out << "\nx(var_dep)_k = \n" << *x.sub_view(var_dep);
00151   }
00152   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00153     if( var_indep.size() )
00154       out << "\nx(var_indep)_k = \n" << *x.sub_view(var_indep);
00155   }
00156 
00157   // If c_k is not updated then we must compute it
00158   bool recalc_c = true;
00159 
00160   if( !s.c().updated_k(0) ) {
00161     s.c().set_k(0);
00162     recalc_c = true;
00163   }
00164   else {
00165     recalc_c = false;
00166   }
00167     
00168   // Get references to Z, Y, Uz and Uy
00169   MatrixOp
00170     &Z_k  = s.Z().set_k(0),
00171     &Y_k  = s.Y().set_k(0),
00172     *Uz_k = (m > r) ? &s.Uz().set_k(0) : NULL,
00173     *Uy_k = (m > r) ? &s.Uy().set_k(0) : NULL;
00174   MatrixIdentConcatStd
00175     &cZ_k = dyn_cast<MatrixIdentConcatStd>(Z_k);
00176   // Release any references to D in Y or Uy
00177    uninitialize_Y_Uy(&Y_k,Uy_k);
00178   // If Z has not been initialized or Z.D is being shared by someone else we need to reconstruct Z.D
00179   bool reconstruct_Z_D = (cZ_k.rows() == n || cZ_k.cols() != n-r || cZ_k.D_ptr().count() > 1);
00180   MatrixIdentConcatStd::D_ptr_t
00181     D_ptr = Teuchos::null;
00182   if( reconstruct_Z_D )
00183     D_ptr = nlp.factory_D()->create();
00184   else
00185     D_ptr = cZ_k.D_ptr();
00186 
00187   // Compute all the quantities.
00188   const bool supports_Gf = nlp.supports_Gf();
00189   Teuchos::RefCountPtr<MatrixOp>
00190     GcU = (m > r) ? nlp.factory_GcU()->create() : Teuchos::null; // ToDo: Reuse GcU somehow? 
00191   VectorMutable
00192     &py_k  = s.py().set_k(0);
00193   nlp.unset_quantities();
00194   nlp.calc_point(
00195     x                                                     // x
00196     ,!s.f().updated_k(0) ? &s.f().set_k(0) : NULL         // f
00197     ,&s.c().get_k(0)                                      // c
00198     ,recalc_c                                             // recalc_c
00199     ,supports_Gf?&s.Gf().set_k(0):NULL                    // Gf
00200     ,&py_k                                                // -inv(C)*c
00201     ,&s.rGf().set_k(0)                                    // rGf
00202     ,GcU.get()                                            // GcU
00203     ,const_cast<MatrixOp*>(D_ptr.get())                   // -inv(C)*N
00204     ,Uz_k                                                 // Uz
00205     );
00206   s.equ_decomp(   nlp.con_decomp()   );
00207   s.equ_undecomp( nlp.con_undecomp() );
00208 
00209   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00210     out << "\nQuantities computed directly from NLPDirect nlp object ...\n";
00211     out << "\nf_k           = " << s.f().get_k(0);
00212     out << "\n||c_k||inf    = " << s.c().get_k(0).norm_inf();
00213     if(supports_Gf)
00214       out << "\n||Gf_k||inf   = " << s.Gf().get_k(0).norm_inf();
00215     out << "\n||py_k||inf   = " << s.py().get_k(0).norm_inf();
00216     out << "\n||rGf_k||inf  = " << s.rGf().get_k(0).norm_inf();
00217     out << std::endl;
00218   }
00219 
00220   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00221     out << "\nc_k  = \n" << s.c().get_k(0);
00222     if(supports_Gf)
00223       out << "\nGf_k = \n" << s.Gf().get_k(0);
00224     out << "\npy_k = \n" << s.py().get_k(0);
00225     out << "\nrGf_k = \n" << s.rGf().get_k(0);
00226     out << std::endl;
00227   }
00228   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00229     out << "\nD = -inv(C)*N = \n" << *D_ptr;
00230     out << std::endl;
00231   }
00232 
00233   if(algo.algo_cntr().check_results()) {
00234     assert_print_nan_inf(s.f().get_k(0),   "f_k",true,&out); 
00235     assert_print_nan_inf(s.c().get_k(0),   "c_k",true,&out); 
00236     if(supports_Gf)
00237       assert_print_nan_inf(s.Gf().get_k(0),  "Gf_k",true,&out); 
00238     assert_print_nan_inf(s.py().get_k(0),  "py_k",true,&out);
00239     assert_print_nan_inf(s.rGf().get_k(0), "rGf_k",true,&out);
00240   }
00241 
00242   // Check the derivatives if we are checking the results
00243   if(   fd_deriv_testing() == FD_TEST
00244     || ( fd_deriv_testing() == FD_DEFAULT && algo.algo_cntr().check_results() )  )
00245   {
00246     
00247     if( olevel >= PRINT_ALGORITHM_STEPS ) {
00248       out << "\n*** Checking derivatives by finite differences\n";
00249     }
00250     
00251     const bool has_bounds = nlp.num_bounded_x() > 0;
00252     const bool nlp_passed = deriv_tester().finite_diff_check(
00253       &nlp
00254       ,x
00255       ,has_bounds ? &nlp.xl() : (const Vector*)NULL
00256       ,has_bounds ? &nlp.xu() : (const Vector*)NULL
00257       ,&s.c().get_k(0)
00258       ,supports_Gf?&s.Gf().get_k(0):NULL
00259       ,&s.py().get_k(0)
00260       ,&s.rGf().get_k(0)
00261       ,GcU.get()
00262       ,D_ptr.get()
00263       ,Uz_k
00264       ,olevel >= PRINT_VECTORS
00265       ,( olevel >= PRINT_ALGORITHM_STEPS ) ? &out : (std::ostream*)NULL
00266       );
00267     TEST_FOR_EXCEPTION(
00268       !nlp_passed, TestFailed
00269       ,"EvalNewPointTailoredApproach_Step::do_step(...) : Error, "
00270       "the tests of the nlp derivatives failed!" );
00271   }
00272 
00273   if( reconstruct_Z_D ) {
00274     //
00275     // Z = [      D     ] space_xD
00276     //     [      I     ] space_xI
00277     //        space_xI
00278     //
00279     cZ_k.initialize(
00280       nlp.space_x()                                          // space_cols
00281       ,nlp.space_x()->sub_space(nlp.var_indep())->clone()    // space_rows
00282       ,MatrixIdentConcatStd::TOP                             // top_or_bottom
00283       ,1.0                                                   // alpha
00284       ,D_ptr                                                 // D_ptr
00285       ,BLAS_Cpp::no_trans                                    // D_trans
00286       );
00287   }
00288 
00289   // Compute py, Y and Uy
00290   calc_py_Y_Uy( nlp, D_ptr, &py_k, &Y_k, Uy_k, olevel, out ); 
00291 
00292   // Compute Ypy = Y*py
00293   VectorMutable
00294     &Ypy_k  = s.Ypy().set_k(0);
00295   V_MtV( &Ypy_k, Y_k, BLAS_Cpp::no_trans, py_k );
00296 
00297   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00298     out << "\nQuantities after computing final py and Ypy ...\n";
00299     out << "\n||py_k||inf   = " << s.py().get_k(0).norm_inf();
00300     out << "\n||Ypy_k||inf  = " << s.Ypy().get_k(0).norm_inf();
00301     out << std::endl;
00302   }
00303 
00304   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00305     out << "\npy_k = \n"  << s.py().get_k(0);
00306     out << "\nYpy_k = \n" << s.Ypy().get_k(0);
00307     out << std::endl;
00308   }
00309 
00310   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00311     out << "\nZ_k = \n" << s.Z().get_k(0);
00312     out << "\nY_k = \n" << s.Y().get_k(0);
00313     out << std::endl;
00314   }
00315 
00316   return true;
00317 
00318 }
00319 
00320 void EvalNewPointTailoredApproach_Step::print_step(
00321   const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type
00322   ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
00323   ) const
00324 {
00325   out
00326     << L << "*** Evaluate the new point for the \"Tailored Approach\"\n"
00327     << L << "if nlp is not initialized then initialize the nlp\n"
00328     << L << "if x is not updated for any k then set x_k = xinit\n"
00329     << L << "if Gf is supported then set Gf_k = Gf(x_k) <: space_x\n"
00330     << L << "For Gc(:,equ_decomp) = [ C' ; N' ] <: space_x:space_c(equ_decomp) compute:\n"
00331     << L << "  py_k = -inv(C)*c_k\n"
00332     << L << "  D = -inv(C)*N <: R^(n x (n-m))\n"
00333     << L << "  rGf_k = Gf_k(var_indep) + D'*Gf_k(var_dep)\n"
00334     << L << "  Z_k = [ D ; I ] <: R^(n x (n-m))\n"
00335     << L << "  if m > r Uz_k = Gc(var_indep,equ_undecomp)' + Gc(var_dep,equ_undecomp)'*D\n"
00336     << L << "if ( (fd_deriv_testing==FD_TEST)\n"
00337     << L << "    or (fd_deriv_testing==FD_DEFAULT and check_results==true\n"
00338     << L << "  ) then\n"
00339     << L << "  check Gf_k, py_k, rGf_k, D, Uz (if m > r) and Vz (if mI > 0) by finite differences.\n"
00340     << L << "end\n";
00341   print_calc_py_Y_Uy( out, L );
00342   out
00343     << L << "Ypy_k = Y_k * py_k\n"
00344     << L << "if c_k is not updated c_k = c(x_k) <: space_c\n"
00345     << L << "if mI > 0 and h_k is not updated h_k = h(x_k) <: space_h\n"
00346     << L << "if f_k is not updated f_k = f(x_k) <: REAL\n"
00347     ;
00348 }
00349 
00350 } // end namespace MoochoPack

Generated on Thu Sep 18 12:35:17 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1