source: branches/GP.Symbols (TimeLag, Diff, Integral)/HeuristicLab.Problems.DataAnalysis/3.3/Symbolic/SimpleArithmeticExpressionInterpreter.cs @ 5060

Last change on this file since 5060 was 5060, checked in by mkommend, 11 years ago

Implemented VariableCondition symbol (ticket #1325).

File size: 14.1 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2010 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.Common;
26using HeuristicLab.Core;
27using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
28using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding.Compiler;
29using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding.Symbols;
30using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
31using HeuristicLab.Problems.DataAnalysis.Symbolic.Symbols;
32
33namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
34  [StorableClass]
35  [Item("SimpleArithmeticExpressionInterpreter", "Interpreter for arithmetic symbolic expression trees including function calls.")]
36  // not thread safe!
37  public sealed class SimpleArithmeticExpressionInterpreter : NamedItem, ISymbolicExpressionTreeInterpreter {
38    private class OpCodes {
39      public const byte Add = 1;
40      public const byte Sub = 2;
41      public const byte Mul = 3;
42      public const byte Div = 4;
43
44      public const byte Sin = 5;
45      public const byte Cos = 6;
46      public const byte Tan = 7;
47
48      public const byte Log = 8;
49      public const byte Exp = 9;
50
51      public const byte IfThenElse = 10;
52
53      public const byte GT = 11;
54      public const byte LT = 12;
55
56      public const byte AND = 13;
57      public const byte OR = 14;
58      public const byte NOT = 15;
59
60
61      public const byte Average = 16;
62
63      public const byte Call = 17;
64
65      public const byte Variable = 18;
66      public const byte LagVariable = 19;
67      public const byte Constant = 20;
68      public const byte Arg = 21;
69
70      public const byte TimeLag = 22;
71      public const byte Integral = 23;
72      public const byte Derivative = 24;
73
74      public const byte VariableCondition = 25;
75    }
76
77    private Dictionary<Type, byte> symbolToOpcode = new Dictionary<Type, byte>() {
78      { typeof(Addition), OpCodes.Add },
79      { typeof(Subtraction), OpCodes.Sub },
80      { typeof(Multiplication), OpCodes.Mul },
81      { typeof(Division), OpCodes.Div },
82      { typeof(Sine), OpCodes.Sin },
83      { typeof(Cosine), OpCodes.Cos },
84      { typeof(Tangent), OpCodes.Tan },
85      { typeof(Logarithm), OpCodes.Log },
86      { typeof(Exponential), OpCodes.Exp },
87      { typeof(IfThenElse), OpCodes.IfThenElse },
88      { typeof(GreaterThan), OpCodes.GT },
89      { typeof(LessThan), OpCodes.LT },
90      { typeof(And), OpCodes.AND },
91      { typeof(Or), OpCodes.OR },
92      { typeof(Not), OpCodes.NOT},
93      { typeof(Average), OpCodes.Average},
94      { typeof(InvokeFunction), OpCodes.Call },
95      { typeof(HeuristicLab.Problems.DataAnalysis.Symbolic.Symbols.Variable), OpCodes.Variable },
96      { typeof(LaggedVariable), OpCodes.LagVariable },
97      { typeof(Constant), OpCodes.Constant },
98      { typeof(Argument), OpCodes.Arg },
99      { typeof(TimeLag), OpCodes.TimeLag},
100      { typeof(Integral), OpCodes.Integral},
101      { typeof(Derivative), OpCodes.Derivative},
102      { typeof(VariableCondition),OpCodes.VariableCondition}
103    };
104    private const int ARGUMENT_STACK_SIZE = 1024;
105
106    private Dataset dataset;
107    private int row;
108    private Instruction[] code;
109    private int pc;
110    private double[] argumentStack = new double[ARGUMENT_STACK_SIZE];
111    private int argStackPointer;
112
113    public override bool CanChangeName {
114      get { return false; }
115    }
116    public override bool CanChangeDescription {
117      get { return false; }
118    }
119
120    [StorableConstructor]
121    private SimpleArithmeticExpressionInterpreter(bool deserializing) : base(deserializing) { }
122    private SimpleArithmeticExpressionInterpreter(SimpleArithmeticExpressionInterpreter original, Cloner cloner) : base(original, cloner) { }
123
124    public override IDeepCloneable Clone(Cloner cloner) {
125      return new SimpleArithmeticExpressionInterpreter(this, cloner);
126    }
127
128    public SimpleArithmeticExpressionInterpreter()
129      : base() {
130    }
131
132    public IEnumerable<double> GetSymbolicExpressionTreeValues(SymbolicExpressionTree tree, Dataset dataset, IEnumerable<int> rows) {
133      this.dataset = dataset;
134      var compiler = new SymbolicExpressionTreeCompiler();
135      compiler.AddInstructionPostProcessingHook(PostProcessInstruction);
136      code = compiler.Compile(tree, MapSymbolToOpCode);
137      foreach (var row in rows) {
138        this.row = row;
139        pc = 0;
140        argStackPointer = 0;
141        yield return Evaluate();
142      }
143    }
144
145    private Instruction PostProcessInstruction(Instruction instr) {
146      if (instr.opCode == OpCodes.Variable) {
147        var variableTreeNode = instr.dynamicNode as VariableTreeNode;
148        instr.iArg0 = (ushort)dataset.GetVariableIndex(variableTreeNode.VariableName);
149      } else if (instr.opCode == OpCodes.LagVariable) {
150        var variableTreeNode = instr.dynamicNode as LaggedVariableTreeNode;
151        instr.iArg0 = (ushort)dataset.GetVariableIndex(variableTreeNode.VariableName);
152      } else if (instr.opCode == OpCodes.VariableCondition) {
153        var variableConditionTreeNode = instr.dynamicNode as VariableConditionTreeNode;
154        instr.iArg0 = (ushort)dataset.GetVariableIndex(variableConditionTreeNode.VariableName);
155      }
156      return instr;
157    }
158
159    private byte MapSymbolToOpCode(SymbolicExpressionTreeNode treeNode) {
160      if (symbolToOpcode.ContainsKey(treeNode.Symbol.GetType()))
161        return symbolToOpcode[treeNode.Symbol.GetType()];
162      else
163        throw new NotSupportedException("Symbol: " + treeNode.Symbol);
164    }
165
166    private double Evaluate() {
167      Instruction currentInstr = code[pc++];
168      switch (currentInstr.opCode) {
169        case OpCodes.Add: {
170            double s = Evaluate();
171            for (int i = 1; i < currentInstr.nArguments; i++) {
172              s += Evaluate();
173            }
174            return s;
175          }
176        case OpCodes.Sub: {
177            double s = Evaluate();
178            for (int i = 1; i < currentInstr.nArguments; i++) {
179              s -= Evaluate();
180            }
181            if (currentInstr.nArguments == 1) s = -s;
182            return s;
183          }
184        case OpCodes.Mul: {
185            double p = Evaluate();
186            for (int i = 1; i < currentInstr.nArguments; i++) {
187              p *= Evaluate();
188            }
189            return p;
190          }
191        case OpCodes.Div: {
192            double p = Evaluate();
193            for (int i = 1; i < currentInstr.nArguments; i++) {
194              p /= Evaluate();
195            }
196            if (currentInstr.nArguments == 1) p = 1.0 / p;
197            return p;
198          }
199        case OpCodes.Average: {
200            double sum = Evaluate();
201            for (int i = 1; i < currentInstr.nArguments; i++) {
202              sum += Evaluate();
203            }
204            return sum / currentInstr.nArguments;
205          }
206        case OpCodes.Cos: {
207            return Math.Cos(Evaluate());
208          }
209        case OpCodes.Sin: {
210            return Math.Sin(Evaluate());
211          }
212        case OpCodes.Tan: {
213            return Math.Tan(Evaluate());
214          }
215        case OpCodes.Exp: {
216            return Math.Exp(Evaluate());
217          }
218        case OpCodes.Log: {
219            return Math.Log(Evaluate());
220          }
221        case OpCodes.IfThenElse: {
222            double condition = Evaluate();
223            double result;
224            if (condition > 0.0) {
225              result = Evaluate(); SkipBakedCode();
226            } else {
227              SkipBakedCode(); result = Evaluate();
228            }
229            return result;
230          }
231        case OpCodes.AND: {
232            double result = Evaluate();
233            for (int i = 1; i < currentInstr.nArguments; i++) {
234              if (result <= 0.0) SkipBakedCode();
235              else {
236                result = Evaluate();
237              }
238            }
239            return result <= 0.0 ? -1.0 : 1.0;
240          }
241        case OpCodes.OR: {
242            double result = Evaluate();
243            for (int i = 1; i < currentInstr.nArguments; i++) {
244              if (result > 0.0) SkipBakedCode();
245              else {
246                result = Evaluate();
247              }
248            }
249            return result > 0.0 ? 1.0 : -1.0;
250          }
251        case OpCodes.NOT: {
252            return -Evaluate();
253          }
254        case OpCodes.GT: {
255            double x = Evaluate();
256            double y = Evaluate();
257            if (x > y) return 1.0;
258            else return -1.0;
259          }
260        case OpCodes.LT: {
261            double x = Evaluate();
262            double y = Evaluate();
263            if (x < y) return 1.0;
264            else return -1.0;
265          }
266        case OpCodes.Call: {
267            // evaluate sub-trees
268            // push on argStack in reverse order
269            for (int i = 0; i < currentInstr.nArguments; i++) {
270              argumentStack[argStackPointer + currentInstr.nArguments - i] = Evaluate();
271            }
272            argStackPointer += currentInstr.nArguments;
273
274            // save the pc
275            int nextPc = pc;
276            // set pc to start of function 
277            pc = currentInstr.iArg0;
278            // evaluate the function
279            double v = Evaluate();
280
281            // decrease the argument stack pointer by the number of arguments pushed
282            // to set the argStackPointer back to the original location
283            argStackPointer -= currentInstr.nArguments;
284
285            // restore the pc => evaluation will continue at point after my subtrees 
286            pc = nextPc;
287            return v;
288          }
289        case OpCodes.Arg: {
290            return argumentStack[argStackPointer - currentInstr.iArg0];
291          }
292        case OpCodes.Variable: {
293            var variableTreeNode = currentInstr.dynamicNode as VariableTreeNode;
294            return dataset[row, currentInstr.iArg0] * variableTreeNode.Weight;
295          }
296        case OpCodes.LagVariable: {
297            var lagVariableTreeNode = currentInstr.dynamicNode as LaggedVariableTreeNode;
298            int actualRow = row + lagVariableTreeNode.Lag;
299            if (actualRow < 0 || actualRow >= dataset.Rows) throw new ArgumentException("Out of range access to dataset row: " + row);
300            return dataset[actualRow, currentInstr.iArg0] * lagVariableTreeNode.Weight;
301          }
302        case OpCodes.Constant: {
303            var constTreeNode = currentInstr.dynamicNode as ConstantTreeNode;
304            return constTreeNode.Value;
305          }
306        case OpCodes.TimeLag: {
307            var timeLagTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
308            if (row + timeLagTreeNode.Lag < 0 || row + timeLagTreeNode.Lag >= dataset.Rows)
309              return double.NaN;
310
311            row += timeLagTreeNode.Lag;
312            double result = Evaluate();
313            row -= timeLagTreeNode.Lag;
314            return result;
315          }
316        case OpCodes.Integral: {
317            int nextPc = pc;
318            var timeLagTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
319            if (row + timeLagTreeNode.Lag < 0 || row + timeLagTreeNode.Lag >= dataset.Rows)
320              return double.NaN;
321            double sum = 0.0;
322            if (timeLagTreeNode.IterateNodesPrefix().OfType<VariableTreeNode>().Any()) {
323              for (int i = 0; i < Math.Abs(timeLagTreeNode.Lag); i++) {
324                row += Math.Sign(timeLagTreeNode.Lag);
325                sum += Evaluate();
326                pc = nextPc;
327              }
328              row -= timeLagTreeNode.Lag;
329              sum += Evaluate();
330            } else sum = Math.Abs(timeLagTreeNode.Lag) * Evaluate();
331            return sum;
332          }
333
334        //mkommend: derivate calculation taken from:
335        //http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
336        //one sided smooth differentiatior, N = 4
337        // y' = 1/8h (f_i + 2f_i-1, -2 f_i-3 - f_i-4)
338        case OpCodes.Derivative: {
339            if (row - 4 < 0) return double.NaN;
340            int nextPc = pc;
341            double f_0 = Evaluate(); row--;
342            pc = nextPc;
343            double f_1 = Evaluate(); row -= 2;
344            pc = nextPc;
345            double f_3 = Evaluate(); row--;
346            pc = nextPc;
347            double f_4 = Evaluate();
348            row += 4;
349
350            return (f_0 + 2 * f_1 - 2 * f_3 - f_4) / 8; // h = 1
351          }
352
353        //mkommend: this symbol uses the logistic function f(x) = 1 / (1 + e^(-alpha * x) )
354        //to determine the relative amounts of the true and false branch see http://en.wikipedia.org/wiki/Logistic_function
355        case OpCodes.VariableCondition: {
356            var variableConditionTreeNode = (VariableConditionTreeNode)currentInstr.dynamicNode;
357            double variableValue = dataset[row, currentInstr.iArg0];
358            double x = variableValue - variableConditionTreeNode.Threshold;
359            double p = 1 / (1 + Math.Exp(-variableConditionTreeNode.Slope * x));
360
361            double trueBranch = Evaluate();
362            double falseBranch = Evaluate();
363
364            return trueBranch * p + falseBranch * (1 - p);
365          }
366        default: throw new NotSupportedException();
367      }
368    }
369
370    // skips a whole branch
371    private void SkipBakedCode() {
372      int i = 1;
373      while (i > 0) {
374        i += code[pc++].nArguments;
375        i--;
376      }
377    }
378  }
379}
Note: See TracBrowser for help on using the repository browser.