source: trunk/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Converters/DerivativeCalculator.cs

Last change on this file was 17902, checked in by gkronber, 4 months ago

#3073 merged reintegration branch to trunk

File size: 12.7 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using System;
23using System.Collections.Generic;
24using System.Linq;
25using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
26
27namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
28  public static class DerivativeCalculator {
29    public static ISymbolicExpressionTree Derive(ISymbolicExpressionTree tree, string variableName) {
30      if (tree.Root.SubtreeCount != 1)
31        throw new NotImplementedException("Derive is not implemented for symbolic expressions with automatically defined functions (ADF)");
32      if (tree.Root.GetSubtree(0).SubtreeCount != 1)
33        throw new NotImplementedException("Derive is not implemented for multi-variate symbolic expressions");
34      var mainBranch = tree.Root.GetSubtree(0).GetSubtree(0);
35      var root = new ProgramRootSymbol().CreateTreeNode();
36      root.AddSubtree(new StartSymbol().CreateTreeNode());
37      var dTree = TreeSimplifier.GetSimplifiedTree(Derive(mainBranch, variableName));
38      //var dTree = Derive(mainBranch, variableName);
39      root.GetSubtree(0).AddSubtree(dTree);
40      return new SymbolicExpressionTree(root);
41    }
42
43    private static readonly Constant constantSy = new Constant();
44    private static readonly Addition addSy = new Addition();
45    private static readonly Subtraction subSy = new Subtraction();
46    private static readonly Multiplication mulSy = new Multiplication();
47    private static readonly Division divSy = new Division();
48    private static readonly Cosine cosSy = new Cosine();
49    private static readonly Square sqrSy = new Square();
50    private static readonly Absolute absSy = new Absolute();
51    private static readonly SquareRoot sqrtSy = new SquareRoot();
52
53    public static ISymbolicExpressionTreeNode Derive(ISymbolicExpressionTreeNode branch, string variableName) {
54      if (branch.Symbol is Constant) {
55        return CreateConstant(0.0);
56      }
57      if (branch.Symbol is Variable) {
58        var varNode = branch as VariableTreeNode;
59        if (varNode.VariableName == variableName) {
60          return CreateConstant(varNode.Weight);
61        } else {
62          return CreateConstant(0.0);
63        }
64      }
65      if (branch.Symbol is Addition) {
66        var sum = addSy.CreateTreeNode();
67        foreach (var subTree in branch.Subtrees) {
68          sum.AddSubtree(Derive(subTree, variableName));
69        }
70        return sum;
71      }
72      if (branch.Symbol is Subtraction) {
73        var sum = subSy.CreateTreeNode();
74        foreach (var subTree in branch.Subtrees) {
75          sum.AddSubtree(Derive(subTree, variableName));
76        }
77        return sum;
78      }
79      if (branch.Symbol is Multiplication) {
80        // (f * g)' = f'*g + f*g'
81        // for multiple factors: (f * g * h)' = ((f*g) * h)' = (f*g)' * h + (f*g) * h'
82
83        if (branch.SubtreeCount >= 2) {
84          var f = (ISymbolicExpressionTreeNode)branch.GetSubtree(0).Clone();
85          var fprime = Derive(f, variableName);
86          for (int i = 1; i < branch.SubtreeCount; i++) {
87            var g = (ISymbolicExpressionTreeNode)branch.GetSubtree(i).Clone();
88            var fg = Product((ISymbolicExpressionTreeNode)f.Clone(), (ISymbolicExpressionTreeNode)g.Clone());
89            var gPrime = Derive(g, variableName);
90            var fgPrime = Sum(Product(fprime, g), Product(gPrime, f));
91            // prepare for next iteration
92            f = fg;
93            fprime = fgPrime;
94          }
95          return fprime;
96        } else
97          // multiplication with only one argument has no effect -> derive the argument
98          return Derive(branch.GetSubtree(0), variableName);
99      }
100      if (branch.Symbol is Division) {
101        // (f/g)' = (f'g - g'f) / g²
102        if (branch.SubtreeCount == 1) {
103          var g = (ISymbolicExpressionTreeNode)branch.GetSubtree(0).Clone();
104          var gPrime = Product(CreateConstant(-1.0), Derive(g, variableName));
105          var sqrNode = new Square().CreateTreeNode();
106          sqrNode.AddSubtree(g);
107          return Div(gPrime, sqrNode);
108        } else {
109          // for two subtrees:
110          // (f/g)' = (f'g - fg')/g²
111
112          // if there are more than 2 subtrees
113          // div(x,y,z) is interpretered as (x/y)/z
114          // which is the same as x / (y*z)
115
116          // --> make a product of all but the first subtree and differentiate as for the 2-argument case above
117          var f = (ISymbolicExpressionTreeNode)branch.GetSubtree(0).Clone();
118          var g = Product(branch.Subtrees.Skip(1).Select(n => (ISymbolicExpressionTreeNode)n.Clone()));
119          var fprime = Derive(f, variableName);
120          var gprime = Derive(g, variableName);
121          var gSqr = Square(g);
122          return Div(Subtract(Product(fprime, g), Product(f, gprime)), gSqr);
123        }
124      }
125      if (branch.Symbol is Logarithm) {
126        var f = (ISymbolicExpressionTreeNode)branch.GetSubtree(0).Clone();
127        return Product(Div(CreateConstant(1.0), f), Derive(f, variableName));
128      }
129      if (branch.Symbol is Exponential) {
130        var f = (ISymbolicExpressionTreeNode)branch.Clone();
131        return Product(f, Derive(branch.GetSubtree(0), variableName));
132      }
133      if (branch.Symbol is Square) {
134        var f = (ISymbolicExpressionTreeNode)branch.GetSubtree(0).Clone();
135        return Product(Product(CreateConstant(2.0), f), Derive(f, variableName));
136      }
137      if (branch.Symbol is SquareRoot) {
138        var f = (ISymbolicExpressionTreeNode)branch.Clone();
139        var u = (ISymbolicExpressionTreeNode)branch.GetSubtree(0).Clone();
140        return Product(Div(CreateConstant(1.0), Product(CreateConstant(2.0), f)), Derive(u, variableName));
141      }
142      if (branch.Symbol is CubeRoot) {
143        var f = (ISymbolicExpressionTreeNode)branch.Clone();
144        var u = (ISymbolicExpressionTreeNode)branch.GetSubtree(0).Clone();
145        return Product(Div(CreateConstant(1.0), Product(CreateConstant(3.0), Square(f))), Derive(u, variableName));  // 1/3 1/cbrt(f(x))^2 d/dx f(x)
146      }
147      if (branch.Symbol is Cube) {
148        var f = (ISymbolicExpressionTreeNode)branch.GetSubtree(0).Clone();
149        return Product(Product(CreateConstant(3.0), Square(f)), Derive(f, variableName));
150      }
151      if (branch.Symbol is Power) {
152        // HL evaluators handle power strangely (exponent is rounded to an integer)
153        // here we only support the case when the exponent is a constant integer
154        var exponent = branch.GetSubtree(1) as ConstantTreeNode;
155        if (exponent != null && Math.Truncate(exponent.Value) == exponent.Value) {
156          var newPower = (ISymbolicExpressionTreeNode)branch.Clone();
157          var f = (ISymbolicExpressionTreeNode)newPower.GetSubtree(0).Clone();
158          var newExponent = (ConstantTreeNode)newPower.GetSubtree(1);
159          newExponent.Value -= 1;
160          return Product(Product(CreateConstant(exponent.Value), newPower), Derive(f, variableName));
161        } else throw new NotSupportedException("Cannot derive non-integer powers");
162      }
163      if (branch.Symbol is Absolute) {
164        var f = (ISymbolicExpressionTreeNode)branch.GetSubtree(0).Clone();
165        var absf = Abs((ISymbolicExpressionTreeNode)f.Clone());
166        return Product(Div(f, absf), Derive(f, variableName));
167      }
168      if (branch.Symbol is AnalyticQuotient) {
169        // aq(a(x), b(x)) = a(x) / sqrt(b(x)²+1)
170        var a = (ISymbolicExpressionTreeNode)branch.GetSubtree(0).Clone();
171        var b = (ISymbolicExpressionTreeNode)branch.GetSubtree(1).Clone();
172
173        var definition = Div(a, SquareRoot(Sum(Square(b), CreateConstant(1.0))));
174        return Derive(definition, variableName);
175      }
176      if (branch.Symbol is Sine) {
177        var u = (ISymbolicExpressionTreeNode)branch.GetSubtree(0).Clone();
178        var cos = (new Cosine()).CreateTreeNode();
179        cos.AddSubtree(u);
180        return Product(cos, Derive(u, variableName));
181      }
182      if (branch.Symbol is Cosine) {
183        var u = (ISymbolicExpressionTreeNode)branch.GetSubtree(0).Clone();
184        var sin = (new Sine()).CreateTreeNode();
185        sin.AddSubtree(u);
186        return Product(CreateConstant(-1.0), Product(sin, Derive(u, variableName)));
187      }
188      if (branch.Symbol is Tangent) {
189        // tan(x)' = 1 / cos²(x)
190        var fxp = Derive(branch.GetSubtree(0), variableName);
191        var u = (ISymbolicExpressionTreeNode)branch.GetSubtree(0).Clone();
192        return Div(fxp, Square(Cosine(u)));
193      }
194      if (branch.Symbol is HyperbolicTangent) {
195        // tanh(f(x))' = f(x)'sech²(f(x)) = f(x)'(1 - tanh²(f(x)))
196        var fxp = Derive(branch.GetSubtree(0), variableName);
197        var tanh = (ISymbolicExpressionTreeNode)branch.Clone();
198        return Product(fxp, Subtract(CreateConstant(1.0), Square(tanh)));
199      }
200      throw new NotSupportedException(string.Format("Symbol {0} is not supported.", branch.Symbol));
201    }
202
203
204    private static ISymbolicExpressionTreeNode Product(ISymbolicExpressionTreeNode f, ISymbolicExpressionTreeNode g) {
205      var product = mulSy.CreateTreeNode();
206      product.AddSubtree(f);
207      product.AddSubtree(g);
208      return product;
209    }
210    private static ISymbolicExpressionTreeNode Product(IEnumerable<ISymbolicExpressionTreeNode> fs) {
211      var product = mulSy.CreateTreeNode();
212      foreach (var f in fs) product.AddSubtree(f);
213      return product;
214    }
215    private static ISymbolicExpressionTreeNode Div(ISymbolicExpressionTreeNode f, ISymbolicExpressionTreeNode g) {
216      var div = divSy.CreateTreeNode();
217      div.AddSubtree(f);
218      div.AddSubtree(g);
219      return div;
220    }
221
222    private static ISymbolicExpressionTreeNode Sum(ISymbolicExpressionTreeNode f, ISymbolicExpressionTreeNode g) {
223      var sum = addSy.CreateTreeNode();
224      sum.AddSubtree(f);
225      sum.AddSubtree(g);
226      return sum;
227    }
228    private static ISymbolicExpressionTreeNode Subtract(ISymbolicExpressionTreeNode f, ISymbolicExpressionTreeNode g) {
229      var sum = subSy.CreateTreeNode();
230      sum.AddSubtree(f);
231      sum.AddSubtree(g);
232      return sum;
233    }
234    private static ISymbolicExpressionTreeNode Cosine(ISymbolicExpressionTreeNode f) {
235      var cos = cosSy.CreateTreeNode();
236      cos.AddSubtree(f);
237      return cos;
238    }
239    private static ISymbolicExpressionTreeNode Abs(ISymbolicExpressionTreeNode f) {
240      var abs = absSy.CreateTreeNode();
241      abs.AddSubtree(f);
242      return abs;
243    }
244    private static ISymbolicExpressionTreeNode Square(ISymbolicExpressionTreeNode f) {
245      var sqr = sqrSy.CreateTreeNode();
246      sqr.AddSubtree(f);
247      return sqr;
248    }
249    private static ISymbolicExpressionTreeNode SquareRoot(ISymbolicExpressionTreeNode f) {
250      var sqrt = sqrtSy.CreateTreeNode();
251      sqrt.AddSubtree(f);
252      return sqrt;
253    }
254
255    private static ISymbolicExpressionTreeNode CreateConstant(double v) {
256      var constNode = (ConstantTreeNode)constantSy.CreateTreeNode();
257      constNode.Value = v;
258      return constNode;
259    }
260
261    public static bool IsCompatible(ISymbolicExpressionTree tree) {
262      var containsUnknownSymbol = (
263        from n in tree.Root.GetSubtree(0).IterateNodesPrefix()
264        where
265          !(n.Symbol is Variable) &&
266          !(n.Symbol is Constant) &&
267          !(n.Symbol is Addition) &&
268          !(n.Symbol is Subtraction) &&
269          !(n.Symbol is Multiplication) &&
270          !(n.Symbol is Division) &&
271          !(n.Symbol is Logarithm) &&
272          !(n.Symbol is Exponential) &&
273          !(n.Symbol is Square) &&
274          !(n.Symbol is SquareRoot) &&
275          !(n.Symbol is Cube) &&
276          !(n.Symbol is CubeRoot) &&
277          !(n.Symbol is Power) &&
278          !(n.Symbol is Absolute) &&
279          !(n.Symbol is AnalyticQuotient) &&
280          !(n.Symbol is HyperbolicTangent) &&
281          !(n.Symbol is Sine) &&
282          !(n.Symbol is Cosine) &&
283          !(n.Symbol is Tangent) &&
284          !(n.Symbol is StartSymbol)
285        select n).Any();
286      return !containsUnknownSymbol;
287    }
288  }
289}
Note: See TracBrowser for help on using the repository browser.