PlayaObjectiveBase.cpp

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   /* get a FD step appropriate to this problem */
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     /* set point to xPlus */
00043     if (isLocal)
00044     {
00045       tmp = x[localIndex];
00046       x[localIndex] = tmp + h;
00047     }
00048 
00050     eval(x, fPlus);
00051 
00052     /* set point to xMinus */
00053     if (isLocal)
00054     {
00055       x[localIndex] = tmp - h;
00056     }
00057 
00058     /* eval at xMinus */
00059     eval(x, fMinus);
00060       
00061     /* check error */
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       /* restore point */
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 }