1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2  * -------------------------------------------------------------------------- *
3  *                                   Lepton                                   *
4  * -------------------------------------------------------------------------- *
5  * This is part of the Lepton expression parser originating from              *
6  * Simbios, the NIH National Center for Physics-Based Simulation of           *
7  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
8  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
9  *                                                                            *
10  * Portions copyright (c) 2013-2016 Stanford University and the Authors.      *
11  * Authors: Peter Eastman                                                     *
12  * Contributors:                                                              *
13  *                                                                            *
14  * Permission is hereby granted, free of charge, to any person obtaining a    *
15  * copy of this software and associated documentation files (the "Software"), *
16  * to deal in the Software without restriction, including without limitation  *
17  * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
18  * and/or sell copies of the Software, and to permit persons to whom the      *
19  * Software is furnished to do so, subject to the following conditions:       *
20  *                                                                            *
21  * The above copyright notice and this permission notice shall be included in *
22  * all copies or substantial portions of the Software.                        *
23  *                                                                            *
24  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
25  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
26  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
27  * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
28  * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
29  * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
30  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
31  * -------------------------------------------------------------------------- *
32 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
33 
34 /* -------------------------------------------------------------------------- *
35  *                                   lepton                                   *
36  * -------------------------------------------------------------------------- *
37  * This is part of the lepton expression parser originating from              *
38  * Simbios, the NIH National Center for Physics-Based Simulation of           *
39  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
40  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
41  *                                                                            *
42  * Portions copyright (c) 2009-2015 Stanford University and the Authors.      *
43  * Authors: Peter Eastman                                                     *
44  * Contributors:                                                              *
45  *                                                                            *
46  * Permission is hereby granted, free of charge, to any person obtaining a    *
47  * copy of this software and associated documentation files (the "Software"), *
48  * to deal in the Software without restriction, including without limitation  *
49  * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
50  * and/or sell copies of the Software, and to permit persons to whom the      *
51  * Software is furnished to do so, subject to the following conditions:       *
52  *                                                                            *
53  * The above copyright notice and this permission notice shall be included in *
54  * all copies or substantial portions of the Software.                        *
55  *                                                                            *
56  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
57  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
58  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
59  * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
60  * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
61  * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
62  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
63  * -------------------------------------------------------------------------- */
64 
65 #include "Operation.h"
66 #include "ExpressionTreeNode.h"
67 #include "MSVC_erfc.h"
68 
69 #ifndef M_PI
70 #define M_PI           3.14159265358979323846
71 #endif
72 
73 namespace PLMD {
74 using namespace lepton;
75 using namespace std;
76 
evaluate(double * args,const map<string,double> & variables) const77 double Operation::Erf::evaluate(double* args, const map<string, double>& variables) const {
78     return erf(args[0]);
79 }
80 
evaluate(double * args,const map<string,double> & variables) const81 double Operation::Erfc::evaluate(double* args, const map<string, double>& variables) const {
82     return erfc(args[0]);
83 }
84 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const85 ExpressionTreeNode Operation::Constant::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
86     return ExpressionTreeNode(new Operation::Constant(0.0));
87 }
88 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const89 ExpressionTreeNode Operation::Variable::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
90     if (variable == name)
91         return ExpressionTreeNode(new Operation::Constant(1.0));
92     return ExpressionTreeNode(new Operation::Constant(0.0));
93 }
94 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const95 ExpressionTreeNode Operation::Custom::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
96     if (function->getNumArguments() == 0)
97         return ExpressionTreeNode(new Operation::Constant(0.0));
98     ExpressionTreeNode result = ExpressionTreeNode(new Operation::Multiply(), ExpressionTreeNode(new Operation::Custom(*this, 0), children), childDerivs[0]);
99     for (int i = 1; i < getNumArguments(); i++) {
100         result = ExpressionTreeNode(new Operation::Add(),
101                                     result,
102                                     ExpressionTreeNode(new Operation::Multiply(), ExpressionTreeNode(new Operation::Custom(*this, i), children), childDerivs[i]));
103     }
104     return result;
105 }
106 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const107 ExpressionTreeNode Operation::Add::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
108     return ExpressionTreeNode(new Operation::Add(), childDerivs[0], childDerivs[1]);
109 }
110 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const111 ExpressionTreeNode Operation::Subtract::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
112     return ExpressionTreeNode(new Operation::Subtract(), childDerivs[0], childDerivs[1]);
113 }
114 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const115 ExpressionTreeNode Operation::Multiply::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
116     return ExpressionTreeNode(new Operation::Add(),
117                               ExpressionTreeNode(new Operation::Multiply(), children[0], childDerivs[1]),
118                               ExpressionTreeNode(new Operation::Multiply(), children[1], childDerivs[0]));
119 }
120 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const121 ExpressionTreeNode Operation::Divide::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
122     return ExpressionTreeNode(new Operation::Divide(),
123                               ExpressionTreeNode(new Operation::Subtract(),
124                                                  ExpressionTreeNode(new Operation::Multiply(), children[1], childDerivs[0]),
125                                                  ExpressionTreeNode(new Operation::Multiply(), children[0], childDerivs[1])),
126                               ExpressionTreeNode(new Operation::Square(), children[1]));
127 }
128 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const129 ExpressionTreeNode Operation::Power::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
130     return ExpressionTreeNode(new Operation::Add(),
131                               ExpressionTreeNode(new Operation::Multiply(),
132                                                  ExpressionTreeNode(new Operation::Multiply(),
133                                                                     children[1],
134                                                                     ExpressionTreeNode(new Operation::Power(),
135                                                                                        children[0], ExpressionTreeNode(new Operation::AddConstant(-1.0), children[1]))),
136                                                  childDerivs[0]),
137                               ExpressionTreeNode(new Operation::Multiply(),
138                                                  ExpressionTreeNode(new Operation::Multiply(),
139                                                                     ExpressionTreeNode(new Operation::Log(), children[0]),
140                                                                     ExpressionTreeNode(new Operation::Power(), children[0], children[1])),
141                                                  childDerivs[1]));
142 }
143 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const144 ExpressionTreeNode Operation::Negate::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
145     return ExpressionTreeNode(new Operation::Negate(), childDerivs[0]);
146 }
147 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const148 ExpressionTreeNode Operation::Sqrt::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
149     return ExpressionTreeNode(new Operation::Multiply(),
150                               ExpressionTreeNode(new Operation::MultiplyConstant(0.5),
151                                                  ExpressionTreeNode(new Operation::Reciprocal(),
152                                                                     ExpressionTreeNode(new Operation::Sqrt(), children[0]))),
153                               childDerivs[0]);
154 }
155 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const156 ExpressionTreeNode Operation::Exp::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
157     return ExpressionTreeNode(new Operation::Multiply(),
158                               ExpressionTreeNode(new Operation::Exp(), children[0]),
159                               childDerivs[0]);
160 }
161 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const162 ExpressionTreeNode Operation::Log::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
163     return ExpressionTreeNode(new Operation::Multiply(),
164                               ExpressionTreeNode(new Operation::Reciprocal(), children[0]),
165                               childDerivs[0]);
166 }
167 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const168 ExpressionTreeNode Operation::Sin::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
169     return ExpressionTreeNode(new Operation::Multiply(),
170                               ExpressionTreeNode(new Operation::Cos(), children[0]),
171                               childDerivs[0]);
172 }
173 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const174 ExpressionTreeNode Operation::Cos::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
175     return ExpressionTreeNode(new Operation::Multiply(),
176                               ExpressionTreeNode(new Operation::Negate(),
177                                                  ExpressionTreeNode(new Operation::Sin(), children[0])),
178                               childDerivs[0]);
179 }
180 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const181 ExpressionTreeNode Operation::Sec::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
182     return ExpressionTreeNode(new Operation::Multiply(),
183                               ExpressionTreeNode(new Operation::Multiply(),
184                                                  ExpressionTreeNode(new Operation::Sec(), children[0]),
185                                                  ExpressionTreeNode(new Operation::Tan(), children[0])),
186                               childDerivs[0]);
187 }
188 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const189 ExpressionTreeNode Operation::Csc::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
190     return ExpressionTreeNode(new Operation::Multiply(),
191                               ExpressionTreeNode(new Operation::Negate(),
192                                                  ExpressionTreeNode(new Operation::Multiply(),
193                                                                     ExpressionTreeNode(new Operation::Csc(), children[0]),
194                                                                     ExpressionTreeNode(new Operation::Cot(), children[0]))),
195                               childDerivs[0]);
196 }
197 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const198 ExpressionTreeNode Operation::Tan::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
199     return ExpressionTreeNode(new Operation::Multiply(),
200                               ExpressionTreeNode(new Operation::Square(),
201                                                  ExpressionTreeNode(new Operation::Sec(), children[0])),
202                               childDerivs[0]);
203 }
204 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const205 ExpressionTreeNode Operation::Cot::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
206     return ExpressionTreeNode(new Operation::Multiply(),
207                               ExpressionTreeNode(new Operation::Negate(),
208                                                  ExpressionTreeNode(new Operation::Square(),
209                                                                     ExpressionTreeNode(new Operation::Csc(), children[0]))),
210                               childDerivs[0]);
211 }
212 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const213 ExpressionTreeNode Operation::Asin::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
214     return ExpressionTreeNode(new Operation::Multiply(),
215                               ExpressionTreeNode(new Operation::Reciprocal(),
216                                                  ExpressionTreeNode(new Operation::Sqrt(),
217                                                                     ExpressionTreeNode(new Operation::Subtract(),
218                                                                                        ExpressionTreeNode(new Operation::Constant(1.0)),
219                                                                                        ExpressionTreeNode(new Operation::Square(), children[0])))),
220                               childDerivs[0]);
221 }
222 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const223 ExpressionTreeNode Operation::Acos::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
224     return ExpressionTreeNode(new Operation::Multiply(),
225                               ExpressionTreeNode(new Operation::Negate(),
226                                                  ExpressionTreeNode(new Operation::Reciprocal(),
227                                                                     ExpressionTreeNode(new Operation::Sqrt(),
228                                                                                        ExpressionTreeNode(new Operation::Subtract(),
229                                                                                                           ExpressionTreeNode(new Operation::Constant(1.0)),
230                                                                                                           ExpressionTreeNode(new Operation::Square(), children[0]))))),
231                               childDerivs[0]);
232 }
233 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const234 ExpressionTreeNode Operation::Atan::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
235     return ExpressionTreeNode(new Operation::Multiply(),
236                               ExpressionTreeNode(new Operation::Reciprocal(),
237                                                  ExpressionTreeNode(new Operation::AddConstant(1.0),
238                                                                     ExpressionTreeNode(new Operation::Square(), children[0]))),
239                               childDerivs[0]);
240 }
241 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const242 ExpressionTreeNode Operation::Sinh::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
243     return ExpressionTreeNode(new Operation::Multiply(),
244                               ExpressionTreeNode(new Operation::Cosh(),
245                                                  children[0]),
246                               childDerivs[0]);
247 }
248 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const249 ExpressionTreeNode Operation::Cosh::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
250     return ExpressionTreeNode(new Operation::Multiply(),
251                               ExpressionTreeNode(new Operation::Sinh(),
252                                                  children[0]),
253                               childDerivs[0]);
254 }
255 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const256 ExpressionTreeNode Operation::Tanh::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
257     return ExpressionTreeNode(new Operation::Multiply(),
258                               ExpressionTreeNode(new Operation::Subtract(),
259                                                  ExpressionTreeNode(new Operation::Constant(1.0)),
260                                                  ExpressionTreeNode(new Operation::Square(),
261                                                                     ExpressionTreeNode(new Operation::Tanh(), children[0]))),
262                               childDerivs[0]);
263 }
264 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const265 ExpressionTreeNode Operation::Erf::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
266     return ExpressionTreeNode(new Operation::Multiply(),
267                               ExpressionTreeNode(new Operation::Multiply(),
268                                                  ExpressionTreeNode(new Operation::Constant(2.0/sqrt(M_PI))),
269                                                  ExpressionTreeNode(new Operation::Exp(),
270                                                                     ExpressionTreeNode(new Operation::Negate(),
271                                                                                        ExpressionTreeNode(new Operation::Square(), children[0])))),
272                               childDerivs[0]);
273 }
274 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const275 ExpressionTreeNode Operation::Erfc::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
276     return ExpressionTreeNode(new Operation::Multiply(),
277                               ExpressionTreeNode(new Operation::Multiply(),
278                                                  ExpressionTreeNode(new Operation::Constant(-2.0/sqrt(M_PI))),
279                                                  ExpressionTreeNode(new Operation::Exp(),
280                                                                     ExpressionTreeNode(new Operation::Negate(),
281                                                                                        ExpressionTreeNode(new Operation::Square(), children[0])))),
282                               childDerivs[0]);
283 }
284 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const285 ExpressionTreeNode Operation::Step::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
286     return ExpressionTreeNode(new Operation::Delta(),children[0]);
287 }
288 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const289 ExpressionTreeNode Operation::Delta::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
290     return ExpressionTreeNode(new Operation::Nandelta(), children[0]);
291 }
292 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const293 ExpressionTreeNode Operation::Nandelta::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
294     return ExpressionTreeNode(new Operation::Nandelta(), children[0]);
295 }
296 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const297 ExpressionTreeNode Operation::Square::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
298     return ExpressionTreeNode(new Operation::Multiply(),
299                               ExpressionTreeNode(new Operation::MultiplyConstant(2.0),
300                                                  children[0]),
301                               childDerivs[0]);
302 }
303 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const304 ExpressionTreeNode Operation::Cube::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
305     return ExpressionTreeNode(new Operation::Multiply(),
306                               ExpressionTreeNode(new Operation::MultiplyConstant(3.0),
307                                                  ExpressionTreeNode(new Operation::Square(), children[0])),
308                               childDerivs[0]);
309 }
310 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const311 ExpressionTreeNode Operation::Reciprocal::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
312     return ExpressionTreeNode(new Operation::Multiply(),
313                               ExpressionTreeNode(new Operation::Negate(),
314                                                  ExpressionTreeNode(new Operation::Reciprocal(),
315                                                                     ExpressionTreeNode(new Operation::Square(), children[0]))),
316                               childDerivs[0]);
317 }
318 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const319 ExpressionTreeNode Operation::AddConstant::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
320     return childDerivs[0];
321 }
322 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const323 ExpressionTreeNode Operation::MultiplyConstant::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
324     return ExpressionTreeNode(new Operation::MultiplyConstant(value),
325                               childDerivs[0]);
326 }
327 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const328 ExpressionTreeNode Operation::PowerConstant::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
329     return ExpressionTreeNode(new Operation::Multiply(),
330                               ExpressionTreeNode(new Operation::MultiplyConstant(value),
331                                                  ExpressionTreeNode(new Operation::PowerConstant(value-1),
332                                                                     children[0])),
333                               childDerivs[0]);
334 }
335 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const336 ExpressionTreeNode Operation::Min::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
337     ExpressionTreeNode step(new Operation::Step(),
338                             ExpressionTreeNode(new Operation::Subtract(), children[0], children[1]));
339     return ExpressionTreeNode(new Operation::Subtract(),
340                               ExpressionTreeNode(new Operation::Multiply(), childDerivs[1], step),
341                               ExpressionTreeNode(new Operation::Multiply(), childDerivs[0],
342                                                  ExpressionTreeNode(new Operation::AddConstant(-1), step)));
343 }
344 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const345 ExpressionTreeNode Operation::Max::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
346     ExpressionTreeNode step(new Operation::Step(),
347                             ExpressionTreeNode(new Operation::Subtract(), children[0], children[1]));
348     return ExpressionTreeNode(new Operation::Subtract(),
349                               ExpressionTreeNode(new Operation::Multiply(), childDerivs[0], step),
350                               ExpressionTreeNode(new Operation::Multiply(), childDerivs[1],
351                                                  ExpressionTreeNode(new Operation::AddConstant(-1), step)));
352 }
353 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const354 ExpressionTreeNode Operation::Abs::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
355     ExpressionTreeNode step(new Operation::Step(), children[0]);
356     return ExpressionTreeNode(new Operation::Multiply(),
357                               childDerivs[0],
358                               ExpressionTreeNode(new Operation::AddConstant(-1),
359                                                  ExpressionTreeNode(new Operation::MultiplyConstant(2), step)));
360 }
361 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const362 ExpressionTreeNode Operation::Floor::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
363     return ExpressionTreeNode(new Operation::Constant(0.0));
364 }
365 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const366 ExpressionTreeNode Operation::Ceil::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
367     return ExpressionTreeNode(new Operation::Constant(0.0));
368 }
369 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const370 ExpressionTreeNode Operation::Select::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
371     vector<ExpressionTreeNode> derivChildren;
372     derivChildren.push_back(children[0]);
373     derivChildren.push_back(childDerivs[1]);
374     derivChildren.push_back(childDerivs[2]);
375     return ExpressionTreeNode(new Operation::Select(), derivChildren);
376 }
377 
378 #define LEPTON_CONST(x) ExpressionTreeNode(new Operation::Constant(x))
379 #define LEPTON_OP1(name,x) ExpressionTreeNode(new Operation::name(),x)
380 #define LEPTON_OP2(name,x,y) ExpressionTreeNode(new Operation::name(),x,y)
381 #define LEPTON_ADD_CONST(x,y) ExpressionTreeNode(new Operation::AddConstant(x),y)
382 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const383 ExpressionTreeNode Operation::Acot::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
384     return
385       LEPTON_OP2(Multiply,
386         LEPTON_OP1(Negate,
387           LEPTON_OP1(Reciprocal,
388             LEPTON_ADD_CONST(1.0,
389               LEPTON_OP1(Square,children[0])
390             )
391           )
392         ),
393         childDerivs[0]
394       );
395 }
396 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const397 ExpressionTreeNode Operation::Asec::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
398     return
399       LEPTON_OP2(Multiply,
400         LEPTON_OP1(Reciprocal,
401           LEPTON_OP2(Multiply,
402             LEPTON_OP1(Abs,children[0]),
403             LEPTON_OP1(Sqrt,
404               LEPTON_OP2(Subtract,
405                 LEPTON_OP1(Square,children[0]),
406                 LEPTON_CONST(1.0)
407               )
408             )
409           )
410         ),
411         childDerivs[0]
412       );
413 }
414 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const415 ExpressionTreeNode Operation::Acsc::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
416     return
417     LEPTON_OP2(Multiply,
418       LEPTON_OP1(Negate,
419         LEPTON_OP1(Reciprocal,
420           LEPTON_OP2(Multiply,
421             LEPTON_OP1(Abs,children[0]),
422             LEPTON_OP1(Sqrt,
423               LEPTON_OP2(Subtract,
424                 LEPTON_OP1(Square,children[0]),
425                 LEPTON_CONST(1.0)
426               )
427             )
428           )
429         )
430       ),
431       childDerivs[0]
432     );
433 }
434 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const435 ExpressionTreeNode Operation::Coth::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
436     return
437     LEPTON_OP2(Multiply,
438       LEPTON_OP2(Subtract,
439         LEPTON_CONST(1.0),
440         LEPTON_OP1(Square,
441           LEPTON_OP1(Coth,children[0])
442         )
443       ),
444       childDerivs[0]
445     );
446 }
447 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const448 ExpressionTreeNode Operation::Sech::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
449     return
450     LEPTON_OP2(Multiply,
451       LEPTON_OP1(Negate,
452         LEPTON_OP2(Multiply,
453           LEPTON_OP1(Tanh,children[0]),
454           LEPTON_OP1(Sech,children[0])
455         )
456       ),
457       childDerivs[0]
458     );
459 }
460 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const461 ExpressionTreeNode Operation::Csch::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
462     return
463     LEPTON_OP2(Multiply,
464       LEPTON_OP1(Negate,
465         LEPTON_OP2(Multiply,
466           LEPTON_OP1(Coth,children[0]),
467           LEPTON_OP1(Csch,children[0])
468         )
469       ),
470       childDerivs[0]
471     );
472 }
473 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const474 ExpressionTreeNode Operation::Acosh::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
475     return
476     LEPTON_OP2(Multiply,
477       LEPTON_OP1(Reciprocal,
478         LEPTON_OP1(Sqrt,
479           LEPTON_OP2(Subtract,
480             LEPTON_OP1(Square,children[0]),
481             LEPTON_CONST(1.0)
482           )
483         )
484       ),
485       childDerivs[0]
486     );
487 }
488 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const489 ExpressionTreeNode Operation::Atanh::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
490     return
491     LEPTON_OP2(Multiply,
492       LEPTON_OP1(Reciprocal,
493         LEPTON_OP2(Subtract,
494           LEPTON_CONST(1.0),
495           LEPTON_OP1(Square,children[0])
496         )
497       ),
498       childDerivs[0]
499     );
500 }
501 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const502 ExpressionTreeNode Operation::Asinh::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
503     return
504     LEPTON_OP2(Multiply,
505       LEPTON_OP1(Reciprocal,
506         LEPTON_OP1(Sqrt,
507           LEPTON_ADD_CONST(1.0,
508             LEPTON_OP1(Square,children[0])
509           )
510         )
511       ),
512       childDerivs[0]
513     );
514 }
515 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const516 ExpressionTreeNode Operation::Acoth::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
517     return
518     LEPTON_OP2(Multiply,
519       LEPTON_OP1(Reciprocal,
520         LEPTON_OP2(Subtract,
521           LEPTON_CONST(1.0),
522           LEPTON_OP1(Square,children[0])
523         )
524       ),
525       childDerivs[0]
526     );
527 }
528 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const529 ExpressionTreeNode Operation::Asech::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
530     return
531     LEPTON_OP2(Multiply,
532       LEPTON_OP1(Negate,
533         LEPTON_OP1(Reciprocal,
534           LEPTON_OP2(Multiply,
535             children[0],
536             LEPTON_OP2(Multiply,
537               LEPTON_ADD_CONST(1.0,
538                 children[0]
539               ),
540               LEPTON_OP1(Sqrt,
541                 LEPTON_OP2(Divide,
542                   LEPTON_OP2(Subtract,
543                     LEPTON_CONST(1.0),
544                     children[0]
545                   ),
546                   LEPTON_ADD_CONST(1.0,
547                     children[0]
548                   )
549                 )
550               )
551             )
552           )
553         )
554       ),
555       childDerivs[0]
556     );
557 }
558 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const559 ExpressionTreeNode Operation::Acsch::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
560     return
561     LEPTON_OP2(Multiply,
562       LEPTON_OP1(Negate,
563         LEPTON_OP1(Reciprocal,
564           LEPTON_OP2(Multiply,
565             LEPTON_OP1(Square,children[0]),
566             LEPTON_OP1(Sqrt,
567               LEPTON_ADD_CONST(1.0,
568                 LEPTON_OP1(Reciprocal,
569                   LEPTON_OP1(Square,children[0])
570                 )
571               )
572             )
573           )
574         )
575       ),
576       childDerivs[0]
577     );
578 }
579 
differentiate(const std::vector<ExpressionTreeNode> & children,const std::vector<ExpressionTreeNode> & childDerivs,const std::string & variable) const580 ExpressionTreeNode Operation::Atan2::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
581     return
582     LEPTON_OP2(Divide,
583       LEPTON_OP2(Subtract,
584         LEPTON_OP2(Multiply, children[1], childDerivs[0]),
585         LEPTON_OP2(Multiply, children[0], childDerivs[1])
586       ),
587       LEPTON_OP2(Add,
588         LEPTON_OP1(Square, children[0]),
589         LEPTON_OP1(Square, children[1])
590       )
591     );
592 }
593 
594 }
595