Sundance 2.1

00001 /* @HEADER@ */
00002 // ************************************************************************
00003 // 
00004 //                              Sundance
00005 //                 Copyright (2005) Sandia Corporation
00006 // 
00007 // Copyright (year first published) Sandia Corporation.  Under the terms 
00008 // of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 
00009 // retains certain rights in this software.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //                                                                                 
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA                                                                                
00025 // Questions? Contact Kevin Long (krlong@sandia.gov), 
00026 // Sandia National Laboratories, Livermore, California, USA
00027 // 
00028 // ************************************************************************
00029 /* @HEADER@ */
00030 
00031 #include "SundanceUnaryFunctor.hpp"
00032 #include "PlayaTabs.hpp"
00033 #include "SundanceOut.hpp"
00034 
00035 using namespace Sundance;
00036 using namespace Teuchos;
00037 
00038 void UnaryFunctor::eval1(const double* const x, 
00039   int nx, 
00040   double* f, 
00041   double* df_dx) const
00042 {
00043   evalFDDerivs1(x, nx, f, df_dx);
00044 }
00045 
00046 
00047 void UnaryFunctor::eval2(const double* const x, 
00048   int nx, 
00049   double* f, 
00050   double* df_dx,
00051   double* d2f_dxx) const
00052 {
00053   evalFDDerivs2(x, nx, f, df_dx, d2f_dxx);
00054 }
00055 
00056 void UnaryFunctor::eval3(const double* const x, 
00057   int nx, 
00058   double* f, 
00059   double* df_dx,
00060   double* d2f_dxx,
00061   double* d3f_dxxx) const
00062 {
00063   evalFDDerivs3(x, nx, f, df_dx, d2f_dxx, d3f_dxxx);
00064 }
00065 
00066 
00067 void UnaryFunctor::evalFDDerivs1(const double* const x, 
00068   int nx, 
00069   double* f, 
00070   double* df_dx) const
00071 {
00072   eval0(x, nx, f);
00073 
00074   double w1 = 1.0/h_/12.0;
00075 
00076   for (int i=0; i<nx; i++)
00077   {
00078     double xPlus1 = x[i] + h_;
00079     double xMinus1 = x[i] - h_;
00080     double fPlus1;
00081     double fMinus1;
00082     double xPlus2 = x[i] + 2.0*h_;
00083     double xMinus2 = x[i] - 2.0*h_;
00084     double fPlus2;
00085     double fMinus2;
00086     eval0(&xPlus1, 1, &fPlus1);
00087     eval0(&xMinus1, 1, &fMinus1);
00088     eval0(&xPlus2, 1, &fPlus2);
00089     eval0(&xMinus2, 1, &fMinus2);
00090     df_dx[i] = w1*( 8.0*(fPlus1-fMinus1) - (fPlus2 - fMinus2) );
00091   }
00092 }
00093 
00094 void UnaryFunctor::evalFDDerivs2(const double* const x, 
00095   int nx, 
00096   double* f, 
00097   double* df_dx,
00098   double* d2f_dxx) const
00099 {
00100   eval0(x, nx, f);
00101 
00102   double w1 = 1.0/h_/12.0;
00103   double w2 = w1/h_;
00104 
00105   for (int i=0; i<nx; i++)
00106   {
00107     double xPlus1 = x[i] + h_;
00108     double xMinus1 = x[i] - h_;
00109     double fPlus1;
00110     double fMinus1;
00111     double xPlus2 = x[i] + 2.0*h_;
00112     double xMinus2 = x[i] - 2.0*h_;
00113     double fPlus2;
00114     double fMinus2;
00115     eval0(&xPlus1, 1, &fPlus1);
00116     eval0(&xMinus1, 1, &fMinus1);
00117     eval0(&xPlus2, 1, &fPlus2);
00118     eval0(&xMinus2, 1, &fMinus2);
00119     df_dx[i] = w1*( 8.0*(fPlus1-fMinus1) - (fPlus2 - fMinus2) );
00120     d2f_dxx[i] = w2*(16.0*(fPlus1 + fMinus1) 
00121       - 30.0*f[i] - (fPlus2+fMinus2));
00122   }
00123 }
00124 
00125 void UnaryFunctor::evalFDDerivs3(const double* const x, 
00126   int nx, 
00127   double* f, 
00128   double* df_dx,
00129   double* d2f_dxx,
00130   double* d3f_dxxx) const
00131 {
00132   eval0(x, nx, f);
00133 
00134   double w1 = 1.0/h_/12.0;
00135   double w2 = w1/h_;
00136   double w3 = 0.5/h_/h_/h_;
00137 
00138   for (int i=0; i<nx; i++)
00139   {
00140     double xPlus1 = x[i] + h_;
00141     double xMinus1 = x[i] - h_;
00142     double fPlus1;
00143     double fMinus1;
00144     double xPlus2 = x[i] + 2.0*h_;
00145     double xMinus2 = x[i] - 2.0*h_;
00146     double fPlus2;
00147     double fMinus2;
00148     eval0(&xPlus1, 1, &fPlus1);
00149     eval0(&xMinus1, 1, &fMinus1);
00150     eval0(&xPlus2, 1, &fPlus2);
00151     eval0(&xMinus2, 1, &fMinus2);
00152     df_dx[i] = w1*( 8.0*(fPlus1-fMinus1) - (fPlus2 - fMinus2) );
00153     d2f_dxx[i] = w2*(16.0*(fPlus1 + fMinus1) 
00154       - 30.0*f[i] - (fPlus2+fMinus2));
00155     d3f_dxxx[i] = w3*(fPlus2 - 2.0*fPlus1 + 2.0*fMinus1 - fMinus2);
00156   }
00157 }
00158 
00159 bool UnaryFunctor::testDerivs(const double& x, const double& tol) const
00160 {
00161   Tabs tabs;
00162   Out::os() << tabs << std::endl << tabs 
00163             << "comparing exact derivs to FD derivs for functor " 
00164             << name() << std::endl;
00165 
00166   double fExact;
00167   double fxExact;
00168   double fxxExact;
00169 
00170   double fFD;
00171   double fxFD;
00172   double fxxFD;
00173 
00174   bool isOK = true;
00175 
00176   /* test first differentiation */
00177   Out::os() << tabs << "computing first derivatives at x=" << x << std::endl;
00178   {
00179     Tabs tabs1;
00180 
00181     eval1(&x, 1, &fExact, &fxExact);
00182     Out::os() << tabs1 << "Exact: f=" << fExact << " df_dx=" << fxExact << std::endl; 
00183     evalFDDerivs1(&x, 1, &fFD, &fxFD);
00184     Out::os() << tabs1 << "FD:    f=" << fFD << " df_dx=" << fxFD << std::endl; 
00185 
00186     double fError = fabs(fFD - fExact)/(fabs(fExact) + h_);
00187     double fxError = fabs(fxFD - fxExact)/(fabs(fxExact)+h_);
00188     {
00189       Tabs tabs2;
00190       Out::os() << tabs2 << "| f-f_FD |=" << fError 
00191                 << " | df-df_FD |=" << fxError << std::endl;
00192       if (fError > tol) 
00193       {
00194         Out::os() << tabs << "ERROR: function value mismatch!" << std::endl;
00195         isOK = false;
00196       }
00197       if (fxError > tol) 
00198       {
00199         Out::os() << tabs << "ERROR: first derivative mismatch!" << std::endl;
00200         isOK = false;
00201       }
00202     }
00203   }
00204   
00205   Out::os() << tabs << std::endl;
00206 
00207   
00208   /* test first and second differentiation */
00209   Out::os() << tabs << "computing first and second derivatives at x=" << x << std::endl;
00210   {
00211     Tabs tabs1;
00212 
00213     eval2(&x, 1, &fExact, &fxExact, &fxxExact);
00214     Out::os() << tabs1 << "Exact: f=" << fExact 
00215               << " df_dx=" << fxExact 
00216               << " d2f_dx2=" << fxxExact 
00217               << std::endl; 
00218 
00219     evalFDDerivs2(&x, 1, &fFD, &fxFD, &fxxFD);
00220     Out::os() << tabs1 << "FD:    f=" << fFD << " df_dx=" << fxFD  
00221               << " d2f_dx2=" << fxxFD << std::endl; 
00222 
00223     double fError = fabs(fFD - fExact)/(fabs(fExact) + h_);
00224     double fxError = fabs(fxFD - fxExact)/(fabs(fxExact)+h_);
00225     double fxxError = fabs(fxxFD - fxxExact)/(fabs(fxxExact)+h_);
00226     {
00227       Tabs tabs2;
00228       Out::os() << tabs2 << "| f-f_FD |=" << fError 
00229                 << " | df-df_FD |=" << fxError 
00230                 << " | d2f-d2f_FD |=" << fxxError << std::endl;
00231       if (fError > tol) 
00232       {
00233         Out::os() << tabs << "ERROR: function value mismatch!" << std::endl;
00234         isOK = false;
00235       }
00236       if (fxError > tol) 
00237       {
00238         Out::os() << tabs << "ERROR: first derivative mismatch!" << std::endl;
00239         isOK = false;
00240       }
00241       if (fxxError > tol) 
00242       {
00243         Out::os() << tabs << "ERROR: second derivative mismatch!" << std::endl;
00244         isOK = false;
00245       }
00246     }
00247   }
00248 
00249   Out::os() << tabs << std::endl;
00250   return isOK;
00251 }
00252 
00253 bool UnaryFunctor::testInvalidValue(const double& badValue) const 
00254 {
00255   Tabs tabs;
00256   bool detectedError = false;
00257 
00258   Out::os() << std::endl << tabs 
00259             << "testing exception detection for bad value x=" << badValue
00260             << " for function  " << name() << std::endl;
00261   try
00262   {
00263     testDerivs(badValue, 1.0);
00264   }
00265   catch(std::exception& e)
00266   {
00267     detectedError = true;
00268   }
00269 
00270   if (!detectedError) 
00271   {
00272     Out::os() << tabs << "ERROR: missed detection of an exception at x=" << badValue
00273               << " for function " << name() << std::endl;
00274   }
00275   else
00276   {
00277     Out::os() << tabs << "the error was detected!" << std::endl;
00278   }
00279   return detectedError;
00280 }
00281 
00282 
00283 bool UnaryFunctor::test(int nx, const double& tol) const
00284 {
00285   bool isOK = true;
00286 
00287   double a = -sqrt(1.5);
00288   double b = sqrt(2.0);
00289 
00290   /* first test a sample of points in the domain */
00291   if (domain()->hasLowerBound())
00292   {
00293     a = domain()->lowerBound();
00294   }
00295   if (domain()->hasUpperBound())
00296   {
00297     b = domain()->upperBound();
00298   }
00299 
00300   double c = (b-a)/((double) nx+1);
00301   for (int i=1; i<=nx; i++)
00302   {
00303     double x = a + i*c;
00304     Out::os() << "testing at " << x << std::endl;
00305     if (domain()->hasExcludedPoint() && fabs(x-domain()->excludedPoint())<1.0e-14)
00306     {
00307       continue;
00308     }
00309     isOK =  testDerivs(x, tol) && isOK ;
00310   }
00311 
00312   /* Test for exception detection at the excluded point, if any */
00313   if (domain()->hasExcludedPoint())
00314   {
00315     Out::os() << "testing detection of excluded point x=" 
00316               << domain()->excludedPoint() << std::endl;
00317     isOK =  testInvalidValue(domain()->excludedPoint()) && isOK ;
00318   }
00319 
00320   /* Test for exception detection at a point below the lower bound */
00321   if (domain()->hasLowerBound())
00322   {
00323     Out::os() << "testing exception below lower bound x=" 
00324               << domain()->lowerBound() << std::endl;
00325     isOK =  testInvalidValue(domain()->lowerBound() - 0.1) && isOK ;
00326   }
00327   
00328   /* Test for exception detection at a point above the upper bound */
00329   if (domain()->hasUpperBound())
00330   {
00331     Out::os() << "testing exception above lower bound x=" 
00332               << domain()->upperBound() << std::endl;
00333 
00334     isOK =  testInvalidValue(domain()->upperBound() + 0.1) && isOK ;
00335   }
00336   
00337   return isOK;
00338   
00339 
00340 }
SundanceUnaryFunctor.cpp
