MoochoPack : Framework for Large-Scale Optimization Algorithms Version of the Day
MoochoPack_EvalNewPointTailoredApproach_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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <ostream>
00043 
00044 #include "MoochoPack_EvalNewPointTailoredApproach_Step.hpp"
00045 #include "MoochoPack_Exceptions.hpp"
00046 #include "MoochoPack_moocho_algo_conversion.hpp"
00047 #include "IterationPack_print_algorithm_step.hpp"
00048 #include "ConstrainedOptPack_MatrixIdentConcatStd.hpp"
00049 #include "NLPInterfacePack_NLPDirect.hpp"
00050 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00051 #include "AbstractLinAlgPack_VectorMutable.hpp"
00052 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00053 #include "AbstractLinAlgPack_VectorOut.hpp"
00054 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00055 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00056 #include "Teuchos_dyn_cast.hpp"
00057 #include "Teuchos_Assert.hpp"
00058 
00059 namespace MoochoPack {
00060 
00061 EvalNewPointTailoredApproach_Step::EvalNewPointTailoredApproach_Step(
00062   const deriv_tester_ptr_t      &deriv_tester
00063   ,const bounds_tester_ptr_t    &bounds_tester
00064   , EFDDerivTesting             fd_deriv_testing
00065   )
00066   :deriv_tester_(deriv_tester)
00067   ,bounds_tester_(bounds_tester)
00068   ,fd_deriv_testing_(fd_deriv_testing)
00069 {}
00070 
00071 bool EvalNewPointTailoredApproach_Step::do_step(
00072   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00073   ,poss_type assoc_step_poss
00074   )
00075 {
00076 
00077   using Teuchos::RCP;
00078   using Teuchos::rcp;
00079   using Teuchos::dyn_cast;
00080   using AbstractLinAlgPack::assert_print_nan_inf;
00081   using LinAlgOpPack::V_MtV;
00082   using IterationPack::print_algorithm_step;
00083 
00084   NLPAlgo             &algo   = rsqp_algo(_algo);
00085   NLPAlgoState        &s      = algo.rsqp_state();
00086   NLPDirect           &nlp    = dyn_cast<NLPDirect>(algo.nlp());
00087 
00088   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00089   EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_journal_output_level();
00090   std::ostream& out = algo.track().journal_out();
00091 
00092   // print step header.
00093   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00094     using IterationPack::print_algorithm_step;
00095     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00096   }
00097 
00098   if(!nlp.is_initialized())
00099     nlp.initialize(algo.algo_cntr().check_results());
00100 
00101   Teuchos::VerboseObjectTempState<NLP>
00102     nlpOutputTempState(
00103       rcp(&nlp,false), Teuchos::getFancyOStream(rcp(&out,false)),
00104       convertToVerbLevel(olevel) );
00105 
00106   const Range1D
00107     var_dep = nlp.var_dep(),
00108     var_indep = nlp.var_indep();
00109 
00110   s.var_dep(var_dep);
00111   s.var_indep(var_indep);
00112 
00113   const size_type
00114     n  = nlp.n(),
00115     m  = nlp.m(),
00116     r  = var_dep.size();
00117 
00118   TEUCHOS_TEST_FOR_EXCEPTION(
00119     m > r, TestFailed
00120     ,"EvalNewPointTailoredApproach_Step::do_step(...) : Error, "
00121     "Undecomposed equalities are supported yet!" );
00122 
00123   IterQuantityAccess<VectorMutable>
00124     &x_iq = s.x();
00125 
00126   if( x_iq.last_updated() == IterQuantity::NONE_UPDATED ) {
00127     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00128       out << "\nx is not updated for any k so set x_k = nlp.xinit() ...\n";
00129     }
00130     x_iq.set_k(0) = nlp.xinit();
00131   }
00132 
00133   // Validate x
00134   if(algo.algo_cntr().check_results()) {
00135     assert_print_nan_inf(
00136       x_iq.get_k(0), "x_k", true
00137       , int(olevel) >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL
00138       );
00139     if( nlp.num_bounded_x() > 0 ) {
00140       if(!bounds_tester().check_in_bounds(
00141            int(olevel)  >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL
00142            ,int(olevel) >= int(PRINT_VECTORS)                // print_all_warnings
00143            ,int(olevel) >= int(PRINT_ITERATION_QUANTITIES)  // print_vectors
00144            ,nlp.xl(),        "xl"
00145            ,nlp.xu(),        "xu"
00146            ,x_iq.get_k(0),   "x_k"
00147            ))
00148       {
00149         TEUCHOS_TEST_FOR_EXCEPTION(
00150           true, TestFailed
00151           ,"EvalNewPointTailoredApproach_Step::do_step(...) : Error, "
00152           "the variables bounds xl <= x_k <= xu where violated!" );
00153       }
00154     }
00155   }
00156 
00157   Vector &x = x_iq.get_k(0);
00158 
00159   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00160     out << "\n||x_k||inf            = " << x.norm_inf();
00161     if( var_dep.size() )
00162       out << "\n||x(var_dep)_k||inf   = " << x.sub_view(var_dep)->norm_inf();
00163     if( var_indep.size() )
00164       out << "\n||x(var_indep)_k||inf = " << x.sub_view(var_indep)->norm_inf();
00165     out << std::endl;
00166   }
00167   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00168     out << "\nx_k = \n" << x;
00169     if( var_dep.size() )
00170       out << "\nx(var_dep)_k = \n" << *x.sub_view(var_dep);
00171   }
00172   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00173     if( var_indep.size() )
00174       out << "\nx(var_indep)_k = \n" << *x.sub_view(var_indep);
00175   }
00176 
00177   // If c_k is not updated then we must compute it
00178   bool recalc_c = true;
00179 
00180   if( !s.c().updated_k(0) ) {
00181     s.c().set_k(0);
00182     recalc_c = true;
00183   }
00184   else {
00185     recalc_c = false;
00186   }
00187     
00188   // Get references to Z, Y, Uz and Uy
00189   MatrixOp
00190     &Z_k  = s.Z().set_k(0),
00191     &Y_k  = s.Y().set_k(0),
00192     *Uz_k = (m > r) ? &s.Uz().set_k(0) : NULL,
00193     *Uy_k = (m > r) ? &s.Uy().set_k(0) : NULL;
00194   MatrixIdentConcatStd
00195     &cZ_k = dyn_cast<MatrixIdentConcatStd>(Z_k);
00196   // Release any references to D in Y or Uy
00197    uninitialize_Y_Uy(&Y_k,Uy_k);
00198   // If Z has not been initialized or Z.D is being shared by someone else we need to reconstruct Z.D
00199   bool reconstruct_Z_D = (cZ_k.rows() == n || cZ_k.cols() != n-r || cZ_k.D_ptr().count() > 1);
00200   MatrixIdentConcatStd::D_ptr_t
00201     D_ptr = Teuchos::null;
00202   if( reconstruct_Z_D )
00203     D_ptr = nlp.factory_D()->create();
00204   else
00205     D_ptr = cZ_k.D_ptr();
00206 
00207   // Compute all the quantities.
00208   const bool supports_Gf = nlp.supports_Gf();
00209   Teuchos::RCP<MatrixOp>
00210     GcU = (m > r) ? nlp.factory_GcU()->create() : Teuchos::null; // ToDo: Reuse GcU somehow? 
00211   VectorMutable
00212     &py_k  = s.py().set_k(0);
00213   nlp.unset_quantities();
00214   nlp.calc_point(
00215     x                                                     // x
00216     ,!s.f().updated_k(0) ? &s.f().set_k(0) : NULL         // f
00217     ,&s.c().get_k(0)                                      // c
00218     ,recalc_c                                             // recalc_c
00219     ,supports_Gf?&s.Gf().set_k(0):NULL                    // Gf
00220     ,&py_k                                                // -inv(C)*c
00221     ,&s.rGf().set_k(0)                                    // rGf
00222     ,GcU.get()                                            // GcU
00223     ,const_cast<MatrixOp*>(D_ptr.get())                   // -inv(C)*N
00224     ,Uz_k                                                 // Uz
00225     );
00226   s.equ_decomp(   nlp.con_decomp()   );
00227   s.equ_undecomp( nlp.con_undecomp() );
00228 
00229   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00230     out << "\nQuantities computed directly from NLPDirect nlp object ...\n";
00231     out << "\nf_k                    = " << s.f().get_k(0);
00232     out << "\n||c_k||inf             = " << s.c().get_k(0).norm_inf();
00233     if(supports_Gf) {
00234       const Vector &Gf = s.Gf().get_k(0);
00235       out << "\n||Gf_k||inf            = " << Gf.norm_inf();
00236       if (var_dep.size())
00237         out << "\n||Gf(var_dep)_k||inf   = " << Gf.sub_view(var_dep)->norm_inf();
00238       if (var_indep.size())
00239         out << "\n||Gf(var_indep)_k||inf = " << Gf.sub_view(var_indep)->norm_inf();
00240     }
00241     out << "\n||py_k||inf            = " << s.py().get_k(0).norm_inf();
00242     out << "\n||rGf_k||inf           = " << s.rGf().get_k(0).norm_inf();
00243     out << std::endl;
00244   }
00245 
00246   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00247     out << "\nrGf_k = \n" << s.rGf().get_k(0);
00248     out << std::endl;
00249   }
00250 
00251   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00252     out << "\nc_k  = \n" << s.c().get_k(0);
00253     if(supports_Gf)
00254       out << "\nGf_k = \n" << s.Gf().get_k(0);
00255     out << "\npy_k = \n" << s.py().get_k(0);
00256     out << std::endl;
00257   }
00258 
00259   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00260     out << "\nD = -inv(C)*N = \n" << *D_ptr;
00261     out << std::endl;
00262   }
00263 
00264   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00265     out << "Printing column norms of D = -inv(C)*N:\n";
00266     RCP<VectorMutable>
00267       e_i = D_ptr->space_rows().create_member();
00268     RCP<VectorMutable>
00269       D_i = D_ptr->space_cols().create_member();
00270     *e_i = 0.0;
00271     for( int i = 1; i <= (n-r); ++i ) {
00272       e_i->set_ele(i,1.0);
00273       V_MtV( &*D_i, *D_ptr, BLAS_Cpp::no_trans, *e_i );
00274       e_i->set_ele(i,0.0);
00275       out << "  ||D(:,"<<i<<")||_2 = " << D_i->norm_2() << "\n";
00276     }
00277   }
00278 
00279   if(algo.algo_cntr().check_results()) {
00280     assert_print_nan_inf(s.f().get_k(0),   "f_k",true,&out); 
00281     assert_print_nan_inf(s.c().get_k(0),   "c_k",true,&out); 
00282     if(supports_Gf)
00283       assert_print_nan_inf(s.Gf().get_k(0),  "Gf_k",true,&out); 
00284     assert_print_nan_inf(s.py().get_k(0),  "py_k",true,&out);
00285     assert_print_nan_inf(s.rGf().get_k(0), "rGf_k",true,&out);
00286   }
00287 
00288   // Check the derivatives if we are checking the results
00289   if(   fd_deriv_testing() == FD_TEST
00290     || ( fd_deriv_testing() == FD_DEFAULT && algo.algo_cntr().check_results() )  )
00291   {
00292     
00293     if( olevel >= PRINT_ALGORITHM_STEPS ) {
00294       out << "\n*** Checking derivatives by finite differences\n";
00295     }
00296     
00297     const bool has_bounds = nlp.num_bounded_x() > 0;
00298     const bool nlp_passed = deriv_tester().finite_diff_check(
00299       &nlp
00300       ,x
00301       ,has_bounds ? &nlp.xl() : (const Vector*)NULL
00302       ,has_bounds ? &nlp.xu() : (const Vector*)NULL
00303       ,&s.c().get_k(0)
00304       ,supports_Gf?&s.Gf().get_k(0):NULL
00305       ,&s.py().get_k(0)
00306       ,&s.rGf().get_k(0)
00307       ,GcU.get()
00308       ,D_ptr.get()
00309       ,Uz_k
00310       ,olevel >= PRINT_VECTORS
00311       ,( olevel >= PRINT_ALGORITHM_STEPS ) ? &out : (std::ostream*)NULL
00312       );
00313     TEUCHOS_TEST_FOR_EXCEPTION(
00314       !nlp_passed, TestFailed
00315       ,"EvalNewPointTailoredApproach_Step::do_step(...) : Error, "
00316       "the tests of the nlp derivatives failed!" );
00317   }
00318 
00319   if( reconstruct_Z_D ) {
00320     //
00321     // Z = [      D     ] space_xD
00322     //     [      I     ] space_xI
00323     //        space_xI
00324     //
00325     cZ_k.initialize(
00326       nlp.space_x()                                          // space_cols
00327       ,nlp.space_x()->sub_space(nlp.var_indep())->clone()    // space_rows
00328       ,MatrixIdentConcatStd::TOP                             // top_or_bottom
00329       ,1.0                                                   // alpha
00330       ,D_ptr                                                 // D_ptr
00331       ,BLAS_Cpp::no_trans                                    // D_trans
00332       );
00333   }
00334 
00335   // Compute py, Y and Uy
00336   if( olevel >= PRINT_ALGORITHM_STEPS )
00337     out << "\nUpdating py_k, Y_k, and Uy_k given D_k ...\n";
00338   calc_py_Y_Uy( nlp, D_ptr, &py_k, &Y_k, Uy_k, olevel, out ); 
00339 
00340   // Compute Ypy = Y*py
00341   VectorMutable
00342     &Ypy_k  = s.Ypy().set_k(0);
00343   V_MtV( &Ypy_k, Y_k, BLAS_Cpp::no_trans, py_k );
00344 
00345   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00346     out << "\nQuantities after computing final py and Ypy ...\n";
00347     out << "\n||py_k||inf   = " << s.py().get_k(0).norm_inf();
00348     out << "\n||Ypy_k||inf  = " << s.Ypy().get_k(0).norm_inf();
00349     out << std::endl;
00350   }
00351 
00352   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00353     out << "\npy_k = \n"  << s.py().get_k(0);
00354   }
00355 
00356   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00357     if(var_indep.size())
00358       out << "\nYpy(var_indep)_k = \n" << *s.Ypy().get_k(0).sub_view(var_indep);
00359     out << std::endl;
00360   }
00361 
00362   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00363     if(var_dep.size())
00364       out << "\nYpy(var_dep)_k = \n" << *s.Ypy().get_k(0).sub_view(var_dep);
00365     out << "\nYpy_k = \n" << s.Ypy().get_k(0);
00366     out << std::endl;
00367   }
00368 
00369   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00370     out << "\nZ_k = \n" << s.Z().get_k(0);
00371     out << "\nY_k = \n" << s.Y().get_k(0);
00372     out << std::endl;
00373   }
00374 
00375   return true;
00376 
00377 }
00378 
00379 void EvalNewPointTailoredApproach_Step::print_step(
00380   const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type
00381   ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
00382   ) const
00383 {
00384   out
00385     << L << "*** Evaluate the new point for the \"Tailored Approach\"\n"
00386     << L << "if nlp is not initialized then initialize the nlp\n"
00387     << L << "if x is not updated for any k then set x_k = nlp.xinit\n"
00388     << L << "if Gf is supported then set Gf_k = Gf(x_k) <: space_x\n"
00389     << L << "For Gc(:,equ_decomp) = [ C' ; N' ] <: space_x|space_c(equ_decomp) compute:\n"
00390     << L << "  py_k = -inv(C)*c_k\n"
00391     << L << "  D = -inv(C)*N <: R^(n x (n-m))\n"
00392     << L << "  rGf_k = D'*Gf_k(var_dep) + Gf_k(var_indep)\n"
00393     << L << "  Z_k = [ D ; I ] <: R^(n x (n-m))\n"
00394     << L << "  if m > r then\n"
00395     << L << "    Uz_k = Gc(var_indep,equ_undecomp)' + Gc(var_dep,equ_undecomp)'*D\n"
00396     << L << "  end\n"
00397     << L << "if ( (fd_deriv_testing==FD_TEST)\n"
00398     << L << "    or (fd_deriv_testing==FD_DEFAULT and check_results==true)\n"
00399     << L << "  ) then\n"
00400     << L << "  check Gf_k, py_k, rGf_k, D, Uz (if m > r) and Vz (if mI > 0) by finite differences.\n"
00401     << L << "end\n";
00402   print_calc_py_Y_Uy( out, L );
00403   out
00404     << L << "Ypy_k = Y_k * py_k\n"
00405     << L << "if c_k is not updated c_k = c(x_k) <: space_c\n"
00406     << L << "if mI > 0 and h_k is not updated h_k = h(x_k) <: space_h\n"
00407     << L << "if f_k is not updated f_k = f(x_k) <: REAL\n"
00408     ;
00409 }
00410 
00411 } // end namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends