libMesh
autodiff.C
Go to the documentation of this file.
1 #include "libmesh/fparser_ad.hh"
2 
3 #include "libmesh_cppunit.h"
4 
5 class FParserAutodiffTest : public CppUnit::TestCase
6 {
7 public:
8  CPPUNIT_TEST_SUITE ( FParserAutodiffTest );
9 
10  CPPUNIT_TEST ( runTests );
11  CPPUNIT_TEST ( registerDerivativeTest );
12  CPPUNIT_TEST ( registerDerivativeRepeatTest );
13 
14  CPPUNIT_TEST_SUITE_END ();
15 
16 private:
17 
18  class ADTest {
19  public:
20  ADTest(const std::string & _func,
21  double _min, double _max, double _dx = 1e-6,
22  double _reltol = 1e-5, int _steps = 20, double _abstol = 1e-10) :
23  func(_func),
24  min(_min),
25  max(_max),
26  dx(_dx),
27  reltol(_reltol),
28  abstol(_abstol),
29  steps(_steps)
30  {
31  CPPUNIT_ASSERT_MESSAGE ("Failed to parse test function", F.Parse(func, "x") == -1);
32  dF.Parse(func, "x");
33  dFopt.Parse(func, "x");
34  dFaopt.Parse(func, "x");
35 
36  CPPUNIT_ASSERT_MESSAGE ("Failed to take derivative of function", dF.AutoDiff("x") == -1);
37 
38  dFopt.Optimize();
39  CPPUNIT_ASSERT_MESSAGE ("Failed to take derivative of optimized function", dFopt.AutoDiff("x") == -1);
40 
41  dFaopt.SetADFlags(FunctionParserAD::ADAutoOptimize, true);
42  CPPUNIT_ASSERT_MESSAGE ("Failed to take derivative of auto-optimized function", dFaopt.AutoDiff("x") == -1);
43  }
44 
45  bool run()
46  {
47  double x1, x2, vdF, vF1, vF2, fd;
48  for (double x = min; x <= max; x += (max-min) / double(steps))
49  {
50  x1 = x - dx/2.0;
51  x2 = x + dx/2.0;
52 
53  vF1 = F.Eval(&x1);
54  vF2 = F.Eval(&x2);
55  fd = (vF2-vF1) / dx;
56 
57  // CPPUNIT_ASSERT(std::abs(fd - vdF) > tol)
58  // CPPUNIT_ASSERT(std::abs(fd - vdFopt) > tol)
59  vdF = dF.Eval(&x);
60  if (std::abs(vdF) > abstol && std::abs((fd - vdF)/vdF) > reltol && std::abs(fd - vdF)> abstol)
61  {
62  std::cout << "Error in " << func << ": " << fd << "!=" << vdF << " at x=" << x << '\n';
63  return false;
64  }
65 
66  vdF = dFopt.Eval(&x);
67  if (std::abs(vdF) > abstol && std::abs((fd - vdF)/vdF) > reltol && std::abs(fd - vdF)> abstol)
68  {
69  std::cout << "Error in opt " << func << ": " << fd << "!=" << vdF << " at x=" << x << '\n';
70  return false;
71  }
72 
73  vdF = dFaopt.Eval(&x);
74  if (std::abs(vdF) > abstol && std::abs((fd - vdF)/vdF) > reltol && std::abs(fd - vdF)> abstol)
75  {
76  std::cout << "Error in auto opt " << func << ": " << fd << "!=" << vdF << " at x=" << x << '\n';
77  return false;
78  }
79  }
80 
81  return true;
82  }
83 
84  private:
85  std::string func;
86  double min, max, dx, reltol, abstol;
87  int steps;
88 
89  FunctionParserAD F, dF, dFopt, dFaopt;
90  };
91 
92  std::vector<ADTest> tests;
93 
94 public:
95  virtual void setUp()
96  {
97  tests.push_back(ADTest("log(x*x) + log2(2*x) + log10(4*x)", 0.1, 3.0));
98  tests.push_back(ADTest("sin(-x) + cos(2*x) + tan(4*x)", -5.0, 5.0, 1e-7, 1e-5, 100));
99  tests.push_back(ADTest("sinh(-x) + cosh(x/2) + tanh(x/3)", -4.0, 4.0, 0.0001, 1e-5, 100));
100  tests.push_back(ADTest("plog(-x,0.01)", 0.001, 0.05, 0.00001, 1e-5, 100));
101  tests.push_back(ADTest("2 + 4*x + 8*x^2 + 16*x^3 + 32*x^4", -5.0, 5.0, 1e-5,1e-4));
102  tests.push_back(ADTest("1/x^2", 0.01, 2.0, 1e-8));
103  tests.push_back(ADTest("sqrt(x*2)", 0.001, 2.0, 1e-6));
104  tests.push_back(ADTest("abs(x*2)", -1.99, 2.0));
105  tests.push_back(ADTest("asin(-x)", -0.99, 0.99));
106  tests.push_back(ADTest("acos(-x)", -0.99, 0.99));
107  tests.push_back(ADTest("atan(-x)", -99, 99));
108  tests.push_back(ADTest("x*sin(-x)*log(x)*tanh(x)", 0.001, 5, 1e-8));
109  tests.push_back(ADTest("exp(-x) + 2*exp2(x)", -1.0, 2.0));
110  tests.push_back(ADTest("hypot(2*x,1) - hypot(1,4*x)", -10, 10.0));
111  tests.push_back(ADTest("if(x<0, (-x)^3, x^3)", -1.0, 1.0));
112  tests.push_back(ADTest("max(x^2-0.5,0)", -1.5, 1.5));
113  tests.push_back(ADTest("min(x^2-0.5,0)", -1.5, 1.5));
114  tests.push_back(ADTest("atan2(x,1) + atan2(2,x)", -0.99, 0.99));
115  tests.push_back(ADTest("0.767^sin(x)", -1.5, 1.5));
116  tests.push_back(ADTest("A := sin(x) + tanh(x); A + sqrt(A) - x", -1.5, 1.5));
117  tests.push_back(ADTest("3*sin(2*x)*sin(2*x)", -5.0, 5.0, 1e-7, 1e-5, 100));
118  tests.push_back(ADTest("erf(0.5*x)", -2., 2., 1e-6, 1e-5, 100));
119  }
120 
121  void runTests()
122  {
123  const unsigned int ntests = tests.size();
124 
125  unsigned int passed = 0;
126  for (unsigned i = 0; i < ntests; ++i)
127  passed += tests[i].run() ? 1 : 0;
128 
129  CPPUNIT_ASSERT_EQUAL (passed, ntests);
130  }
131 
132  void registerDerivativeTest()
133  {
134  FunctionParserAD R;
135  std::string func = "x*a";
136  R.SetADFlags(FunctionParserAD::ADCacheDerivatives, true);
137 
138  // Parse the input expression into bytecode
139  R.Parse(func, "x,a");
140 
141  // add a new variable y and map it to the da/dx derivative
142  R.AddVariable("y");
143  R.RegisterDerivative("a", "x", "y");
144 
145  // parameter vector
146  double p[3];
147  double & x = p[0];
148  double & a = p[1];
149  double & y = p[2];
150 
151  FunctionParserAD dR(R);
152  CPPUNIT_ASSERT_EQUAL (dR.AutoDiff("x"), -1);
153  dR.Optimize();
154  // dR = a+x*y
155 
156  FunctionParserAD d2R(dR);
157  CPPUNIT_ASSERT_EQUAL (d2R.AutoDiff("x"), -1);
158  d2R.Optimize();
159  // d2R = 2*y
160 
161  // we probe the parsers and check if they agree with the reference solution
162  for (x = -1.0; x < 1.0; x+=0.3726)
163  for (a = -1.0; a < 1.0; a+=0.2642)
164  for (y = -1.0; y < 1.0; y+=0.3156)
165  {
166  LIBMESH_ASSERT_FP_EQUAL(R.Eval(p), x*a, 1.e-12);
167  LIBMESH_ASSERT_FP_EQUAL(dR.Eval(p), a+x*y, 1.e-12);
168  LIBMESH_ASSERT_FP_EQUAL(d2R.Eval(p), 2*y, 1.e-12);
169  }
170  }
171 
172  void registerDerivativeRepeatTest()
173  {
174  // now do the same functional form but with a different mapping to see if the cache
175  // signature was correctly updated
176  FunctionParserAD R;
177  std::string func = "x*a";
178  R.SetADFlags(FunctionParserAD::ADCacheDerivatives, true);
179 
180  // Parse the input expression into bytecode
181  R.Parse(func, "x,a");
182 
183  // add a new variable y but do not map it to the da/dx derivative!
184  R.AddVariable("y");
185  R.RegisterDerivative("a", "x", "a");
186 
187  // parameter vector
188  double p[3];
189  double & x = p[0];
190  double & a = p[1];
191  double & y = p[2];
192 
193  FunctionParserAD dR(R);
194  CPPUNIT_ASSERT_EQUAL (dR.AutoDiff("x"), -1);
195  dR.Optimize();
196  // dR = a + x*a
197 
198  FunctionParserAD d2R(dR);
199  CPPUNIT_ASSERT_EQUAL (d2R.AutoDiff("x"), -1);
200  d2R.Optimize();
201  // d2R = 2*a + x*a
202 
203  // we probe the parsers and check if they agree with the reference solution
204  for (x = -1.0; x < 1.0; x+=0.3726)
205  for (a = -1.0; a < 1.0; a+=0.2642)
206  for (y = -1.0; y < 1.0; y+=0.3156)
207  {
208  LIBMESH_ASSERT_FP_EQUAL(R.Eval(p), x*a, 1.e-12);
209  LIBMESH_ASSERT_FP_EQUAL(dR.Eval(p), a+x*a, 1.e-12);
210  LIBMESH_ASSERT_FP_EQUAL(d2R.Eval(p), 2*a+x*a, 1.e-12);
211  }
212  }
213 };
214 
215 CPPUNIT_TEST_SUITE_REGISTRATION ( FParserAutodiffTest );
FParserAutodiffTest::ADTest::reltol
double reltol
Definition: autodiff.C:86
FParserAutodiffTest::CPPUNIT_TEST_SUITE_END
CPPUNIT_TEST_SUITE_END()
FParserAutodiffTest::ADTest::steps
int steps
Definition: autodiff.C:87
FParserAutodiffTest::ADTest
Definition: autodiff.C:18
FParserAutodiffTest::registerDerivativeRepeatTest
void registerDerivativeRepeatTest()
Definition: autodiff.C:172
FParserAutodiffTest::setUp
virtual void setUp()
Definition: autodiff.C:95
FParserAutodiffTest::ADTest::dF
FunctionParserAD dF
Definition: autodiff.C:89
FParserAutodiffTest::registerDerivativeTest
void registerDerivativeTest()
Definition: autodiff.C:132
FParserAutodiffTest::CPPUNIT_TEST_SUITE
CPPUNIT_TEST_SUITE(FParserAutodiffTest)
std::abs
MetaPhysicL::DualNumber< T, D > abs(const MetaPhysicL::DualNumber< T, D > &in)
CPPUNIT_TEST_SUITE_REGISTRATION
CPPUNIT_TEST_SUITE_REGISTRATION(FParserAutodiffTest)
FParserAutodiffTest::ADTest::ADTest
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:20
FParserAutodiffTest
Definition: autodiff.C:5
FParserAutodiffTest::ADTest::func
std::string func
Definition: autodiff.C:85
FParserAutodiffTest::ADTest::min
double min
Definition: autodiff.C:86
libmesh_cppunit.h
FParserAutodiffTest::runTests
void runTests()
Definition: autodiff.C:121
FParserAutodiffTest::ADTest::dFopt
FunctionParserAD dFopt
Definition: autodiff.C:89
FParserAutodiffTest::ADTest::run
bool run()
Definition: autodiff.C:45
FParserAutodiffTest::ADTest::dFaopt
FunctionParserAD dFaopt
Definition: autodiff.C:89
FParserAutodiffTest::ADTest::F
FunctionParserAD F
Definition: autodiff.C:89
FParserAutodiffTest::CPPUNIT_TEST
CPPUNIT_TEST(runTests)
FParserAutodiffTest::ADTest::max
double max
Definition: autodiff.C:86
FParserAutodiffTest::ADTest::dx
double dx
Definition: autodiff.C:86
FParserAutodiffTest::tests
std::vector< ADTest > tests
Definition: autodiff.C:92
FParserAutodiffTest::ADTest::abstol
double abstol
Definition: autodiff.C:86
Generated on Sat Jan 25 2020 12:07:06 for libMesh by   doxygen 1.8.16