NLPInterfacePack: C++ Interfaces and Implementation for Non-Linear Programs  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
NLPInterfacePack_NLPFirstDerivTester.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <assert.h>
43 #include <math.h>
44 
45 #include <iomanip>
46 #include <sstream>
47 #include <limits>
48 
49 #include "NLPInterfacePack_NLPFirstDerivTester.hpp"
50 #include "NLPInterfacePack_NLPFirstOrder.hpp"
51 #include "AbstractLinAlgPack_MatrixOp.hpp"
52 #include "AbstractLinAlgPack_VectorMutable.hpp"
53 #include "AbstractLinAlgPack_VectorOut.hpp"
54 #include "AbstractLinAlgPack_VectorSpace.hpp"
55 #include "AbstractLinAlgPack_VectorStdOps.hpp"
56 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
57 #include "AbstractLinAlgPack_EtaVector.hpp"
58 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
59 #include "TestingHelperPack_update_success.hpp"
60 #include "Teuchos_Assert.hpp"
61 
62 namespace {
63 template< class T >
64 inline
65 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
66 } // end namespace
67 
68 namespace NLPInterfacePack {
69 
71  const calc_fd_prod_ptr_t &calc_fd_prod
72  ,ETestingMethod fd_testing_method
73  ,size_type num_fd_directions
74  ,value_type warning_tol
75  ,value_type error_tol
76  )
77  :calc_fd_prod_(calc_fd_prod)
78  ,fd_testing_method_(fd_testing_method)
79  ,num_fd_directions_(num_fd_directions)
80  ,warning_tol_(warning_tol)
81  ,error_tol_(error_tol)
82 {}
83 
85  NLP *nlp
86  ,const Vector &xo
87  ,const Vector *xl
88  ,const Vector *xu
89  ,const MatrixOp *Gc
90  ,const Vector *Gf
91  ,bool print_all_warnings
92  ,std::ostream *out
93  ) const
94 {
95  namespace rcp = MemMngPack;
96  using AbstractLinAlgPack::assert_print_nan_inf;
97 
98  const size_type
99  //n = nlp->n(),
100  m = nlp->m();
101 
102  // ///////////////////////////////////
103  // Validate the input
104 
106  !m && Gc, std::invalid_argument
107  ,"NLPFirstDerivTester::finite_diff_check(...) : "
108  "Error, Gc must be NULL if m == 0" );
109 
110  assert_print_nan_inf(xo, "xo",true,out);
111  if(Gf)
112  assert_print_nan_inf(*Gf, "Gf",true,out);
113 
114  bool success = true;
115 
116  try {
117 
118  switch(fd_testing_method_) {
119  case FD_COMPUTE_ALL:
120  return fd_check_all(nlp,xo,xl,xu,Gc,Gf,print_all_warnings,out);
121  case FD_DIRECTIONAL:
122  return fd_directional_check(nlp,xo,xl,xu,Gc,Gf,print_all_warnings,out);
123  default:
125  }
126 
127  } // end try
128  catch( const AbstractLinAlgPack::NaNInfException& except ) {
129  if(out)
130  *out
131  << "Error: found a NaN or Inf. Stoping tests!\n";
132  success = false;
133  }
134 
135  return success; // will never be executed
136 }
137 
138 // private
139 
140 bool NLPFirstDerivTester::fd_check_all(
141  NLP *nlp
142  ,const Vector &xo
143  ,const Vector *xl
144  ,const Vector *xu
145  ,const MatrixOp *Gc
146  ,const Vector *Gf
147  ,bool print_all_warnings
148  ,std::ostream *out
149  ) const
150 {
151  using std::setw;
152  using std::endl;
153  using std::right;
154 
155  namespace rcp = MemMngPack;
159  using AbstractLinAlgPack::assert_print_nan_inf;
161  using LinAlgOpPack::V_StV;
162  using LinAlgOpPack::V_MtV;
163 
164  using TestingHelperPack::update_success;
165 
166  bool success = true;
167 
168  const size_type
169  n = nlp->n();
170  //m = nlp->m();
171 
172  // //////////////////////////////////////////////
173  // Validate the input
174 
176  space_x = nlp->space_x(),
177  space_c = nlp->space_c();
178 
179  const CalcFiniteDiffProd
180  &fd_deriv_prod = this->calc_fd_prod();
181 
182  const value_type
183  //rand_y_l = -1.0, rand_y_u = 1.0,
184  small_num = ::pow(std::numeric_limits<value_type>::epsilon(),0.25);
185 
186  if(out)
187  *out
188  << "\nComparing Gf and/or Gc with finite-difference values "
189  "FDGf and/or FDGc one variable i at a time ...\n";
190 
191  value_type max_Gf_warning_viol = 0.0;
192  int num_Gf_warning_viol = 0;
193  value_type max_Gc_warning_viol = 0.0;
194  int num_Gc_warning_viol = 0;
195 
196  VectorSpace::vec_mut_ptr_t
197  e_i = space_x->create_member(),
198  Gc_i = ( Gc ? space_c->create_member() : Teuchos::null ),
199  FDGc_i = ( Gc ? space_c->create_member() : Teuchos::null );
200  *e_i = 0.0;
201 
202  for( size_type i = 1; i <= n; ++ i ) {
203  EtaVector eta_i(i,n);
204  e_i->set_ele(i,1.0);
205  // Compute exact??? values
206  value_type
207  Gf_i = Gf ? Gf->get_ele(i) : 0.0;
208  if(Gc)
209  V_MtV( Gc_i.get(), *Gc, BLAS_Cpp::trans, eta_i() );
210  // Compute finite difference values
211  value_type
212  FDGf_i;
213  const bool preformed_fd = fd_deriv_prod.calc_deriv_product(
214  xo,xl,xu
215  ,*e_i
216  ,NULL // fo
217  ,NULL // co
218  ,true // check_nan_inf
219  ,nlp
220  ,Gf ? &FDGf_i : NULL
221  ,Gc ? FDGc_i.get() : NULL
222  ,out
223  );
224  if( !preformed_fd ) {
225  if(out)
226  *out
227  << "\nError, the finite difference computation was not preformed due to cramped bounds\n"
228  << "Finite difference test failed!\n" << endl;
229  return false;
230  }
231 
232  // Compare the quantities
233  // Gf
234  assert_print_nan_inf(FDGf_i, "FDGf'*e_i",true,out);
235  const value_type
236  Gf_err = ::fabs( Gf_i - FDGf_i ) / ( ::fabs(Gf_i) + ::fabs(FDGf_i) + small_num );
237  if(out)
238  *out
239  << "\nrel_err(Gf("<<i<<"),FDGf'*e("<<i<<")) = "
240  << "rel_err(" << Gf_i << "," << FDGf_i << ") = "
241  << Gf_err << endl;
242  if( Gf_err >= warning_tol() ) {
243  max_Gf_warning_viol = my_max( max_Gf_warning_viol, Gf_err );
244  ++num_Gf_warning_viol;
245  }
246  if( Gf_err >= error_tol() ) {
247  if(out)
248  *out
249  << "\nError, exceeded Gf_error_tol = " << error_tol() << endl
250  << "Stoping the tests!\n";
251  if(print_all_warnings)
252  *out
253  << "\ne_i =\n" << *e_i
254  << "\nGf_i =\n" << Gf_i << std::endl
255  << "\nFDGf_i =\n" << FDGf_i << std::endl;
256  update_success( false, &success );
257  return false;
258  }
259  // Gc
260  if(Gc) {
261  const value_type
262  sum_Gc_i = sum(*Gc_i),
263  sum_FDGc_i = sum(*FDGc_i);
264  assert_print_nan_inf(sum_FDGc_i, "sum(FDGc'*e_i)",true,out);
265  const value_type
266  calc_err = ::fabs( ( sum_Gc_i - sum_FDGc_i )
267  /( ::fabs(sum_Gc_i) + ::fabs(sum_FDGc_i) + small_num ) );
268  if(out)
269  *out
270  << "\nrel_err(sum(Gc'*e("<<i<<")),sum(FDGc'*e("<<i<<"))) = "
271  << "rel_err(" << sum_Gc_i << "," << sum_FDGc_i << ") = "
272  << calc_err << endl;
273  if( calc_err >= warning_tol() ) {
274  max_Gc_warning_viol = my_max( max_Gc_warning_viol, calc_err );
275  ++num_Gc_warning_viol;
276  }
277  if( calc_err >= error_tol() ) {
278  if(out)
279  *out
280  << "\nError, rel_err(sum(Gc'*e("<<i<<")),sum(FDGc'*e("<<i<<"))) = "
281  << "rel_err(" << sum_Gc_i << "," << sum_FDGc_i << ") = "
282  << calc_err << endl
283  << "exceeded error_tol = " << error_tol() << endl
284  << "Stoping the tests!\n";
285  if(print_all_warnings)
286  *out << "\ne_i =\n" << *e_i
287  << "\nGc_i =\n" << *Gc_i
288  << "\nFDGc_i =\n" << *FDGc_i;
289  update_success( false, &success );
290  return false;
291  }
292  if( calc_err >= warning_tol() ) {
293  if(out)
294  *out
295  << "\nWarning, rel_err(sum(Gc'*e("<<i<<")),sum(FDGc'*e("<<i<<"))) = "
296  << "rel_err(" << sum_Gc_i << "," << sum_FDGc_i << ") = "
297  << calc_err << endl
298  << "exceeded warning_tol = " << warning_tol() << endl;
299  }
300  }
301  e_i->set_ele(i,0.0);
302  }
303  if(out && num_Gf_warning_viol)
304  *out
305  << "\nFor Gf, there were " << num_Gf_warning_viol << " warning tolerance "
306  << "violations out of num_fd_directions = " << num_fd_directions()
307  << " computations of FDGf'*e(i)\n"
308  << "and the maximum violation was " << max_Gf_warning_viol
309  << " > Gf_waring_tol = " << warning_tol() << endl;
310  if(out && num_Gc_warning_viol)
311  *out
312  << "\nFor Gc, there were " << num_Gc_warning_viol << " warning tolerance "
313  << "violations out of num_fd_directions = " << num_fd_directions()
314  << " computations of FDGc'*e(i)\n"
315  << "and the maximum violation was " << max_Gc_warning_viol
316  << " > Gc_waring_tol = " << warning_tol() << endl;
317  if(out)
318  *out
319  << "\nCongradulations! All of the computed errors were within the specified error tolerance!\n";
320 
321  return true;
322 
323 }
324 
325 bool NLPFirstDerivTester::fd_directional_check(
326  NLP *nlp
327  ,const Vector &xo
328  ,const Vector *xl
329  ,const Vector *xu
330  ,const MatrixOp *Gc
331  ,const Vector *Gf
332  ,bool print_all_warnings
333  ,std::ostream *out
334  ) const
335 {
336  using std::setw;
337  using std::endl;
338  using std::right;
339 
340  namespace rcp = MemMngPack;
344  using AbstractLinAlgPack::assert_print_nan_inf;
346  using LinAlgOpPack::V_StV;
347  using LinAlgOpPack::V_MtV;
348 
349  using TestingHelperPack::update_success;
350 
351  bool success = true;
352 
353  //const size_type
354  //n = nlp->n(),
355  //m = nlp->m();
356 
357  // //////////////////////////////////////////////
358  // Validate the input
359 
361  space_x = nlp->space_x(),
362  space_c = nlp->space_c();
363 
364  const CalcFiniteDiffProd
365  &fd_deriv_prod = this->calc_fd_prod();
366 
367  const value_type
368  rand_y_l = -1.0, rand_y_u = 1.0,
369  small_num = ::pow(std::numeric_limits<value_type>::epsilon(),0.25);
370 
371  if(out)
372  *out
373  << "\nComparing directional products Gf'*y and/or Gc'*y with finite difference values "
374  " FDGf'*y and/or FDGc'*y for random y's ...\n";
375 
376  value_type max_Gf_warning_viol = 0.0;
377  int num_Gf_warning_viol = 0;
378 
379  VectorSpace::vec_mut_ptr_t
380  y = space_x->create_member(),
381  Gc_prod = ( Gc ? space_c->create_member() : Teuchos::null ),
382  FDGc_prod = ( Gc ? space_c->create_member() : Teuchos::null );
383 
384  const int num_fd_directions_used = ( num_fd_directions() > 0 ? num_fd_directions() : 1 );
385 
386  for( int direc_i = 1; direc_i <= num_fd_directions_used; ++direc_i ) {
387  if( num_fd_directions() > 0 ) {
388  random_vector( rand_y_l, rand_y_u, y.get() );
389  if(out)
390  *out
391  << "\n****"
392  << "\n**** Random directional vector " << direc_i << " ( ||y||_1 / n = "
393  << (y->norm_1() / y->dim()) << " )"
394  << "\n***\n";
395  }
396  else {
397  *y = 1.0;
398  if(out)
399  *out
400  << "\n****"
401  << "\n**** Ones vector y ( ||y||_1 / n = "<<(y->norm_1()/y->dim())<<" )"
402  << "\n***\n";
403  }
404  // Compute exact??? values
405  value_type
406  Gf_y = Gf ? dot( *Gf, *y ) : 0.0;
407  if(Gc)
408  V_MtV( Gc_prod.get(), *Gc, BLAS_Cpp::trans, *y );
409  // Compute finite difference values
410  value_type
411  FDGf_y;
412  const bool preformed_fd = fd_deriv_prod.calc_deriv_product(
413  xo,xl,xu
414  ,*y
415  ,NULL // fo
416  ,NULL // co
417  ,true // check_nan_inf
418  ,nlp
419  ,Gf ? &FDGf_y : NULL
420  ,Gc ? FDGc_prod.get() : NULL
421  ,out
422  );
423  if( !preformed_fd ) {
424  if(out)
425  *out
426  << "\nError, the finite difference computation was not preformed due to cramped bounds\n"
427  << "Finite difference test failed!\n" << endl;
428  return false;
429  }
430 
431  // Compare the quantities
432  // Gf
433  assert_print_nan_inf(FDGf_y, "FDGf'*y",true,out);
434  const value_type
435  Gf_err = ::fabs( Gf_y - FDGf_y ) / ( ::fabs(Gf_y) + ::fabs(FDGf_y) + small_num );
436  if(out)
437  *out
438  << "\nrel_err(Gf'*y,FDGf'*y) = "
439  << "rel_err(" << Gf_y << "," << FDGf_y << ") = "
440  << Gf_err << endl;
441  if( Gf_err >= warning_tol() ) {
442  max_Gf_warning_viol = my_max( max_Gf_warning_viol, Gf_err );
443  ++num_Gf_warning_viol;
444  }
445  if( Gf_err >= error_tol() ) {
446  if(out)
447  *out
448  << "\nError, exceeded Gf_error_tol = " << error_tol() << endl
449  << "Stoping the tests!\n";
450  return false;
451  }
452  // Gc
453  if(Gc) {
454  const value_type
455  sum_Gc_prod = sum(*Gc_prod),
456  sum_FDGc_prod = sum(*FDGc_prod);
457  assert_print_nan_inf(sum_FDGc_prod, "sum(FDGc'*y)",true,out);
458  const value_type
459  calc_err = ::fabs( ( sum_Gc_prod - sum_FDGc_prod )
460  /( ::fabs(sum_Gc_prod) + ::fabs(sum_FDGc_prod) + small_num ) );
461  if(out)
462  *out
463  << "\nrel_err(sum(Gc'*y),sum(FDGc'*y)) = "
464  << "rel_err(" << sum_Gc_prod << "," << sum_FDGc_prod << ") = "
465  << calc_err << endl;
466  if( calc_err >= error_tol() ) {
467  if(out)
468  *out
469  << "\nError, rel_err(sum(Gc'*y),sum(FDGc'*y)) = "
470  << "rel_err(" << sum_Gc_prod << "," << sum_FDGc_prod << ") = "
471  << calc_err << endl
472  << "exceeded error_tol = " << error_tol() << endl
473  << "Stoping the tests!\n";
474  if(print_all_warnings)
475  *out << "\ny =\n" << *y
476  << "\nGc_prod =\n" << *Gc_prod
477  << "\nFDGc_prod =\n" << *FDGc_prod;
478  update_success( false, &success );
479  return false;
480  }
481  if( calc_err >= warning_tol() ) {
482  if(out)
483  *out
484  << "\nWarning, rel_err(sum(Gc'*y),sum(FDGc'*y)) = "
485  << "rel_err(" << sum_Gc_prod << "," << sum_FDGc_prod << ") = "
486  << calc_err << endl
487  << "exceeded warning_tol = " << warning_tol() << endl;
488  }
489  }
490  }
491  if(out && num_Gf_warning_viol)
492  *out
493  << "\nFor Gf, there were " << num_Gf_warning_viol << " warning tolerance "
494  << "violations out of num_fd_directions = " << num_fd_directions()
495  << " computations of FDGf'*y\n"
496  << "and the maximum violation was " << max_Gf_warning_viol
497  << " > Gf_waring_tol = " << warning_tol() << endl;
498 
499  if(out)
500  *out
501  << "\nCongradulations! All of the computed errors were within the specified error tolerance!\n";
502 
503  return true;
504 }
505 
506 } // end namespace NLPInterfacePack
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Strategy interface for computing the product of the derivatives of the functions of an NLP along give...
virtual size_type n() const
Return the number of variables.
size_t size_type
bool finite_diff_check(NLP *nlp, const Vector &xo, const Vector *xl, const Vector *xu, const MatrixOp *Gc, const Vector *Gf, bool print_all_warnings, std::ostream *out) const
This function takes an NLP object and its computed derivatives and function values and validates the ...
Teuchos::RCP< const VectorSpace > vec_space_ptr_t
NLP interface class {abstract}.
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
NLPFirstDerivTester(const calc_fd_prod_ptr_t &calc_fd_prod=Teuchos::rcp(new CalcFiniteDiffProd()), ETestingMethod fd_testing_method=FD_DIRECTIONAL, size_type num_fd_directions=1, value_type warning_tol=1e-8, value_type error_tol=1e-3)
Constructor.
virtual bool calc_deriv_product(const Vector &xo, const Vector *xl, const Vector *xu, const Vector &v, const value_type *fo, const Vector *co, bool check_nan_inf, NLP *nlp, value_type *Gf_prod, VectorMutable *Gc_prod, std::ostream *out, bool trace=false, bool dump_all=false) const
Compute the directional derivatives by finite differences.
value_type sum(const Vector &v_rhs)
virtual size_type m() const
Return the number of general equality constraints.
void random_vector(value_type l, value_type u, VectorMutable *v)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)