00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
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
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
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
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
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
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
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 }