libMesh
autodiff.C
Go to the documentation of this file.
1 #include "libmesh/fparser_ad.hh"
2 
3 // Ignore unused parameter warnings coming from cppunit headers
4 #include <libmesh/ignore_warnings.h>
5 #include <cppunit/extensions/HelperMacros.h>
6 #include <cppunit/TestCase.h>
7 #include <libmesh/restore_warnings.h>
8 
9 // THE CPPUNIT_TEST_SUITE_END macro expands to code that involves
10 // std::auto_ptr, which in turn produces -Wdeprecated-declarations
11 // warnings. These can be ignored in GCC as long as we wrap the
12 // offending code in appropriate pragmas. We can't get away with a
13 // single ignore_warnings.h inclusion at the beginning of this file,
14 // since the libmesh headers pull in a restore_warnings.h at some
15 // point. We also don't bother restoring warnings at the end of this
16 // file since it's not a header.
17 #include <libmesh/ignore_warnings.h>
18 
19 class FParserAutodiffTest : public CppUnit::TestCase
20 {
21 public:
23 
26 
28 
29 private:
30 
31  class ADTest {
32  public:
33  ADTest(const std::string & _func,
34  double _min, double _max, double _dx = 1e-6,
35  double _reltol = 1e-5, int _steps = 20, double _abstol = 1e-10) :
36  func(_func),
37  min(_min),
38  max(_max),
39  dx(_dx),
40  reltol(_reltol),
41  abstol(_abstol),
42  steps(_steps)
43  {
44  CPPUNIT_ASSERT_MESSAGE ("Failed to parse test function", F.Parse(func, "x") == -1);
45  dF.Parse(func, "x");
46  dFopt.Parse(func, "x");
47  dFaopt.Parse(func, "x");
48 
49  CPPUNIT_ASSERT_MESSAGE ("Failed to take derivative of function", dF.AutoDiff("x") == -1);
50 
51  dFopt.Optimize();
52  CPPUNIT_ASSERT_MESSAGE ("Failed to take derivative of optimized function", dFopt.AutoDiff("x") == -1);
53 
54  dFaopt.SetADFlags(FunctionParserAD::ADAutoOptimize, true);
55  CPPUNIT_ASSERT_MESSAGE ("Failed to take derivative of auto-optimized function", dFaopt.AutoDiff("x") == -1);
56  }
57 
58  bool run()
59  {
60  double x1, x2, vdF, vF1, vF2, fd;
61  for (double x = min; x <= max; x += (max-min) / double(steps))
62  {
63  x1 = x - dx/2.0;
64  x2 = x + dx/2.0;
65 
66  vF1 = F.Eval(&x1);
67  vF2 = F.Eval(&x2);
68  fd = (vF2-vF1) / dx;
69 
70  // CPPUNIT_ASSERT(std::abs(fd - vdF) > tol)
71  // CPPUNIT_ASSERT(std::abs(fd - vdFopt) > tol)
72  vdF = dF.Eval(&x);
73  if (std::abs(vdF) > abstol && std::abs((fd - vdF)/vdF) > reltol && std::abs(fd - vdF)> abstol)
74  {
75  std::cout << "Error in " << func << ": " << fd << "!=" << vdF << " at x=" << x << '\n';
76  return false;
77  }
78 
79  vdF = dFopt.Eval(&x);
80  if (std::abs(vdF) > abstol && std::abs((fd - vdF)/vdF) > reltol && std::abs(fd - vdF)> abstol)
81  {
82  std::cout << "Error in opt " << func << ": " << fd << "!=" << vdF << " at x=" << x << '\n';
83  return false;
84  }
85 
86  vdF = dFaopt.Eval(&x);
87  if (std::abs(vdF) > abstol && std::abs((fd - vdF)/vdF) > reltol && std::abs(fd - vdF)> abstol)
88  {
89  std::cout << "Error in auto opt " << func << ": " << fd << "!=" << vdF << " at x=" << x << '\n';
90  return false;
91  }
92  }
93 
94  return true;
95  }
96 
97  private:
98  std::string func;
99  double min, max, dx, reltol, abstol;
100  int steps;
101 
102  FunctionParserAD F, dF, dFopt, dFaopt;
103  };
104 
105  std::vector<ADTest> tests;
106 
107 public:
108  virtual void setUp()
109  {
110  tests.push_back(ADTest("log(x*x) + log2(2*x) + log10(4*x)", 0.1, 3.0));
111  tests.push_back(ADTest("sin(-x) + cos(2*x) + tan(4*x)", -5.0, 5.0, 1e-7, 1e-5, 100));
112  tests.push_back(ADTest("sinh(-x) + cosh(x/2) + tanh(x/3)", -4.0, 4.0, 0.0001, 1e-5, 100));
113  tests.push_back(ADTest("plog(-x,0.01)", 0.001, 0.05, 0.00001, 1e-5, 100));
114  tests.push_back(ADTest("2 + 4*x + 8*x^2 + 16*x^3 + 32*x^4", -5.0, 5.0, 1e-5,1e-4));
115  tests.push_back(ADTest("1/x^2", 0.01, 2.0, 1e-8));
116  tests.push_back(ADTest("sqrt(x*2)", 0.001, 2.0, 1e-6));
117  tests.push_back(ADTest("abs(x*2)", -1.99, 2.0));
118  tests.push_back(ADTest("asin(-x)", -0.99, 0.99));
119  tests.push_back(ADTest("acos(-x)", -0.99, 0.99));
120  tests.push_back(ADTest("atan(-x)", -99, 99));
121  tests.push_back(ADTest("x*sin(-x)*log(x)*tanh(x)", 0.001, 5, 1e-8));
122  tests.push_back(ADTest("exp(-x) + 2*exp2(x)", -1.0, 2.0));
123  tests.push_back(ADTest("hypot(2*x,1) - hypot(1,4*x)", -10, 10.0));
124  tests.push_back(ADTest("if(x<0, (-x)^3, x^3)", -1.0, 1.0));
125  tests.push_back(ADTest("max(x^2-0.5,0)", -1.5, 1.5));
126  tests.push_back(ADTest("min(x^2-0.5,0)", -1.5, 1.5));
127  tests.push_back(ADTest("atan2(x,1) + atan2(2,x)", -0.99, 0.99));
128  tests.push_back(ADTest("0.767^sin(x)", -1.5, 1.5));
129  tests.push_back(ADTest("A := sin(x) + tanh(x); A + sqrt(A) - x", -1.5, 1.5));
130  tests.push_back(ADTest("3*sin(2*x)*sin(2*x)", -5.0, 5.0, 1e-7, 1e-5, 100));
131  tests.push_back(ADTest("erf(0.5*x)", -2., 2., 1e-6, 1e-5, 100));
132  }
133 
134  void runTests()
135  {
136  const unsigned int ntests = tests.size();
137 
138  unsigned int passed = 0;
139  for (unsigned i = 0; i < ntests; ++i)
140  passed += tests[i].run() ? 1 : 0;
141 
142  CPPUNIT_ASSERT_EQUAL (passed, ntests);
143  }
144 
146  {
147  FunctionParserAD R;
148  std::string func = "x*a";
149 
150  // Parse the input expression into bytecode
151  R.Parse(func, "x,a");
152 
153  // add a new variable y and map it to the da/dx derivative
154  R.AddVariable("y");
155  R.RegisterDerivative("a", "x", "y");
156 
157  // parameter vector
158  double p[3];
159  double & x = p[0];
160  double & a = p[1];
161  double & y = p[2];
162 
163  FunctionParserAD dR(R);
164  CPPUNIT_ASSERT_EQUAL (dR.AutoDiff("x"), -1);
165  dR.Optimize();
166  // dR = a+x*y
167 
168  FunctionParserAD d2R(dR);
169  CPPUNIT_ASSERT_EQUAL (d2R.AutoDiff("x"), -1);
170  d2R.Optimize();
171  // d2R = 2*y
172 
173  // we probe the parsers and check if they agree with the reference solution
174  for (x = -1.0; x < 1.0; x+=0.3726)
175  for (a = -1.0; a < 1.0; a+=0.2642)
176  for (y = -1.0; y < 1.0; y+=0.3156)
177  {
178  CPPUNIT_ASSERT_DOUBLES_EQUAL(R.Eval(p), x*a, 1.e-12);
179  CPPUNIT_ASSERT_DOUBLES_EQUAL(dR.Eval(p), a+x*y, 1.e-12);
180  CPPUNIT_ASSERT_DOUBLES_EQUAL(d2R.Eval(p), 2*y, 1.e-12);
181  }
182  }
183 };
184 
FunctionParserAD F
Definition: autodiff.C:102
double abs(double a)
std::vector< ADTest > tests
Definition: autodiff.C:105
virtual void setUp()
Definition: autodiff.C:108
void registerDerivativeTest()
Definition: autodiff.C:145
PetscErrorCode Vec x
FunctionParserAD dF
Definition: autodiff.C:102
CPPUNIT_TEST(runTests)
CPPUNIT_TEST_SUITE_REGISTRATION(FParserAutodiffTest)
FunctionParserAD dFopt
Definition: autodiff.C:102
ADTest(const std::string &_func, double _min, double _max, double _dx=1e-6, double _reltol=1e-5, int _steps=20, double _abstol=1e-10)
Definition: autodiff.C:33
CPPUNIT_TEST_SUITE(FParserAutodiffTest)
FunctionParserAD dFaopt
Definition: autodiff.C:102