Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/SymbolicDataAnalysisExpressionTreeInterpreter.cs @ 7695

Last change on this file since 7695 was 7695, checked in by gkronber, 12 years ago

#1810 merged patch to add square and square root function symbols by mkommend
File size: 19.8 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2012 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 HeuristicLab.Common;
25using HeuristicLab.Core;
26using HeuristicLab.Data;
27using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
28using HeuristicLab.Parameters;
29using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
30
31namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
32  [StorableClass]
33  [Item("SymbolicDataAnalysisExpressionTreeInterpreter", "Interpreter for symbolic expression trees including automatically defined functions.")]
34  public sealed class SymbolicDataAnalysisExpressionTreeInterpreter : ParameterizedNamedItem, ISymbolicDataAnalysisExpressionTreeInterpreter {
35    private const string CheckExpressionsWithIntervalArithmeticParameterName = "CheckExpressionsWithIntervalArithmetic";
36    private const string EvaluatedSolutionsParameterName = "EvaluatedSolutions";
37    #region private classes
38    private class InterpreterState {
39      private double[] argumentStack;
40      private int argumentStackPointer;
41      private Instruction[] code;
42      private int pc;
43      public int ProgramCounter {
44        get { return pc; }
45        set { pc = value; }
46      }
47      internal InterpreterState(Instruction[] code, int argumentStackSize) {
48        this.code = code;
49        this.pc = 0;
50        if (argumentStackSize > 0) {
51          this.argumentStack = new double[argumentStackSize];
52        }
53        this.argumentStackPointer = 0;
54      }
55
56      internal void Reset() {
57        this.pc = 0;
58        this.argumentStackPointer = 0;
59      }
60
61      internal Instruction NextInstruction() {
62        return code[pc++];
63      }
64      private void Push(double val) {
65        argumentStack[argumentStackPointer++] = val;
66      }
67      private double Pop() {
68        return argumentStack[--argumentStackPointer];
69      }
70
71      internal void CreateStackFrame(double[] argValues) {
72        // push in reverse order to make indexing easier
73        for (int i = argValues.Length - 1; i >= 0; i--) {
74          argumentStack[argumentStackPointer++] = argValues[i];
75        }
76        Push(argValues.Length);
77      }
78
79      internal void RemoveStackFrame() {
80        int size = (int)Pop();
81        argumentStackPointer -= size;
82      }
83
84      internal double GetStackFrameValue(ushort index) {
85        // layout of stack:
86        // [0]   <- argumentStackPointer
87        // [StackFrameSize = N + 1]
88        // [Arg0] <- argumentStackPointer - 2 - 0
89        // [Arg1] <- argumentStackPointer - 2 - 1
90        // [...]
91        // [ArgN] <- argumentStackPointer - 2 - N
92        // <Begin of stack frame>
93        return argumentStack[argumentStackPointer - index - 2];
94      }
95    }
96    private class OpCodes {
97      public const byte Add = 1;
98      public const byte Sub = 2;
99      public const byte Mul = 3;
100      public const byte Div = 4;
101
102      public const byte Sin = 5;
103      public const byte Cos = 6;
104      public const byte Tan = 7;
105
106      public const byte Log = 8;
107      public const byte Exp = 9;
108
109      public const byte IfThenElse = 10;
110
111      public const byte GT = 11;
112      public const byte LT = 12;
113
114      public const byte AND = 13;
115      public const byte OR = 14;
116      public const byte NOT = 15;
117
118
119      public const byte Average = 16;
120
121      public const byte Call = 17;
122
123      public const byte Variable = 18;
124      public const byte LagVariable = 19;
125      public const byte Constant = 20;
126      public const byte Arg = 21;
127
128      public const byte Power = 22;
129      public const byte Root = 23;
130      public const byte TimeLag = 24;
131      public const byte Integral = 25;
132      public const byte Derivative = 26;
133
134      public const byte VariableCondition = 27;
135
136      public const byte Square = 28;
137      public const byte SquareRoot = 29;
138    }
139    #endregion
140
141    private Dictionary<Type, byte> symbolToOpcode = new Dictionary<Type, byte>() {
142      { typeof(Addition), OpCodes.Add },
143      { typeof(Subtraction), OpCodes.Sub },
144      { typeof(Multiplication), OpCodes.Mul },
145      { typeof(Division), OpCodes.Div },
146      { typeof(Sine), OpCodes.Sin },
147      { typeof(Cosine), OpCodes.Cos },
148      { typeof(Tangent), OpCodes.Tan },
149      { typeof(Logarithm), OpCodes.Log },
150      { typeof(Exponential), OpCodes.Exp },
151      { typeof(IfThenElse), OpCodes.IfThenElse },
152      { typeof(GreaterThan), OpCodes.GT },
153      { typeof(LessThan), OpCodes.LT },
154      { typeof(And), OpCodes.AND },
155      { typeof(Or), OpCodes.OR },
156      { typeof(Not), OpCodes.NOT},
157      { typeof(Average), OpCodes.Average},
158      { typeof(InvokeFunction), OpCodes.Call },
159      { typeof(Variable), OpCodes.Variable },
160      { typeof(LaggedVariable), OpCodes.LagVariable },
161      { typeof(Constant), OpCodes.Constant },
162      { typeof(Argument), OpCodes.Arg },
163      { typeof(Power),OpCodes.Power},
164      { typeof(Root),OpCodes.Root},
165      { typeof(TimeLag), OpCodes.TimeLag},
166      { typeof(Integral), OpCodes.Integral},
167      { typeof(Derivative), OpCodes.Derivative},
168      { typeof(VariableCondition),OpCodes.VariableCondition},
169      { typeof(Square),OpCodes.Square},
170      {typeof(SquareRoot),OpCodes.SquareRoot}
171    };
172
173    public override bool CanChangeName {
174      get { return false; }
175    }
176    public override bool CanChangeDescription {
177      get { return false; }
178    }
179
180    #region parameter properties
181    public IValueParameter<BoolValue> CheckExpressionsWithIntervalArithmeticParameter {
182      get { return (IValueParameter<BoolValue>)Parameters[CheckExpressionsWithIntervalArithmeticParameterName]; }
183    }
184
185    public IValueParameter<IntValue> EvaluatedSolutionsParameter {
186      get { return (IValueParameter<IntValue>)Parameters[EvaluatedSolutionsParameterName]; }
187    }
188    #endregion
189
190    #region properties
191    public BoolValue CheckExpressionsWithIntervalArithmetic {
192      get { return CheckExpressionsWithIntervalArithmeticParameter.Value; }
193      set { CheckExpressionsWithIntervalArithmeticParameter.Value = value; }
194    }
195
196    public IntValue EvaluatedSolutions {
197      get { return EvaluatedSolutionsParameter.Value; }
198      set { EvaluatedSolutionsParameter.Value = value; }
199    }
200    #endregion
201
202    [StorableConstructor]
203    private SymbolicDataAnalysisExpressionTreeInterpreter(bool deserializing) : base(deserializing) { }
204    private SymbolicDataAnalysisExpressionTreeInterpreter(SymbolicDataAnalysisExpressionTreeInterpreter original, Cloner cloner) : base(original, cloner) { }
205    public override IDeepCloneable Clone(Cloner cloner) {
206      return new SymbolicDataAnalysisExpressionTreeInterpreter(this, cloner);
207    }
208
209    public SymbolicDataAnalysisExpressionTreeInterpreter()
210      : base("SymbolicDataAnalysisExpressionTreeInterpreter", "Interpreter for symbolic expression trees including automatically defined functions.") {
211      Parameters.Add(new ValueParameter<BoolValue>(CheckExpressionsWithIntervalArithmeticParameterName, "Switch that determines if the interpreter checks the validity of expressions with interval arithmetic before evaluating the expression.", new BoolValue(false)));
212      Parameters.Add(new ValueParameter<IntValue>(EvaluatedSolutionsParameterName, "A counter for the total number of solutions the interpreter has evaluated", new IntValue(0)));
213    }
214
215    [StorableHook(HookType.AfterDeserialization)]
216    private void AfterDeserialization() {
217      if (!Parameters.ContainsKey(EvaluatedSolutionsParameterName))
218        Parameters.Add(new ValueParameter<IntValue>(EvaluatedSolutionsParameterName, "A counter for the total number of solutions the interpreter has evaluated", new IntValue(0)));
219    }
220
221    #region IStatefulItem
222    public void InitializeState() {
223      EvaluatedSolutions.Value = 0;
224    }
225
226    public void ClearState() {
227    }
228    #endregion
229
230    public IEnumerable<double> GetSymbolicExpressionTreeValues(ISymbolicExpressionTree tree, Dataset dataset, IEnumerable<int> rows) {
231      if (CheckExpressionsWithIntervalArithmetic.Value)
232        throw new NotSupportedException("Interval arithmetic is not yet supported in the symbolic data analysis interpreter.");
233      EvaluatedSolutions.Value++; // increment the evaluated solutions counter
234      var compiler = new SymbolicExpressionTreeCompiler();
235      Instruction[] code = compiler.Compile(tree, MapSymbolToOpCode);
236      int necessaryArgStackSize = 0;
237      for (int i = 0; i < code.Length; i++) {
238        Instruction instr = code[i];
239        if (instr.opCode == OpCodes.Variable) {
240          var variableTreeNode = instr.dynamicNode as VariableTreeNode;
241          instr.iArg0 = dataset.GetReadOnlyDoubleValues(variableTreeNode.VariableName);
242          code[i] = instr;
243        } else if (instr.opCode == OpCodes.LagVariable) {
244          var laggedVariableTreeNode = instr.dynamicNode as LaggedVariableTreeNode;
245          instr.iArg0 = dataset.GetReadOnlyDoubleValues(laggedVariableTreeNode.VariableName);
246          code[i] = instr;
247        } else if (instr.opCode == OpCodes.VariableCondition) {
248          var variableConditionTreeNode = instr.dynamicNode as VariableConditionTreeNode;
249          instr.iArg0 = dataset.GetReadOnlyDoubleValues(variableConditionTreeNode.VariableName);
250        } else if (instr.opCode == OpCodes.Call) {
251          necessaryArgStackSize += instr.nArguments + 1;
252        }
253      }
254      var state = new InterpreterState(code, necessaryArgStackSize);
255
256      foreach (var rowEnum in rows) {
257        int row = rowEnum;
258        state.Reset();
259        yield return Evaluate(dataset, ref row, state);
260      }
261    }
262
263    private double Evaluate(Dataset dataset, ref int row, InterpreterState state) {
264      Instruction currentInstr = state.NextInstruction();
265      switch (currentInstr.opCode) {
266        case OpCodes.Add: {
267            double s = Evaluate(dataset, ref row, state);
268            for (int i = 1; i < currentInstr.nArguments; i++) {
269              s += Evaluate(dataset, ref row, state);
270            }
271            return s;
272          }
273        case OpCodes.Sub: {
274            double s = Evaluate(dataset, ref row, state);
275            for (int i = 1; i < currentInstr.nArguments; i++) {
276              s -= Evaluate(dataset, ref row, state);
277            }
278            if (currentInstr.nArguments == 1) s = -s;
279            return s;
280          }
281        case OpCodes.Mul: {
282            double p = Evaluate(dataset, ref row, state);
283            for (int i = 1; i < currentInstr.nArguments; i++) {
284              p *= Evaluate(dataset, ref row, state);
285            }
286            return p;
287          }
288        case OpCodes.Div: {
289            double p = Evaluate(dataset, ref row, state);
290            for (int i = 1; i < currentInstr.nArguments; i++) {
291              p /= Evaluate(dataset, ref row, state);
292            }
293            if (currentInstr.nArguments == 1) p = 1.0 / p;
294            return p;
295          }
296        case OpCodes.Average: {
297            double sum = Evaluate(dataset, ref row, state);
298            for (int i = 1; i < currentInstr.nArguments; i++) {
299              sum += Evaluate(dataset, ref row, state);
300            }
301            return sum / currentInstr.nArguments;
302          }
303        case OpCodes.Cos: {
304            return Math.Cos(Evaluate(dataset, ref row, state));
305          }
306        case OpCodes.Sin: {
307            return Math.Sin(Evaluate(dataset, ref row, state));
308          }
309        case OpCodes.Tan: {
310            return Math.Tan(Evaluate(dataset, ref row, state));
311          }
312        case OpCodes.Square: {
313            return Math.Pow(Evaluate(dataset, ref row, state), 2);
314          }
315        case OpCodes.Power: {
316            double x = Evaluate(dataset, ref row, state);
317            double y = Math.Round(Evaluate(dataset, ref row, state));
318            return Math.Pow(x, y);
319          }
320        case OpCodes.SquareRoot: {
321            return Math.Sqrt(Evaluate(dataset, ref row, state));
322          }
323        case OpCodes.Root: {
324            double x = Evaluate(dataset, ref row, state);
325            double y = Math.Round(Evaluate(dataset, ref row, state));
326            return Math.Pow(x, 1 / y);
327          }
328        case OpCodes.Exp: {
329            return Math.Exp(Evaluate(dataset, ref row, state));
330          }
331        case OpCodes.Log: {
332            return Math.Log(Evaluate(dataset, ref row, state));
333          }
334        case OpCodes.IfThenElse: {
335            double condition = Evaluate(dataset, ref row, state);
336            double result;
337            if (condition > 0.0) {
338              result = Evaluate(dataset, ref row, state); SkipInstructions(state);
339            } else {
340              SkipInstructions(state); result = Evaluate(dataset, ref row, state);
341            }
342            return result;
343          }
344        case OpCodes.AND: {
345            double result = Evaluate(dataset, ref row, state);
346            for (int i = 1; i < currentInstr.nArguments; i++) {
347              if (result > 0.0) result = Evaluate(dataset, ref row, state);
348              else {
349                SkipInstructions(state);
350              }
351            }
352            return result > 0.0 ? 1.0 : -1.0;
353          }
354        case OpCodes.OR: {
355            double result = Evaluate(dataset, ref row, state);
356            for (int i = 1; i < currentInstr.nArguments; i++) {
357              if (result <= 0.0) result = Evaluate(dataset, ref row, state);
358              else {
359                SkipInstructions(state);
360              }
361            }
362            return result > 0.0 ? 1.0 : -1.0;
363          }
364        case OpCodes.NOT: {
365            return Evaluate(dataset, ref row, state) > 0.0 ? -1.0 : 1.0;
366          }
367        case OpCodes.GT: {
368            double x = Evaluate(dataset, ref row, state);
369            double y = Evaluate(dataset, ref row, state);
370            if (x > y) return 1.0;
371            else return -1.0;
372          }
373        case OpCodes.LT: {
374            double x = Evaluate(dataset, ref row, state);
375            double y = Evaluate(dataset, ref row, state);
376            if (x < y) return 1.0;
377            else return -1.0;
378          }
379        case OpCodes.TimeLag: {
380            var timeLagTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
381            row += timeLagTreeNode.Lag;
382            double result = Evaluate(dataset, ref row, state);
383            row -= timeLagTreeNode.Lag;
384            return result;
385          }
386        case OpCodes.Integral: {
387            int savedPc = state.ProgramCounter;
388            var timeLagTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
389            double sum = 0.0;
390            for (int i = 0; i < Math.Abs(timeLagTreeNode.Lag); i++) {
391              row += Math.Sign(timeLagTreeNode.Lag);
392              sum += Evaluate(dataset, ref row, state);
393              state.ProgramCounter = savedPc;
394            }
395            row -= timeLagTreeNode.Lag;
396            sum += Evaluate(dataset, ref row, state);
397            return sum;
398          }
399
400        //mkommend: derivate calculation taken from:
401        //http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
402        //one sided smooth differentiatior, N = 4
403        // y' = 1/8h (f_i + 2f_i-1, -2 f_i-3 - f_i-4)
404        case OpCodes.Derivative: {
405            int savedPc = state.ProgramCounter;
406            double f_0 = Evaluate(dataset, ref row, state); row--;
407            state.ProgramCounter = savedPc;
408            double f_1 = Evaluate(dataset, ref row, state); row -= 2;
409            state.ProgramCounter = savedPc;
410            double f_3 = Evaluate(dataset, ref row, state); row--;
411            state.ProgramCounter = savedPc;
412            double f_4 = Evaluate(dataset, ref row, state);
413            row += 4;
414
415            return (f_0 + 2 * f_1 - 2 * f_3 - f_4) / 8; // h = 1
416          }
417        case OpCodes.Call: {
418            // evaluate sub-trees
419            double[] argValues = new double[currentInstr.nArguments];
420            for (int i = 0; i < currentInstr.nArguments; i++) {
421              argValues[i] = Evaluate(dataset, ref row, state);
422            }
423            // push on argument values on stack
424            state.CreateStackFrame(argValues);
425
426            // save the pc
427            int savedPc = state.ProgramCounter;
428            // set pc to start of function 
429            state.ProgramCounter = (ushort)currentInstr.iArg0;
430            // evaluate the function
431            double v = Evaluate(dataset, ref row, state);
432
433            // delete the stack frame
434            state.RemoveStackFrame();
435
436            // restore the pc => evaluation will continue at point after my subtrees 
437            state.ProgramCounter = savedPc;
438            return v;
439          }
440        case OpCodes.Arg: {
441            return state.GetStackFrameValue((ushort)currentInstr.iArg0);
442          }
443        case OpCodes.Variable: {
444            if (row < 0 || row >= dataset.Rows)
445              return double.NaN;
446            var variableTreeNode = (VariableTreeNode)currentInstr.dynamicNode;
447            return ((IList<double>)currentInstr.iArg0)[row] * variableTreeNode.Weight;
448          }
449        case OpCodes.LagVariable: {
450            var laggedVariableTreeNode = (LaggedVariableTreeNode)currentInstr.dynamicNode;
451            int actualRow = row + laggedVariableTreeNode.Lag;
452            if (actualRow < 0 || actualRow >= dataset.Rows)
453              return double.NaN;
454            return ((IList<double>)currentInstr.iArg0)[actualRow] * laggedVariableTreeNode.Weight;
455          }
456        case OpCodes.Constant: {
457            var constTreeNode = currentInstr.dynamicNode as ConstantTreeNode;
458            return constTreeNode.Value;
459          }
460
461        //mkommend: this symbol uses the logistic function f(x) = 1 / (1 + e^(-alpha * x) )
462        //to determine the relative amounts of the true and false branch see http://en.wikipedia.org/wiki/Logistic_function
463        case OpCodes.VariableCondition: {
464            if (row < 0 || row >= dataset.Rows)
465              return double.NaN;
466            var variableConditionTreeNode = (VariableConditionTreeNode)currentInstr.dynamicNode;
467            double variableValue = ((IList<double>)currentInstr.iArg0)[row];
468            double x = variableValue - variableConditionTreeNode.Threshold;
469            double p = 1 / (1 + Math.Exp(-variableConditionTreeNode.Slope * x));
470
471            double trueBranch = Evaluate(dataset, ref row, state);
472            double falseBranch = Evaluate(dataset, ref row, state);
473
474            return trueBranch * p + falseBranch * (1 - p);
475          }
476        default: throw new NotSupportedException();
477      }
478    }
479
480    private byte MapSymbolToOpCode(ISymbolicExpressionTreeNode treeNode) {
481      if (symbolToOpcode.ContainsKey(treeNode.Symbol.GetType()))
482        return symbolToOpcode[treeNode.Symbol.GetType()];
483      else
484        throw new NotSupportedException("Symbol: " + treeNode.Symbol);
485    }
486
487    // skips a whole branch
488    private void SkipInstructions(InterpreterState state) {
489      int i = 1;
490      while (i > 0) {
491        i += state.NextInstruction().nArguments;
492        i--;
493      }
494    }
495  }
496}
Note: See TracBrowser for help on using the repository browser.