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

Generated on Tue Oct 20 12:51:46 2009 for MOOCHO (Single Doxygen Collection) by doxygen 1.4.7