00001 #include "PlayaObjectiveBase.hpp"
00002 #include "PlayaOut.hpp"
00003 #include "PlayaTabs.hpp"
00004 #include "PlayaGhostView.hpp"
00005 #include "PlayaLinearCombinationImpl.hpp"
00006
00007
00008
00009 namespace Playa
00010 {
00011 using std::endl;
00012 using std::ostream;
00013 using std::setw;
00014
00015 bool ObjectiveBase::fdCheck(const Vector<double>& xIn,
00016 double tol,
00017 int verbosity) const
00018 {
00019 double f0, fPlus, fMinus;
00020
00021 double h = fdStep();
00022
00023 PLAYA_MSG1(verbosity, "running fdCheck with stepsize=" << h);
00024
00025 Vector<double> x = xIn.copy();
00026 Vector<double> x0 = x.copy();
00027 Vector<double> gf = x.copy();
00028 evalGrad(xIn, f0, gf);
00029
00030 int nTot = x.space().dim();
00031 int n = x.space().numLocalElements();
00032 int lowestIndex = x.space().baseGlobalNaturalIndex();
00033
00034 Array<double> df_dx(n);
00035
00036 for (int globalIndex=0; globalIndex<nTot; globalIndex++)
00037 {
00038 double tmp=0.0;
00039 bool isLocal = globalIndex >= lowestIndex
00040 && globalIndex < (lowestIndex+n);
00041 int localIndex = globalIndex - lowestIndex;
00042
00043 if (isLocal)
00044 {
00045 tmp = x[localIndex];
00046 x[localIndex] = tmp + h;
00047 }
00048
00049
00050 eval(x, fPlus);
00051
00052
00053 if (isLocal)
00054 {
00055 x[localIndex] = tmp - h;
00056 }
00057
00058
00059 eval(x, fMinus);
00060
00061
00062 if (isLocal)
00063 {
00064 df_dx[localIndex] = (fPlus - fMinus)/2.0/h;
00065 PLAYA_MSG2(verbosity,
00066 "g=" << setw(5) << globalIndex << ", l=" << setw(5) << localIndex
00067 << " f0=" << setw(12) << f0
00068 << " fPlus=" << setw(12) << fPlus
00069 << " fMinus=" << setw(12) << fMinus << " df_dx="
00070 << setw(12) << df_dx[localIndex]);
00071 PLAYA_MSG3(verbosity, "i " << globalIndex << " x_i=" << tmp
00072 << " f(x)=" << f0
00073 << " f(x+h)=" << fPlus
00074 << " f(x-h)=" << fMinus
00075 );
00076
00077 x[localIndex] = tmp;
00078 }
00079 }
00080
00081 double localMaxErr = 0.0;
00082
00083 VectorSpace<double> space = x.space();
00084 for (int i=0; i<space.numLocalElements(); i++)
00085 {
00086 double num = fabs(df_dx[i]-gf[i]);
00087 double den = fabs(df_dx[i]) + fabs(gf[i]) + 1.0e-14;
00088 double r = 0.0;
00089 if (fabs(den) > 1.0e-16) r = num/den;
00090 else r = 1.0;
00091 PLAYA_MSG2(verbosity, "i=" << setw(5) << i
00092 << " FD=" << setw(16) << df_dx[i]
00093 << " grad=" << setw(16) << gf[i]
00094 << " r=" << setw(16) << r);
00095 if (localMaxErr < r) localMaxErr = r;
00096 }
00097 PLAYA_MSG2(verbosity, "local max error = " << localMaxErr);
00098 double maxErr = localMaxErr;
00099 PLAYA_MSG1(verbosity, "fd check: max error = " << maxErr);
00100 return maxErr <= tol;
00101 }
00102
00103 }