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

Last change on this file since 5052 was 5052, checked in by mkommend, 12 years ago

Removed unused caching code (ticket #1256).

File size: 13.0 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
75    private Dictionary<Type, byte> symbolToOpcode = new Dictionary<Type, byte>() {
76      { typeof(Addition), OpCodes.Add },
77      { typeof(Subtraction), OpCodes.Sub },
78      { typeof(Multiplication), OpCodes.Mul },
79      { typeof(Division), OpCodes.Div },
80      { typeof(Sine), OpCodes.Sin },
81      { typeof(Cosine), OpCodes.Cos },
82      { typeof(Tangent), OpCodes.Tan },
83      { typeof(Logarithm), OpCodes.Log },
84      { typeof(Exponential), OpCodes.Exp },
85      { typeof(IfThenElse), OpCodes.IfThenElse },
86      { typeof(GreaterThan), OpCodes.GT },
87      { typeof(LessThan), OpCodes.LT },
88      { typeof(And), OpCodes.AND },
89      { typeof(Or), OpCodes.OR },
90      { typeof(Not), OpCodes.NOT},
91      { typeof(Average), OpCodes.Average},
92      { typeof(InvokeFunction), OpCodes.Call },
93      { typeof(HeuristicLab.Problems.DataAnalysis.Symbolic.Symbols.Variable), OpCodes.Variable },
94      { typeof(LaggedVariable), OpCodes.LagVariable },
95      { typeof(Constant), OpCodes.Constant },
96      { typeof(Argument), OpCodes.Arg },
97      { typeof(TimeLag), OpCodes.TimeLag},
98      { typeof(Integral), OpCodes.Integral},
99      { typeof(Derivative), OpCodes.Derivative}
100    };
101    private const int ARGUMENT_STACK_SIZE = 1024;
102
103    private Dataset dataset;
104    private int row;
105    private Instruction[] code;
106    private int pc;
107    private double[] argumentStack = new double[ARGUMENT_STACK_SIZE];
108    private int argStackPointer;
109
110    public override bool CanChangeName {
111      get { return false; }
112    }
113    public override bool CanChangeDescription {
114      get { return false; }
115    }
116
117    [StorableConstructor]
118    private SimpleArithmeticExpressionInterpreter(bool deserializing) : base(deserializing) { }
119    private SimpleArithmeticExpressionInterpreter(SimpleArithmeticExpressionInterpreter original, Cloner cloner) : base(original, cloner) { }
120
121    public override IDeepCloneable Clone(Cloner cloner) {
122      return new SimpleArithmeticExpressionInterpreter(this, cloner);
123    }
124
125    public SimpleArithmeticExpressionInterpreter()
126      : base() {
127    }
128
129    public IEnumerable<double> GetSymbolicExpressionTreeValues(SymbolicExpressionTree tree, Dataset dataset, IEnumerable<int> rows) {
130      this.dataset = dataset;
131      var compiler = new SymbolicExpressionTreeCompiler();
132      compiler.AddInstructionPostProcessingHook(PostProcessInstruction);
133      code = compiler.Compile(tree, MapSymbolToOpCode);
134      foreach (var row in rows) {
135        this.row = row;
136        pc = 0;
137        argStackPointer = 0;
138        yield return Evaluate();
139      }
140    }
141
142    private Instruction PostProcessInstruction(Instruction instr) {
143      if (instr.opCode == OpCodes.Variable) {
144        var variableTreeNode = instr.dynamicNode as VariableTreeNode;
145        instr.iArg0 = (ushort)dataset.GetVariableIndex(variableTreeNode.VariableName);
146      } else if (instr.opCode == OpCodes.LagVariable) {
147        var variableTreeNode = instr.dynamicNode as LaggedVariableTreeNode;
148        instr.iArg0 = (ushort)dataset.GetVariableIndex(variableTreeNode.VariableName);
149      }
150      return instr;
151    }
152
153    private byte MapSymbolToOpCode(SymbolicExpressionTreeNode treeNode) {
154      if (symbolToOpcode.ContainsKey(treeNode.Symbol.GetType()))
155        return symbolToOpcode[treeNode.Symbol.GetType()];
156      else
157        throw new NotSupportedException("Symbol: " + treeNode.Symbol);
158    }
159
160    private double Evaluate() {
161      Instruction currentInstr = code[pc++];
162      switch (currentInstr.opCode) {
163        case OpCodes.Add: {
164            double s = Evaluate();
165            for (int i = 1; i < currentInstr.nArguments; i++) {
166              s += Evaluate();
167            }
168            return s;
169          }
170        case OpCodes.Sub: {
171            double s = Evaluate();
172            for (int i = 1; i < currentInstr.nArguments; i++) {
173              s -= Evaluate();
174            }
175            if (currentInstr.nArguments == 1) s = -s;
176            return s;
177          }
178        case OpCodes.Mul: {
179            double p = Evaluate();
180            for (int i = 1; i < currentInstr.nArguments; i++) {
181              p *= Evaluate();
182            }
183            return p;
184          }
185        case OpCodes.Div: {
186            double p = Evaluate();
187            for (int i = 1; i < currentInstr.nArguments; i++) {
188              p /= Evaluate();
189            }
190            if (currentInstr.nArguments == 1) p = 1.0 / p;
191            return p;
192          }
193        case OpCodes.Average: {
194            double sum = Evaluate();
195            for (int i = 1; i < currentInstr.nArguments; i++) {
196              sum += Evaluate();
197            }
198            return sum / currentInstr.nArguments;
199          }
200        case OpCodes.Cos: {
201            return Math.Cos(Evaluate());
202          }
203        case OpCodes.Sin: {
204            return Math.Sin(Evaluate());
205          }
206        case OpCodes.Tan: {
207            return Math.Tan(Evaluate());
208          }
209        case OpCodes.Exp: {
210            return Math.Exp(Evaluate());
211          }
212        case OpCodes.Log: {
213            return Math.Log(Evaluate());
214          }
215        case OpCodes.IfThenElse: {
216            double condition = Evaluate();
217            double result;
218            if (condition > 0.0) {
219              result = Evaluate(); SkipBakedCode();
220            } else {
221              SkipBakedCode(); result = Evaluate();
222            }
223            return result;
224          }
225        case OpCodes.AND: {
226            double result = Evaluate();
227            for (int i = 1; i < currentInstr.nArguments; i++) {
228              if (result <= 0.0) SkipBakedCode();
229              else {
230                result = Evaluate();
231              }
232            }
233            return result <= 0.0 ? -1.0 : 1.0;
234          }
235        case OpCodes.OR: {
236            double result = Evaluate();
237            for (int i = 1; i < currentInstr.nArguments; i++) {
238              if (result > 0.0) SkipBakedCode();
239              else {
240                result = Evaluate();
241              }
242            }
243            return result > 0.0 ? 1.0 : -1.0;
244          }
245        case OpCodes.NOT: {
246            return -Evaluate();
247          }
248        case OpCodes.GT: {
249            double x = Evaluate();
250            double y = Evaluate();
251            if (x > y) return 1.0;
252            else return -1.0;
253          }
254        case OpCodes.LT: {
255            double x = Evaluate();
256            double y = Evaluate();
257            if (x < y) return 1.0;
258            else return -1.0;
259          }
260        case OpCodes.Call: {
261            // evaluate sub-trees
262            // push on argStack in reverse order
263            for (int i = 0; i < currentInstr.nArguments; i++) {
264              argumentStack[argStackPointer + currentInstr.nArguments - i] = Evaluate();
265            }
266            argStackPointer += currentInstr.nArguments;
267
268            // save the pc
269            int nextPc = pc;
270            // set pc to start of function 
271            pc = currentInstr.iArg0;
272            // evaluate the function
273            double v = Evaluate();
274
275            // decrease the argument stack pointer by the number of arguments pushed
276            // to set the argStackPointer back to the original location
277            argStackPointer -= currentInstr.nArguments;
278
279            // restore the pc => evaluation will continue at point after my subtrees 
280            pc = nextPc;
281            return v;
282          }
283        case OpCodes.Arg: {
284            return argumentStack[argStackPointer - currentInstr.iArg0];
285          }
286        case OpCodes.Variable: {
287            var variableTreeNode = currentInstr.dynamicNode as VariableTreeNode;
288            return dataset[row, currentInstr.iArg0] * variableTreeNode.Weight;
289          }
290        case OpCodes.LagVariable: {
291            var lagVariableTreeNode = currentInstr.dynamicNode as LaggedVariableTreeNode;
292            int actualRow = row + lagVariableTreeNode.Lag;
293            if (actualRow < 0 || actualRow >= dataset.Rows) throw new ArgumentException("Out of range access to dataset row: " + row);
294            return dataset[actualRow, currentInstr.iArg0] * lagVariableTreeNode.Weight;
295          }
296        case OpCodes.Constant: {
297            var constTreeNode = currentInstr.dynamicNode as ConstantTreeNode;
298            return constTreeNode.Value;
299          }
300        case OpCodes.TimeLag: {
301            var timeLagTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
302            if (row + timeLagTreeNode.Lag < 0 || row + timeLagTreeNode.Lag >= dataset.Rows)
303              return double.NaN;
304
305            row += timeLagTreeNode.Lag;
306            double result = Evaluate();
307            row -= timeLagTreeNode.Lag;
308            return result;
309          }
310        case OpCodes.Integral: {
311            int nextPc = pc;
312            var timeLagTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
313            if (row + timeLagTreeNode.Lag < 0 || row + timeLagTreeNode.Lag >= dataset.Rows)
314              return double.NaN;
315            double sum = 0.0;
316            if (timeLagTreeNode.IterateNodesPrefix().OfType<VariableTreeNode>().Any()) {
317              for (int i = 0; i < Math.Abs(timeLagTreeNode.Lag); i++) {
318                row += Math.Sign(timeLagTreeNode.Lag);
319                sum += Evaluate();
320                pc = nextPc;
321              }
322              row -= timeLagTreeNode.Lag;
323              sum += Evaluate();
324            } else sum = Math.Abs(timeLagTreeNode.Lag) * Evaluate();
325            return sum;
326          }
327
328        //mkommend: derivate calculation taken from:
329        //http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
330        //one sided smooth differentiatior, N = 4
331        // y' = 1/8h (f_i + 2f_i-1, -2 f_i-3 - f_i-4)
332        case OpCodes.Derivative: {
333            if (row - 4 < 0) return double.NaN;
334            int nextPc = pc;
335            double f_0 = Evaluate(); row--;
336            pc = nextPc;
337            double f_1 = Evaluate(); row -= 2;
338            pc = nextPc;
339            double f_3 = Evaluate(); row--;
340            pc = nextPc;
341            double f_4 = Evaluate();
342            row += 4;
343
344            return (f_0 + 2 * f_1 - 2 * f_3 - f_4) / 8; // h = 1
345          }
346        default: throw new NotSupportedException();
347      }
348    }
349
350    // skips a whole branch
351    private void SkipBakedCode() {
352      int i = 1;
353      while (i > 0) {
354        i += code[pc++].nArguments;
355        i--;
356      }
357    }
358  }
359}
Note: See TracBrowser for help on using the repository browser.