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

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

Added new symbols for GP (ticket #1256).

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