Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.Problems.DataAnalysis/3.3/Symbolic/SimpleArithmeticExpressionInterpreter.cs @ 5377

Last change on this file since 5377 was 5373, checked in by mkommend, 14 years ago

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