NLPInterfacePack: C++ Interfaces and Implementation for Non-Linear Programs Version of the Day
NLPInterfacePack_NLPFirstDerivTester.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 <assert.h>
00043 #include <math.h>
00044 
00045 #include <iomanip>
00046 #include <sstream>
00047 #include <limits>
00048 
00049 #include "NLPInterfacePack_NLPFirstDerivTester.hpp"
00050 #include "NLPInterfacePack_NLPFirstOrder.hpp"
00051 #include "AbstractLinAlgPack_MatrixOp.hpp"
00052 #include "AbstractLinAlgPack_VectorMutable.hpp"
00053 #include "AbstractLinAlgPack_VectorOut.hpp"
00054 #include "AbstractLinAlgPack_VectorSpace.hpp"
00055 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00056 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00057 #include "AbstractLinAlgPack_EtaVector.hpp"
00058 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00059 #include "TestingHelperPack_update_success.hpp"
00060 #include "Teuchos_Assert.hpp"
00061 
00062 namespace {
00063 template< class T >
00064 inline
00065 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00066 } // end namespace
00067 
00068 namespace NLPInterfacePack {
00069 
00070 NLPFirstDerivTester::NLPFirstDerivTester(
00071   const calc_fd_prod_ptr_t  &calc_fd_prod
00072   ,ETestingMethod           fd_testing_method
00073   ,size_type                num_fd_directions
00074   ,value_type               warning_tol
00075   ,value_type               error_tol
00076   )
00077   :calc_fd_prod_(calc_fd_prod)
00078   ,fd_testing_method_(fd_testing_method)
00079   ,num_fd_directions_(num_fd_directions)
00080   ,warning_tol_(warning_tol)
00081   ,error_tol_(error_tol)
00082 {}
00083 
00084 bool NLPFirstDerivTester::finite_diff_check(
00085   NLP               *nlp
00086   ,const Vector     &xo
00087   ,const Vector     *xl
00088   ,const Vector     *xu
00089   ,const MatrixOp   *Gc
00090   ,const Vector     *Gf
00091   ,bool             print_all_warnings
00092   ,std::ostream     *out
00093   ) const
00094 {
00095   namespace rcp = MemMngPack;
00096   using AbstractLinAlgPack::assert_print_nan_inf;
00097 
00098   const size_type
00099     //n  = nlp->n(),
00100     m  = nlp->m();
00101 
00102   // ///////////////////////////////////
00103   // Validate the input
00104 
00105   TEUCHOS_TEST_FOR_EXCEPTION(
00106     !m && Gc, std::invalid_argument
00107     ,"NLPFirstDerivTester::finite_diff_check(...) : "
00108     "Error, Gc must be NULL if m == 0" );
00109 
00110   assert_print_nan_inf(xo, "xo",true,out);
00111   if(Gf)
00112     assert_print_nan_inf(*Gf, "Gf",true,out); 
00113 
00114   bool success = true;
00115 
00116   try {
00117 
00118   switch(fd_testing_method_) {
00119     case FD_COMPUTE_ALL:
00120       return fd_check_all(nlp,xo,xl,xu,Gc,Gf,print_all_warnings,out);
00121     case FD_DIRECTIONAL:
00122       return fd_directional_check(nlp,xo,xl,xu,Gc,Gf,print_all_warnings,out);
00123     default:
00124       TEUCHOS_TEST_FOR_EXCEPT(true);
00125   }
00126 
00127   } // end try
00128   catch( const AbstractLinAlgPack::NaNInfException& except ) {
00129     if(out)
00130       *out
00131         << "Error: found a NaN or Inf.  Stoping tests!\n";
00132     success = false;
00133   }
00134   
00135   return success; // will never be executed
00136 }
00137 
00138 // private
00139 
00140 bool NLPFirstDerivTester::fd_check_all(
00141   NLP               *nlp
00142   ,const Vector     &xo
00143   ,const Vector     *xl
00144   ,const Vector     *xu
00145   ,const MatrixOp   *Gc
00146   ,const Vector     *Gf
00147   ,bool             print_all_warnings
00148   ,std::ostream     *out
00149   ) const
00150 {
00151   using std::setw;
00152   using std::endl;
00153   using std::right;
00154 
00155   namespace rcp = MemMngPack;
00156   using AbstractLinAlgPack::sum;
00157   using AbstractLinAlgPack::dot;
00158   using AbstractLinAlgPack::Vp_StV;
00159   using AbstractLinAlgPack::assert_print_nan_inf;
00160   using AbstractLinAlgPack::random_vector;
00161   using LinAlgOpPack::V_StV;
00162   using LinAlgOpPack::V_MtV;
00163 
00164    using TestingHelperPack::update_success;
00165 
00166   bool success = true;
00167 
00168   const size_type
00169     n  = nlp->n();
00170   //m  = nlp->m();
00171 
00172   // //////////////////////////////////////////////
00173   // Validate the input
00174 
00175   NLP::vec_space_ptr_t
00176     space_x = nlp->space_x(),
00177     space_c = nlp->space_c();
00178 
00179   const CalcFiniteDiffProd
00180     &fd_deriv_prod = this->calc_fd_prod();
00181 
00182   const value_type
00183     //rand_y_l = -1.0, rand_y_u = 1.0,
00184     small_num = ::pow(std::numeric_limits<value_type>::epsilon(),0.25);
00185 
00186   if(out)
00187     *out
00188       << "\nComparing Gf and/or Gc with finite-difference values "
00189         "FDGf and/or FDGc one variable i at a time ...\n";
00190 
00191   value_type  max_Gf_warning_viol = 0.0;
00192   int         num_Gf_warning_viol = 0;
00193   value_type  max_Gc_warning_viol = 0.0;
00194   int         num_Gc_warning_viol = 0;
00195 
00196   VectorSpace::vec_mut_ptr_t
00197     e_i       = space_x->create_member(),
00198     Gc_i      = ( Gc ? space_c->create_member()  : Teuchos::null ),
00199     FDGc_i    = ( Gc ? space_c->create_member()  : Teuchos::null );
00200   *e_i = 0.0;
00201 
00202   for( size_type i = 1; i <= n; ++ i ) {
00203     EtaVector eta_i(i,n);
00204     e_i->set_ele(i,1.0);
00205     // Compute exact??? values
00206     value_type
00207       Gf_i = Gf ? Gf->get_ele(i) : 0.0;
00208     if(Gc)
00209       V_MtV( Gc_i.get(), *Gc, BLAS_Cpp::trans, eta_i() );
00210     // Compute finite difference values
00211     value_type
00212       FDGf_i;
00213     const bool preformed_fd = fd_deriv_prod.calc_deriv_product(
00214       xo,xl,xu
00215       ,*e_i
00216       ,NULL // fo
00217       ,NULL // co
00218       ,true // check_nan_inf
00219       ,nlp
00220       ,Gf ? &FDGf_i : NULL
00221       ,Gc ? FDGc_i.get() : NULL
00222       ,out
00223       );
00224     if( !preformed_fd ) {
00225       if(out)
00226         *out
00227           << "\nError, the finite difference computation was not preformed due to cramped bounds\n"
00228           << "Finite difference test failed!\n" << endl;
00229       return false;
00230     }
00231     
00232     // Compare the quantities
00233     // Gf
00234     assert_print_nan_inf(FDGf_i, "FDGf'*e_i",true,out);
00235     const value_type
00236       Gf_err = ::fabs( Gf_i - FDGf_i ) / ( ::fabs(Gf_i) + ::fabs(FDGf_i) + small_num );
00237     if(out)
00238       *out
00239         << "\nrel_err(Gf("<<i<<"),FDGf'*e("<<i<<")) = "
00240         << "rel_err(" << Gf_i << "," << FDGf_i << ") = "
00241         << Gf_err << endl;
00242     if( Gf_err >= warning_tol() ) {
00243       max_Gf_warning_viol = my_max( max_Gf_warning_viol, Gf_err );
00244       ++num_Gf_warning_viol;
00245     }
00246     if( Gf_err >= error_tol() ) {
00247       if(out)
00248         *out
00249           << "\nError, exceeded Gf_error_tol = " << error_tol() << endl
00250           << "Stoping the tests!\n";
00251       if(print_all_warnings)
00252         *out
00253           << "\ne_i =\n"     << *e_i
00254           << "\nGf_i =\n"    << Gf_i << std::endl
00255           << "\nFDGf_i =\n"  << FDGf_i << std::endl;
00256       update_success( false, &success );
00257       return false;
00258     }
00259     // Gc
00260     if(Gc) {
00261       const value_type
00262         sum_Gc_i   = sum(*Gc_i),
00263         sum_FDGc_i = sum(*FDGc_i);
00264       assert_print_nan_inf(sum_FDGc_i, "sum(FDGc'*e_i)",true,out);
00265       const value_type
00266         calc_err = ::fabs( ( sum_Gc_i - sum_FDGc_i )
00267                    /( ::fabs(sum_Gc_i) + ::fabs(sum_FDGc_i) + small_num ) );
00268       if(out)
00269         *out
00270           << "\nrel_err(sum(Gc'*e("<<i<<")),sum(FDGc'*e("<<i<<"))) = "
00271           << "rel_err(" << sum_Gc_i << "," << sum_FDGc_i << ") = "
00272           << calc_err << endl;
00273       if( calc_err >= warning_tol() ) {
00274         max_Gc_warning_viol = my_max( max_Gc_warning_viol, calc_err );
00275         ++num_Gc_warning_viol;
00276       }
00277       if( calc_err >= error_tol() ) {
00278         if(out)
00279           *out
00280             << "\nError, rel_err(sum(Gc'*e("<<i<<")),sum(FDGc'*e("<<i<<"))) = "
00281             << "rel_err(" << sum_Gc_i << "," << sum_FDGc_i << ") = "
00282             << calc_err << endl
00283             << "exceeded error_tol = " << error_tol() << endl
00284             << "Stoping the tests!\n";
00285         if(print_all_warnings)
00286           *out << "\ne_i =\n"     << *e_i
00287              << "\nGc_i =\n"    << *Gc_i
00288              << "\nFDGc_i =\n"  << *FDGc_i;
00289         update_success( false, &success );
00290         return false;
00291       }
00292       if( calc_err >= warning_tol() ) {
00293         if(out)
00294           *out
00295             << "\nWarning, rel_err(sum(Gc'*e("<<i<<")),sum(FDGc'*e("<<i<<"))) = "
00296             << "rel_err(" << sum_Gc_i << "," << sum_FDGc_i << ") = "
00297             << calc_err << endl
00298             << "exceeded warning_tol = " << warning_tol() << endl;
00299       }
00300     }
00301     e_i->set_ele(i,0.0);
00302   }
00303   if(out && num_Gf_warning_viol)
00304     *out
00305       << "\nFor Gf, there were " << num_Gf_warning_viol << " warning tolerance "
00306       << "violations out of num_fd_directions = " << num_fd_directions()
00307       << " computations of FDGf'*e(i)\n"
00308       << "and the maximum violation was " << max_Gf_warning_viol
00309       << " > Gf_waring_tol = " << warning_tol() << endl;
00310   if(out && num_Gc_warning_viol)
00311     *out
00312       << "\nFor Gc, there were " << num_Gc_warning_viol << " warning tolerance "
00313       << "violations out of num_fd_directions = " << num_fd_directions()
00314       << " computations of FDGc'*e(i)\n"
00315       << "and the maximum violation was " << max_Gc_warning_viol
00316       << " > Gc_waring_tol = " << warning_tol() << endl;
00317   if(out)
00318     *out
00319       << "\nCongradulations!  All of the computed errors were within the specified error tolerance!\n";
00320 
00321   return true;
00322 
00323 }
00324 
00325 bool NLPFirstDerivTester::fd_directional_check(
00326   NLP               *nlp
00327   ,const Vector     &xo
00328   ,const Vector     *xl
00329   ,const Vector     *xu
00330   ,const MatrixOp   *Gc
00331   ,const Vector     *Gf
00332   ,bool             print_all_warnings
00333   ,std::ostream     *out
00334   ) const
00335 {
00336   using std::setw;
00337   using std::endl;
00338   using std::right;
00339 
00340   namespace rcp = MemMngPack;
00341   using AbstractLinAlgPack::sum;
00342   using AbstractLinAlgPack::dot;
00343   using AbstractLinAlgPack::Vp_StV;
00344   using AbstractLinAlgPack::assert_print_nan_inf;
00345   using AbstractLinAlgPack::random_vector;
00346   using LinAlgOpPack::V_StV;
00347   using LinAlgOpPack::V_MtV;
00348 
00349    using TestingHelperPack::update_success;
00350 
00351   bool success = true;
00352 
00353   //const size_type
00354   //n  = nlp->n(),
00355   //m  = nlp->m();
00356 
00357   // //////////////////////////////////////////////
00358   // Validate the input
00359 
00360   NLP::vec_space_ptr_t
00361     space_x = nlp->space_x(),
00362     space_c = nlp->space_c();
00363 
00364   const CalcFiniteDiffProd
00365     &fd_deriv_prod = this->calc_fd_prod();
00366 
00367   const value_type
00368     rand_y_l = -1.0, rand_y_u = 1.0,
00369     small_num = ::pow(std::numeric_limits<value_type>::epsilon(),0.25);
00370 
00371   if(out)
00372     *out
00373       << "\nComparing directional products Gf'*y and/or Gc'*y with finite difference values "
00374         " FDGf'*y and/or FDGc'*y for random y's ...\n";
00375 
00376   value_type  max_Gf_warning_viol = 0.0;
00377   int         num_Gf_warning_viol = 0;
00378 
00379   VectorSpace::vec_mut_ptr_t
00380     y         = space_x->create_member(),
00381     Gc_prod   = ( Gc ? space_c->create_member()  : Teuchos::null ),
00382     FDGc_prod = ( Gc ? space_c->create_member()  : Teuchos::null );
00383 
00384   const int num_fd_directions_used = ( num_fd_directions() > 0 ? num_fd_directions() : 1 );
00385 
00386   for( int direc_i = 1; direc_i <= num_fd_directions_used; ++direc_i ) {
00387     if( num_fd_directions() > 0 ) {
00388       random_vector( rand_y_l, rand_y_u, y.get() );
00389       if(out)
00390         *out
00391           << "\n****"
00392           << "\n**** Random directional vector " << direc_i << " ( ||y||_1 / n = "
00393           << (y->norm_1() / y->dim()) << " )"
00394           << "\n***\n";
00395     }
00396     else {
00397       *y = 1.0;
00398       if(out)
00399         *out
00400           << "\n****"
00401           << "\n**** Ones vector y ( ||y||_1 / n = "<<(y->norm_1()/y->dim())<<" )"
00402           << "\n***\n";
00403     }
00404     // Compute exact??? values
00405     value_type
00406       Gf_y = Gf ? dot( *Gf, *y ) : 0.0;
00407     if(Gc)
00408       V_MtV( Gc_prod.get(), *Gc, BLAS_Cpp::trans, *y );
00409     // Compute finite difference values
00410     value_type
00411       FDGf_y;
00412     const bool preformed_fd = fd_deriv_prod.calc_deriv_product(
00413       xo,xl,xu
00414       ,*y
00415       ,NULL // fo
00416       ,NULL // co
00417       ,true // check_nan_inf
00418       ,nlp
00419       ,Gf ? &FDGf_y : NULL
00420       ,Gc ? FDGc_prod.get() : NULL
00421       ,out
00422       );
00423     if( !preformed_fd ) {
00424       if(out)
00425         *out
00426           << "\nError, the finite difference computation was not preformed due to cramped bounds\n"
00427           << "Finite difference test failed!\n" << endl;
00428       return false;
00429     }
00430     
00431     // Compare the quantities
00432     // Gf
00433     assert_print_nan_inf(FDGf_y, "FDGf'*y",true,out);
00434     const value_type
00435       Gf_err = ::fabs( Gf_y - FDGf_y ) / ( ::fabs(Gf_y) + ::fabs(FDGf_y) + small_num );
00436     if(out)
00437       *out
00438         << "\nrel_err(Gf'*y,FDGf'*y) = "
00439         << "rel_err(" << Gf_y << "," << FDGf_y << ") = "
00440         << Gf_err << endl;
00441     if( Gf_err >= warning_tol() ) {
00442       max_Gf_warning_viol = my_max( max_Gf_warning_viol, Gf_err );
00443       ++num_Gf_warning_viol;
00444     }
00445     if( Gf_err >= error_tol() ) {
00446       if(out)
00447         *out
00448           << "\nError, exceeded Gf_error_tol = " << error_tol() << endl
00449           << "Stoping the tests!\n";
00450       return false;
00451     }
00452     // Gc
00453     if(Gc) {
00454       const value_type
00455         sum_Gc_prod   = sum(*Gc_prod),
00456         sum_FDGc_prod = sum(*FDGc_prod);
00457       assert_print_nan_inf(sum_FDGc_prod, "sum(FDGc'*y)",true,out);
00458       const value_type
00459         calc_err = ::fabs( ( sum_Gc_prod - sum_FDGc_prod )
00460                    /( ::fabs(sum_Gc_prod) + ::fabs(sum_FDGc_prod) + small_num ) );
00461       if(out)
00462         *out
00463           << "\nrel_err(sum(Gc'*y),sum(FDGc'*y)) = "
00464           << "rel_err(" << sum_Gc_prod << "," << sum_FDGc_prod << ") = "
00465           << calc_err << endl;
00466       if( calc_err >= error_tol() ) {
00467         if(out)
00468           *out
00469             << "\nError, rel_err(sum(Gc'*y),sum(FDGc'*y)) = "
00470             << "rel_err(" << sum_Gc_prod << "," << sum_FDGc_prod << ") = "
00471             << calc_err << endl
00472             << "exceeded error_tol = " << error_tol() << endl
00473             << "Stoping the tests!\n";
00474         if(print_all_warnings)
00475           *out << "\ny =\n"          << *y
00476              << "\nGc_prod =\n"    << *Gc_prod
00477              << "\nFDGc_prod =\n"  << *FDGc_prod;
00478         update_success( false, &success );
00479         return false;
00480       }
00481       if( calc_err >= warning_tol() ) {
00482         if(out)
00483           *out
00484             << "\nWarning, rel_err(sum(Gc'*y),sum(FDGc'*y)) = "
00485             << "rel_err(" << sum_Gc_prod << "," << sum_FDGc_prod << ") = "
00486             << calc_err << endl
00487             << "exceeded warning_tol = " << warning_tol() << endl;
00488       }
00489     }
00490   }
00491   if(out && num_Gf_warning_viol)
00492     *out
00493       << "\nFor Gf, there were " << num_Gf_warning_viol << " warning tolerance "
00494       << "violations out of num_fd_directions = " << num_fd_directions()
00495       << " computations of FDGf'*y\n"
00496       << "and the maximum violation was " << max_Gf_warning_viol
00497       << " > Gf_waring_tol = " << warning_tol() << endl;
00498 
00499   if(out)
00500     *out
00501       << "\nCongradulations!  All of the computed errors were within the specified error tolerance!\n";
00502 
00503   return true;
00504 }
00505 
00506 } // end namespace NLPInterfacePack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends