source: stable/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Converters/TreeToAutoDiffTermConverter.cs @ 18015

Last change on this file since 18015 was 18015, checked in by gkronber, 2 months ago

#3095: merged r17817 from trunk to stable

File size: 14.9 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 System.Runtime.Serialization;
26using AutoDiff;
27using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
28
29namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
30  public class TreeToAutoDiffTermConverter {
31    public delegate double ParametricFunction(double[] vars, double[] @params);
32
33    public delegate Tuple<double[], double> ParametricFunctionGradient(double[] vars, double[] @params);
34
35    #region helper class
36    public class DataForVariable {
37      public readonly string variableName;
38      public readonly string variableValue; // for factor vars
39      public readonly int lag;
40
41      public DataForVariable(string varName, string varValue, int lag) {
42        this.variableName = varName;
43        this.variableValue = varValue;
44        this.lag = lag;
45      }
46
47      public override bool Equals(object obj) {
48        var other = obj as DataForVariable;
49        if (other == null) return false;
50        return other.variableName.Equals(this.variableName) &&
51               other.variableValue.Equals(this.variableValue) &&
52               other.lag == this.lag;
53      }
54
55      public override int GetHashCode() {
56        return variableName.GetHashCode() ^ variableValue.GetHashCode() ^ lag;
57      }
58    }
59    #endregion
60
61    #region derivations of functions
62    // create function factory for arctangent
63    private static readonly Func<Term, UnaryFunc> arctan = UnaryFunc.Factory(
64      eval: Math.Atan,
65      diff: x => 1 / (1 + x * x));
66
67    private static readonly Func<Term, UnaryFunc> sin = UnaryFunc.Factory(
68      eval: Math.Sin,
69      diff: Math.Cos);
70
71    private static readonly Func<Term, UnaryFunc> cos = UnaryFunc.Factory(
72      eval: Math.Cos,
73      diff: x => -Math.Sin(x));
74
75    private static readonly Func<Term, UnaryFunc> tan = UnaryFunc.Factory(
76      eval: Math.Tan,
77      diff: x => 1 + Math.Tan(x) * Math.Tan(x));
78    private static readonly Func<Term, UnaryFunc> tanh = UnaryFunc.Factory(
79      eval: Math.Tanh,
80      diff: x => 1 - Math.Tanh(x) * Math.Tanh(x));
81    private static readonly Func<Term, UnaryFunc> erf = UnaryFunc.Factory(
82      eval: alglib.errorfunction,
83      diff: x => 2.0 * Math.Exp(-(x * x)) / Math.Sqrt(Math.PI));
84
85    private static readonly Func<Term, UnaryFunc> norm = UnaryFunc.Factory(
86      eval: alglib.normaldistribution,
87      diff: x => -(Math.Exp(-(x * x)) * Math.Sqrt(Math.Exp(x * x)) * x) / Math.Sqrt(2 * Math.PI));
88
89    private static readonly Func<Term, UnaryFunc> abs = UnaryFunc.Factory(
90      eval: Math.Abs,
91      diff: x => Math.Sign(x)
92      );
93
94    private static readonly Func<Term, UnaryFunc> cbrt = UnaryFunc.Factory(
95      eval: x => x < 0 ? -Math.Pow(-x, 1.0 / 3) : Math.Pow(x, 1.0 / 3),
96      diff: x => { var cbrt_x = x < 0 ? -Math.Pow(-x, 1.0 / 3) : Math.Pow(x, 1.0 / 3); return 1.0 / (3 * cbrt_x * cbrt_x); }
97      );
98
99
100
101    #endregion
102
103    public static bool TryConvertToAutoDiff(ISymbolicExpressionTree tree, bool makeVariableWeightsVariable, bool addLinearScalingTerms,
104      out List<DataForVariable> parameters, out double[] initialConstants,
105      out ParametricFunction func,
106      out ParametricFunctionGradient func_grad) {
107
108      // use a transformator object which holds the state (variable list, parameter list, ...) for recursive transformation of the tree
109      var transformator = new TreeToAutoDiffTermConverter(makeVariableWeightsVariable, addLinearScalingTerms);
110      AutoDiff.Term term;
111      try {
112        term = transformator.ConvertToAutoDiff(tree.Root.GetSubtree(0));
113        var parameterEntries = transformator.parameters.ToArray(); // guarantee same order for keys and values
114        var compiledTerm = term.Compile(transformator.variables.ToArray(),
115          parameterEntries.Select(kvp => kvp.Value).ToArray());
116        parameters = new List<DataForVariable>(parameterEntries.Select(kvp => kvp.Key));
117        initialConstants = transformator.initialConstants.ToArray();
118        func = (vars, @params) => compiledTerm.Evaluate(vars, @params);
119        func_grad = (vars, @params) => compiledTerm.Differentiate(vars, @params);
120        return true;
121      } catch (ConversionException) {
122        func = null;
123        func_grad = null;
124        parameters = null;
125        initialConstants = null;
126      }
127      return false;
128    }
129
130    // state for recursive transformation of trees
131    private readonly
132    List<double> initialConstants;
133    private readonly Dictionary<DataForVariable, AutoDiff.Variable> parameters;
134    private readonly List<AutoDiff.Variable> variables;
135    private readonly bool makeVariableWeightsVariable;
136    private readonly bool addLinearScalingTerms;
137
138    private TreeToAutoDiffTermConverter(bool makeVariableWeightsVariable, bool addLinearScalingTerms) {
139      this.makeVariableWeightsVariable = makeVariableWeightsVariable;
140      this.addLinearScalingTerms = addLinearScalingTerms;
141      this.initialConstants = new List<double>();
142      this.parameters = new Dictionary<DataForVariable, AutoDiff.Variable>();
143      this.variables = new List<AutoDiff.Variable>();
144    }
145
146    private AutoDiff.Term ConvertToAutoDiff(ISymbolicExpressionTreeNode node) {
147      if (node.Symbol is Constant) {
148        initialConstants.Add(((ConstantTreeNode)node).Value);
149        var var = new AutoDiff.Variable();
150        variables.Add(var);
151        return var;
152      }
153      if (node.Symbol is Variable || node.Symbol is BinaryFactorVariable) {
154        var varNode = node as VariableTreeNodeBase;
155        var factorVarNode = node as BinaryFactorVariableTreeNode;
156        // factor variable values are only 0 or 1 and set in x accordingly
157        var varValue = factorVarNode != null ? factorVarNode.VariableValue : string.Empty;
158        var par = FindOrCreateParameter(parameters, varNode.VariableName, varValue);
159
160        if (makeVariableWeightsVariable) {
161          initialConstants.Add(varNode.Weight);
162          var w = new AutoDiff.Variable();
163          variables.Add(w);
164          return AutoDiff.TermBuilder.Product(w, par);
165        } else {
166          return varNode.Weight * par;
167        }
168      }
169      if (node.Symbol is FactorVariable) {
170        var factorVarNode = node as FactorVariableTreeNode;
171        var products = new List<Term>();
172        foreach (var variableValue in factorVarNode.Symbol.GetVariableValues(factorVarNode.VariableName)) {
173          var par = FindOrCreateParameter(parameters, factorVarNode.VariableName, variableValue);
174
175          initialConstants.Add(factorVarNode.GetValue(variableValue));
176          var wVar = new AutoDiff.Variable();
177          variables.Add(wVar);
178
179          products.Add(AutoDiff.TermBuilder.Product(wVar, par));
180        }
181        return AutoDiff.TermBuilder.Sum(products);
182      }
183      if (node.Symbol is LaggedVariable) {
184        var varNode = node as LaggedVariableTreeNode;
185        var par = FindOrCreateParameter(parameters, varNode.VariableName, string.Empty, varNode.Lag);
186
187        if (makeVariableWeightsVariable) {
188          initialConstants.Add(varNode.Weight);
189          var w = new AutoDiff.Variable();
190          variables.Add(w);
191          return AutoDiff.TermBuilder.Product(w, par);
192        } else {
193          return varNode.Weight * par;
194        }
195      }
196      if (node.Symbol is Addition) {
197        List<AutoDiff.Term> terms = new List<Term>();
198        foreach (var subTree in node.Subtrees) {
199          terms.Add(ConvertToAutoDiff(subTree));
200        }
201        return AutoDiff.TermBuilder.Sum(terms);
202      }
203      if (node.Symbol is Subtraction) {
204        List<AutoDiff.Term> terms = new List<Term>();
205        for (int i = 0; i < node.SubtreeCount; i++) {
206          AutoDiff.Term t = ConvertToAutoDiff(node.GetSubtree(i));
207          if (i > 0) t = -t;
208          terms.Add(t);
209        }
210        if (terms.Count == 1) return -terms[0];
211        else return AutoDiff.TermBuilder.Sum(terms);
212      }
213      if (node.Symbol is Multiplication) {
214        List<AutoDiff.Term> terms = new List<Term>();
215        foreach (var subTree in node.Subtrees) {
216          terms.Add(ConvertToAutoDiff(subTree));
217        }
218        if (terms.Count == 1) return terms[0];
219        else return terms.Aggregate((a, b) => new AutoDiff.Product(a, b));
220      }
221      if (node.Symbol is Division) {
222        List<AutoDiff.Term> terms = new List<Term>();
223        foreach (var subTree in node.Subtrees) {
224          terms.Add(ConvertToAutoDiff(subTree));
225        }
226        if (terms.Count == 1) return 1.0 / terms[0];
227        else return terms.Aggregate((a, b) => new AutoDiff.Product(a, 1.0 / b));
228      }
229      if (node.Symbol is Absolute) {
230        var x1 = ConvertToAutoDiff(node.GetSubtree(0));
231        return abs(x1);
232      }
233      if (node.Symbol is AnalyticQuotient) {
234        var x1 = ConvertToAutoDiff(node.GetSubtree(0));
235        var x2 = ConvertToAutoDiff(node.GetSubtree(1));
236        return x1 / (TermBuilder.Power(1 + x2 * x2, 0.5));
237      }
238      if (node.Symbol is Logarithm) {
239        return AutoDiff.TermBuilder.Log(
240          ConvertToAutoDiff(node.GetSubtree(0)));
241      }
242      if (node.Symbol is Exponential) {
243        return AutoDiff.TermBuilder.Exp(
244          ConvertToAutoDiff(node.GetSubtree(0)));
245      }
246      if (node.Symbol is Square) {
247        return AutoDiff.TermBuilder.Power(
248          ConvertToAutoDiff(node.GetSubtree(0)), 2.0);
249      }
250      if (node.Symbol is SquareRoot) {
251        return AutoDiff.TermBuilder.Power(
252          ConvertToAutoDiff(node.GetSubtree(0)), 0.5);
253      }
254      if (node.Symbol is Cube) {
255        return AutoDiff.TermBuilder.Power(
256          ConvertToAutoDiff(node.GetSubtree(0)), 3.0);
257      }
258      if (node.Symbol is CubeRoot) {
259        return cbrt(ConvertToAutoDiff(node.GetSubtree(0)));
260      }
261      if (node.Symbol is Power) {
262        var powerNode = node.GetSubtree(1) as ConstantTreeNode;
263        if (powerNode == null)
264          throw new NotSupportedException("Only integer powers are allowed in parameter optimization. Try to use exp() and log() instead of the power symbol.");
265        var intPower = Math.Truncate(powerNode.Value);
266        if (intPower != powerNode.Value)
267          throw new NotSupportedException("Only integer powers are allowed in parameter optimization. Try to use exp() and log() instead of the power symbol.");
268        return AutoDiff.TermBuilder.Power(ConvertToAutoDiff(node.GetSubtree(0)), intPower);
269      }
270      if (node.Symbol is Sine) {
271        return sin(
272          ConvertToAutoDiff(node.GetSubtree(0)));
273      }
274      if (node.Symbol is Cosine) {
275        return cos(
276          ConvertToAutoDiff(node.GetSubtree(0)));
277      }
278      if (node.Symbol is Tangent) {
279        return tan(
280          ConvertToAutoDiff(node.GetSubtree(0)));
281      }
282      if (node.Symbol is HyperbolicTangent) {
283        return tanh(
284          ConvertToAutoDiff(node.GetSubtree(0)));
285      }
286      if (node.Symbol is Erf) {
287        return erf(
288          ConvertToAutoDiff(node.GetSubtree(0)));
289      }
290      if (node.Symbol is Norm) {
291        return norm(
292          ConvertToAutoDiff(node.GetSubtree(0)));
293      }
294      if (node.Symbol is StartSymbol) {
295        if (addLinearScalingTerms) {
296          // scaling variables α, β are given at the beginning of the parameter vector
297          var alpha = new AutoDiff.Variable();
298          var beta = new AutoDiff.Variable();
299          variables.Add(beta);
300          variables.Add(alpha);
301          var t = ConvertToAutoDiff(node.GetSubtree(0));
302          return t * alpha + beta;
303        } else return ConvertToAutoDiff(node.GetSubtree(0));
304      }
305      throw new ConversionException();
306    }
307
308
309    // for each factor variable value we need a parameter which represents a binary indicator for that variable & value combination
310    // each binary indicator is only necessary once. So we only create a parameter if this combination is not yet available
311    private static Term FindOrCreateParameter(Dictionary<DataForVariable, AutoDiff.Variable> parameters,
312      string varName, string varValue = "", int lag = 0) {
313      var data = new DataForVariable(varName, varValue, lag);
314
315      AutoDiff.Variable par = null;
316      if (!parameters.TryGetValue(data, out par)) {
317        // not found -> create new parameter and entries in names and values lists
318        par = new AutoDiff.Variable();
319        parameters.Add(data, par);
320      }
321      return par;
322    }
323
324    public static bool IsCompatible(ISymbolicExpressionTree tree) {
325      var containsUnknownSymbol = (
326        from n in tree.Root.GetSubtree(0).IterateNodesPrefix()
327        where
328          !(n.Symbol is Variable) &&
329          !(n.Symbol is BinaryFactorVariable) &&
330          !(n.Symbol is FactorVariable) &&
331          !(n.Symbol is LaggedVariable) &&
332          !(n.Symbol is Constant) &&
333          !(n.Symbol is Addition) &&
334          !(n.Symbol is Subtraction) &&
335          !(n.Symbol is Multiplication) &&
336          !(n.Symbol is Division) &&
337          !(n.Symbol is Logarithm) &&
338          !(n.Symbol is Exponential) &&
339          !(n.Symbol is SquareRoot) &&
340          !(n.Symbol is Square) &&
341          !(n.Symbol is Sine) &&
342          !(n.Symbol is Cosine) &&
343          !(n.Symbol is Tangent) &&
344          !(n.Symbol is HyperbolicTangent) &&
345          !(n.Symbol is Erf) &&
346          !(n.Symbol is Norm) &&
347          !(n.Symbol is StartSymbol) &&
348          !(n.Symbol is Absolute) &&
349          !(n.Symbol is AnalyticQuotient) &&
350          !(n.Symbol is Cube) &&
351          !(n.Symbol is CubeRoot) &&
352          !(n.Symbol is Power)
353        select n).Any();
354      return !containsUnknownSymbol;
355    }
356    #region exception class
357    [Serializable]
358    public class ConversionException : Exception {
359
360      public ConversionException() {
361      }
362
363      public ConversionException(string message) : base(message) {
364      }
365
366      public ConversionException(string message, Exception inner) : base(message, inner) {
367      }
368
369      protected ConversionException(
370        SerializationInfo info,
371        StreamingContext context) : base(info, context) {
372      }
373    }
374    #endregion
375  }
376}
Note: See TracBrowser for help on using the repository browser.