Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.TimeSeries/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interfaces/SymbolicDataAnalysisExpressionTreeILEmittingInterpreter.cs @ 8431

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

#1081: Moved interpreter specific classes to separate directory.
File size: 45.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 System.Linq;
25using System.Reflection;
26using System.Reflection.Emit;
27using HeuristicLab.Common;
28using HeuristicLab.Core;
29using HeuristicLab.Data;
30using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
31using HeuristicLab.Parameters;
32using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
33
34namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
35  [StorableClass]
36  [Item("SymbolicDataAnalysisExpressionTreeILEmittingInterpreter", "Interpreter for symbolic expression trees.")]
37  public sealed class SymbolicDataAnalysisExpressionTreeILEmittingInterpreter : ParameterizedNamedItem,
38                                                                                ISymbolicDataAnalysisExpressionTreeInterpreter {
39    private static MethodInfo listGetValue =
40      typeof(IList<double>).GetProperty("Item", new Type[] { typeof(int) }).GetGetMethod();
41
42    private static MethodInfo cos = typeof(Math).GetMethod("Cos", new Type[] { typeof(double) });
43    private static MethodInfo sin = typeof(Math).GetMethod("Sin", new Type[] { typeof(double) });
44    private static MethodInfo tan = typeof(Math).GetMethod("Tan", new Type[] { typeof(double) });
45    private static MethodInfo exp = typeof(Math).GetMethod("Exp", new Type[] { typeof(double) });
46    private static MethodInfo log = typeof(Math).GetMethod("Log", new Type[] { typeof(double) });
47    private static MethodInfo power = typeof(Math).GetMethod("Pow", new Type[] { typeof(double), typeof(double) });
48    private static MethodInfo round = typeof(Math).GetMethod("Round", new Type[] { typeof(double) });
49    private static MethodInfo sqrt = typeof(Math).GetMethod("Sqrt", new Type[] { typeof(double) });
50
51    internal delegate double CompiledFunction(int sampleIndex, IList<double>[] columns, IList<double>[] cachedValues, int cacheStartIndex);
52    private static Type thisType = typeof(SymbolicDataAnalysisExpressionTreeILEmittingInterpreter);
53    private static MethodInfo airyA = thisType.GetMethod("AiryA", new Type[] { typeof(double) });
54    private static MethodInfo airyB = thisType.GetMethod("AiryB", new Type[] { typeof(double) });
55    private static MethodInfo gamma = thisType.GetMethod("Gamma", new Type[] { typeof(double) });
56    private static MethodInfo psi = thisType.GetMethod("Psi", new Type[] { typeof(double) });
57    private static MethodInfo dawson = thisType.GetMethod("Dawson", new Type[] { typeof(double) });
58    private static MethodInfo expIntegralEi = thisType.GetMethod("ExpIntegralEi", new Type[] { typeof(double) });
59    private static MethodInfo sinIntegral = thisType.GetMethod("SinIntegral", new Type[] { typeof(double) });
60    private static MethodInfo cosIntegral = thisType.GetMethod("CosIntegral", new Type[] { typeof(double) });
61    private static MethodInfo hypSinIntegral = thisType.GetMethod("HypSinIntegral", new Type[] { typeof(double) });
62    private static MethodInfo hypCosIntegral = thisType.GetMethod("HypCosIntegral", new Type[] { typeof(double) });
63    private static MethodInfo fresnelCosIntegral = thisType.GetMethod("FresnelCosIntegral", new Type[] { typeof(double) });
64    private static MethodInfo fresnelSinIntegral = thisType.GetMethod("FresnelSinIntegral", new Type[] { typeof(double) });
65    private static MethodInfo norm = thisType.GetMethod("Norm", new Type[] { typeof(double) });
66    private static MethodInfo erf = thisType.GetMethod("Erf", new Type[] { typeof(double) });
67    private static MethodInfo bessel = thisType.GetMethod("Bessel", new Type[] { typeof(double) });
68
69   
70
71    private const string CheckExpressionsWithIntervalArithmeticParameterName = "CheckExpressionsWithIntervalArithmetic";
72    private const string EvaluatedSolutionsParameterName = "EvaluatedSolutions";
73
74    #region private classes
75
76    private class InterpreterState {
77      private Instruction[] code;
78      private int pc;
79
80      public int ProgramCounter {
81        get { return pc; }
82        set { pc = value; }
83      }
84
85      private bool inLaggedContext;
86
87      public bool InLaggedContext {
88        get { return inLaggedContext; }
89        set { inLaggedContext = value; }
90      }
91
92      internal InterpreterState(Instruction[] code) {
93        this.inLaggedContext = false;
94        this.code = code;
95        this.pc = 0;
96      }
97
98      internal Instruction NextInstruction() {
99        return code[pc++];
100      }
101    }
102
103    private class OpCodes {
104      public const byte Add = 1;
105      public const byte Sub = 2;
106      public const byte Mul = 3;
107      public const byte Div = 4;
108
109      public const byte Sin = 5;
110      public const byte Cos = 6;
111      public const byte Tan = 7;
112
113      public const byte Log = 8;
114      public const byte Exp = 9;
115
116      public const byte IfThenElse = 10;
117
118      public const byte GT = 11;
119      public const byte LT = 12;
120
121      public const byte AND = 13;
122      public const byte OR = 14;
123      public const byte NOT = 15;
124
125
126      public const byte Average = 16;
127
128      public const byte Call = 17;
129
130      public const byte Variable = 18;
131      public const byte LagVariable = 19;
132      public const byte Constant = 20;
133      public const byte Arg = 21;
134
135      public const byte Power = 22;
136      public const byte Root = 23;
137      public const byte TimeLag = 24;
138      public const byte Integral = 25;
139      public const byte Derivative = 26;
140
141      public const byte VariableCondition = 27;
142
143      public const byte Square = 28;
144      public const byte SquareRoot = 29;
145      public const byte Gamma = 30;
146      public const byte Psi = 31;
147      public const byte Dawson = 32;
148      public const byte ExponentialIntegralEi = 33;
149      public const byte CosineIntegral = 34;
150      public const byte SineIntegral = 35;
151      public const byte HyperbolicCosineIntegral = 36;
152      public const byte HyperbolicSineIntegral = 37;
153      public const byte FresnelCosineIntegral = 38;
154      public const byte FresnelSineIntegral = 39;
155      public const byte AiryA = 40;
156      public const byte AiryB = 41;
157      public const byte Norm = 42;
158      public const byte Erf = 43;
159      public const byte Bessel = 44;
160    }
161
162    #endregion
163
164    private Dictionary<Type, byte> symbolToOpcode = new Dictionary<Type, byte>()
165                                                      {
166                                                        {typeof (Addition), OpCodes.Add},
167                                                        {typeof (Subtraction), OpCodes.Sub},
168                                                        {typeof (Multiplication), OpCodes.Mul},
169                                                        {typeof (Division), OpCodes.Div},
170                                                        {typeof (Sine), OpCodes.Sin},
171                                                        {typeof (Cosine), OpCodes.Cos},
172                                                        {typeof (Tangent), OpCodes.Tan},
173                                                        {typeof (Logarithm), OpCodes.Log},
174                                                        {typeof (Exponential), OpCodes.Exp},
175                                                        {typeof (IfThenElse), OpCodes.IfThenElse},
176                                                        {typeof (GreaterThan), OpCodes.GT},
177                                                        {typeof (LessThan), OpCodes.LT},
178                                                        {typeof (And), OpCodes.AND},
179                                                        {typeof (Or), OpCodes.OR},
180                                                        {typeof (Not), OpCodes.NOT},
181                                                        {typeof (Average), OpCodes.Average},
182                                                        {typeof (InvokeFunction), OpCodes.Call},
183                                                        {
184                                                          typeof (HeuristicLab.Problems.DataAnalysis.Symbolic.Variable),
185                                                          OpCodes.Variable
186                                                          },
187                                                        {typeof (LaggedVariable), OpCodes.LagVariable},
188                                                        {typeof (Constant), OpCodes.Constant},
189                                                        {typeof (Argument), OpCodes.Arg},
190                                                        {typeof (Power), OpCodes.Power},
191                                                        {typeof (Root), OpCodes.Root},
192                                                        {typeof (TimeLag), OpCodes.TimeLag},
193                                                        {typeof (Integral), OpCodes.Integral},
194                                                        {typeof (Derivative), OpCodes.Derivative},
195                                                        {typeof (VariableCondition), OpCodes.VariableCondition},
196                                                        {typeof (Square), OpCodes.Square},
197                                                        {typeof (SquareRoot), OpCodes.SquareRoot},
198                                                        {typeof (Gamma), OpCodes.Gamma},
199                                                        {typeof (Psi), OpCodes.Psi},
200                                                        {typeof (Dawson), OpCodes.Dawson},
201                                                        {typeof (ExponentialIntegralEi), OpCodes.ExponentialIntegralEi},
202                                                        {typeof (CosineIntegral), OpCodes.CosineIntegral},
203                                                        {typeof (SineIntegral), OpCodes.SineIntegral},
204                                                        {
205                                                          typeof (HyperbolicCosineIntegral),
206                                                          OpCodes.HyperbolicCosineIntegral
207                                                          },
208                                                        {
209                                                          typeof (HyperbolicSineIntegral), OpCodes.HyperbolicSineIntegral
210                                                          },
211                                                        {typeof (FresnelCosineIntegral), OpCodes.FresnelCosineIntegral},
212                                                        {typeof (FresnelSineIntegral), OpCodes.FresnelSineIntegral},
213                                                        {typeof (AiryA), OpCodes.AiryA},
214                                                        {typeof (AiryB), OpCodes.AiryB},
215                                                        {typeof (Norm), OpCodes.Norm},
216                                                        {typeof (Erf), OpCodes.Erf},
217                                                        {typeof (Bessel), OpCodes.Bessel}
218                                                      };
219
220    public override bool CanChangeName {
221      get { return false; }
222    }
223
224    public override bool CanChangeDescription {
225      get { return false; }
226    }
227
228    #region parameter properties
229
230    public IValueParameter<BoolValue> CheckExpressionsWithIntervalArithmeticParameter {
231      get { return (IValueParameter<BoolValue>)Parameters[CheckExpressionsWithIntervalArithmeticParameterName]; }
232    }
233
234    public IValueParameter<IntValue> EvaluatedSolutionsParameter {
235      get { return (IValueParameter<IntValue>)Parameters[EvaluatedSolutionsParameterName]; }
236    }
237
238    #endregion
239
240    #region properties
241
242    public BoolValue CheckExpressionsWithIntervalArithmetic {
243      get { return CheckExpressionsWithIntervalArithmeticParameter.Value; }
244      set { CheckExpressionsWithIntervalArithmeticParameter.Value = value; }
245    }
246
247    public IntValue EvaluatedSolutions {
248      get { return EvaluatedSolutionsParameter.Value; }
249      set { EvaluatedSolutionsParameter.Value = value; }
250    }
251
252    #endregion
253
254
255    [StorableConstructor]
256    private SymbolicDataAnalysisExpressionTreeILEmittingInterpreter(bool deserializing)
257      : base(deserializing) {
258    }
259
260    private SymbolicDataAnalysisExpressionTreeILEmittingInterpreter(
261      SymbolicDataAnalysisExpressionTreeILEmittingInterpreter original, Cloner cloner)
262      : base(original, cloner) {
263    }
264
265    public override IDeepCloneable Clone(Cloner cloner) {
266      return new SymbolicDataAnalysisExpressionTreeILEmittingInterpreter(this, cloner);
267    }
268
269    public SymbolicDataAnalysisExpressionTreeILEmittingInterpreter()
270      : base("SymbolicDataAnalysisExpressionTreeILEmittingInterpreter", "Interpreter for symbolic expression trees.") {
271      Parameters.Add(new ValueParameter<BoolValue>(CheckExpressionsWithIntervalArithmeticParameterName,
272                                                   "Switch that determines if the interpreter checks the validity of expressions with interval arithmetic before evaluating the expression.",
273                                                   new BoolValue(false)));
274      Parameters.Add(new ValueParameter<IntValue>(EvaluatedSolutionsParameterName,
275                                                  "A counter for the total number of solutions the interpreter has evaluated",
276                                                  new IntValue(0)));
277    }
278
279    [StorableHook(HookType.AfterDeserialization)]
280    private void AfterDeserialization() {
281      if (!Parameters.ContainsKey(EvaluatedSolutionsParameterName))
282        Parameters.Add(new ValueParameter<IntValue>(EvaluatedSolutionsParameterName,
283                                                    "A counter for the total number of solutions the interpreter has evaluated",
284                                                    new IntValue(0)));
285    }
286
287    #region IStatefulItem
288
289    public void InitializeState() {
290      EvaluatedSolutions.Value = 0;
291    }
292
293    public void ClearState() {
294      EvaluatedSolutions.Value = 0;
295    }
296
297    #endregion
298
299    public IEnumerable<double> GetSymbolicExpressionTreeValues(ISymbolicExpressionTree tree, Dataset dataset,
300                                                               IEnumerable<int> rows) {
301      return GetSymbolicExpressionTreeValues(tree, dataset, new string[] { "#NOTHING#" }, rows);
302    }
303    public IEnumerable<double> GetSymbolicExpressionTreeValues(ISymbolicExpressionTree tree, Dataset dataset, string[] targetVariables, IEnumerable<int> rows) {
304      return GetSymbolicExpressionTreeValues(tree, dataset, targetVariables, rows, 1);
305    }
306    // for each row for each horizon for each target variable one value
307    public IEnumerable<double> GetSymbolicExpressionTreeValues(ISymbolicExpressionTree tree, Dataset dataset, string[] targetVariables, IEnumerable<int> rows, int horizon) {
308      if (CheckExpressionsWithIntervalArithmetic.Value)
309        throw new NotSupportedException(
310          "Interval arithmetic is not yet supported in the symbolic data analysis interpreter.");
311      EvaluatedSolutions.Value++; // increment the evaluated solutions counter
312      var compiler = new SymbolicExpressionTreeCompiler();
313      Instruction[] code = compiler.Compile(tree, MapSymbolToOpCode);
314      int necessaryArgStackSize = 0;
315
316      Dictionary<string, int> doubleVariableNames =
317        dataset.DoubleVariables.Select((x, i) => new { x, i }).ToDictionary(e => e.x, e => e.i);
318
319      for (int i = 0; i < code.Length; i++) {
320        Instruction instr = code[i];
321        if (instr.opCode == OpCodes.Variable) {
322          var variableTreeNode = instr.dynamicNode as VariableTreeNode;
323          instr.iArg0 = doubleVariableNames[variableTreeNode.VariableName];
324          code[i] = instr;
325        } else if (instr.opCode == OpCodes.LagVariable) {
326          var variableTreeNode = instr.dynamicNode as LaggedVariableTreeNode;
327          instr.iArg0 = doubleVariableNames[variableTreeNode.VariableName];
328          code[i] = instr;
329        } else if (instr.opCode == OpCodes.VariableCondition) {
330          var variableConditionTreeNode = instr.dynamicNode as VariableConditionTreeNode;
331          instr.iArg0 = doubleVariableNames[variableConditionTreeNode.VariableName];
332        } else if (instr.opCode == OpCodes.Call) {
333          necessaryArgStackSize += instr.nArguments + 1;
334        }
335      }
336      var state = new InterpreterState(code);
337      Type[] methodArgs = { typeof(int), typeof(IList<double>[]), typeof(IList<double>[]), typeof(int) };
338
339      CompiledFunction[] function = new CompiledFunction[targetVariables.Length];
340      for (int i = 0; i < function.Length; i++) {
341        DynamicMethod testFun = new DynamicMethod("TestFun", typeof(double), methodArgs, typeof(SymbolicDataAnalysisExpressionTreeILEmittingInterpreter).Module);
342        ILGenerator il = testFun.GetILGenerator();
343        CompileInstructions(il, state, dataset);
344        il.Emit(System.Reflection.Emit.OpCodes.Conv_R8);
345        il.Emit(System.Reflection.Emit.OpCodes.Ret);
346        function[i] = (CompiledFunction)testFun.CreateDelegate(typeof(CompiledFunction));
347      }
348      var values = doubleVariableNames.Keys
349        .Select(v => dataset.GetReadOnlyDoubleValues(v))
350        .ToArray();
351      var cachedValues = (from var in doubleVariableNames.Keys
352                          select new double[horizon]).ToArray();
353      foreach (var row in rows) {
354        // init first line of cache
355        int c = 0;
356        foreach (var var in doubleVariableNames.Keys)
357          cachedValues[c++][0] = dataset.GetDoubleValue(var, row);
358        for (int horizonRow = row; horizonRow < row + horizon; horizonRow++) {
359          for (int i = 0; i < function.Length; i++) {
360            var componentProg = function[i](horizonRow, values, cachedValues, row);
361            // set cachedValues for prognosis of future values
362            if (horizon > 1)
363              cachedValues[doubleVariableNames[targetVariables[i]]][horizonRow - row] = componentProg;
364            yield return componentProg;
365          }
366        }
367      }
368    }
369
370    private void CompileInstructions(ILGenerator il, InterpreterState state, Dataset ds) {
371      Instruction currentInstr = state.NextInstruction();
372      int nArgs = currentInstr.nArguments;
373
374      switch (currentInstr.opCode) {
375        case OpCodes.Add: {
376            if (nArgs > 0) {
377              CompileInstructions(il, state, ds);
378            }
379            for (int i = 1; i < nArgs; i++) {
380              CompileInstructions(il, state, ds);
381              il.Emit(System.Reflection.Emit.OpCodes.Add);
382            }
383            return;
384          }
385        case OpCodes.Sub: {
386            if (nArgs == 1) {
387              CompileInstructions(il, state, ds);
388              il.Emit(System.Reflection.Emit.OpCodes.Neg);
389              return;
390            }
391            if (nArgs > 0) {
392              CompileInstructions(il, state, ds);
393            }
394            for (int i = 1; i < nArgs; i++) {
395              CompileInstructions(il, state, ds);
396              il.Emit(System.Reflection.Emit.OpCodes.Sub);
397            }
398            return;
399          }
400        case OpCodes.Mul: {
401            if (nArgs > 0) {
402              CompileInstructions(il, state, ds);
403            }
404            for (int i = 1; i < nArgs; i++) {
405              CompileInstructions(il, state, ds);
406              il.Emit(System.Reflection.Emit.OpCodes.Mul);
407            }
408            return;
409          }
410        case OpCodes.Div: {
411            if (nArgs == 1) {
412              il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0);
413              CompileInstructions(il, state, ds);
414              il.Emit(System.Reflection.Emit.OpCodes.Div);
415              return;
416            }
417            if (nArgs > 0) {
418              CompileInstructions(il, state, ds);
419            }
420            for (int i = 1; i < nArgs; i++) {
421              CompileInstructions(il, state, ds);
422              il.Emit(System.Reflection.Emit.OpCodes.Div);
423            }
424            return;
425          }
426        case OpCodes.Average: {
427            CompileInstructions(il, state, ds);
428            for (int i = 1; i < nArgs; i++) {
429              CompileInstructions(il, state, ds);
430              il.Emit(System.Reflection.Emit.OpCodes.Add);
431            }
432            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, nArgs);
433            il.Emit(System.Reflection.Emit.OpCodes.Div);
434            return;
435          }
436        case OpCodes.Cos: {
437            CompileInstructions(il, state, ds);
438            il.Emit(System.Reflection.Emit.OpCodes.Call, cos);
439            return;
440          }
441        case OpCodes.Sin: {
442            CompileInstructions(il, state, ds);
443            il.Emit(System.Reflection.Emit.OpCodes.Call, sin);
444            return;
445          }
446        case OpCodes.Tan: {
447            CompileInstructions(il, state, ds);
448            il.Emit(System.Reflection.Emit.OpCodes.Call, tan);
449            return;
450          }
451        case OpCodes.Power: {
452            CompileInstructions(il, state, ds);
453            CompileInstructions(il, state, ds);
454            il.Emit(System.Reflection.Emit.OpCodes.Call, round);
455            il.Emit(System.Reflection.Emit.OpCodes.Call, power);
456            return;
457          }
458        case OpCodes.Root: {
459            CompileInstructions(il, state, ds);
460            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // 1 / round(...)
461            CompileInstructions(il, state, ds);
462            il.Emit(System.Reflection.Emit.OpCodes.Call, round);
463            il.Emit(System.Reflection.Emit.OpCodes.Div);
464            il.Emit(System.Reflection.Emit.OpCodes.Call, power);
465            return;
466          }
467        case OpCodes.Exp: {
468            CompileInstructions(il, state, ds);
469            il.Emit(System.Reflection.Emit.OpCodes.Call, exp);
470            return;
471          }
472        case OpCodes.Log: {
473            CompileInstructions(il, state, ds);
474            il.Emit(System.Reflection.Emit.OpCodes.Call, log);
475            return;
476          }
477        case OpCodes.Square: {
478            CompileInstructions(il, state, ds);
479            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 2.0);
480            il.Emit(System.Reflection.Emit.OpCodes.Call, power);
481            return;
482          }
483        case OpCodes.SquareRoot: {
484            CompileInstructions(il, state, ds);
485            il.Emit(System.Reflection.Emit.OpCodes.Call, sqrt);
486            return;
487          }
488        case OpCodes.AiryA: {
489            CompileInstructions(il, state, ds);
490            il.Emit(System.Reflection.Emit.OpCodes.Call, airyA);
491            return;
492          }
493        case OpCodes.AiryB: {
494            CompileInstructions(il, state, ds);
495            il.Emit(System.Reflection.Emit.OpCodes.Call, airyB);
496            return;
497          }
498        case OpCodes.Bessel: {
499            CompileInstructions(il, state, ds);
500            il.Emit(System.Reflection.Emit.OpCodes.Call, bessel);
501            return;
502          }
503        case OpCodes.CosineIntegral: {
504            CompileInstructions(il, state, ds);
505            il.Emit(System.Reflection.Emit.OpCodes.Call, cosIntegral);
506            return;
507          }
508        case OpCodes.Dawson: {
509            CompileInstructions(il, state, ds);
510            il.Emit(System.Reflection.Emit.OpCodes.Call, dawson);
511            return;
512          }
513        case OpCodes.Erf: {
514            CompileInstructions(il, state, ds);
515            il.Emit(System.Reflection.Emit.OpCodes.Call, erf);
516            return;
517          }
518        case OpCodes.ExponentialIntegralEi: {
519            CompileInstructions(il, state, ds);
520            il.Emit(System.Reflection.Emit.OpCodes.Call, expIntegralEi);
521            return;
522          }
523        case OpCodes.FresnelCosineIntegral: {
524            CompileInstructions(il, state, ds);
525            il.Emit(System.Reflection.Emit.OpCodes.Call, fresnelCosIntegral);
526            return;
527          }
528        case OpCodes.FresnelSineIntegral: {
529            CompileInstructions(il, state, ds);
530            il.Emit(System.Reflection.Emit.OpCodes.Call, fresnelSinIntegral);
531            return;
532          }
533        case OpCodes.Gamma: {
534            CompileInstructions(il, state, ds);
535            il.Emit(System.Reflection.Emit.OpCodes.Call, gamma);
536            return;
537          }
538        case OpCodes.HyperbolicCosineIntegral: {
539            CompileInstructions(il, state, ds);
540            il.Emit(System.Reflection.Emit.OpCodes.Call, hypCosIntegral);
541            return;
542          }
543        case OpCodes.HyperbolicSineIntegral: {
544            CompileInstructions(il, state, ds);
545            il.Emit(System.Reflection.Emit.OpCodes.Call, hypSinIntegral);
546            return;
547          }
548        case OpCodes.Norm: {
549            CompileInstructions(il, state, ds);
550            il.Emit(System.Reflection.Emit.OpCodes.Call, norm);
551            return;
552          }
553        case OpCodes.Psi: {
554            CompileInstructions(il, state, ds);
555            il.Emit(System.Reflection.Emit.OpCodes.Call, psi);
556            return;
557          }
558        case OpCodes.SineIntegral: {
559            CompileInstructions(il, state, ds);
560            il.Emit(System.Reflection.Emit.OpCodes.Call, sinIntegral);
561            return;
562          }
563        case OpCodes.IfThenElse: {
564            Label end = il.DefineLabel();
565            Label c1 = il.DefineLabel();
566            CompileInstructions(il, state, ds);
567            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_0); // > 0
568            il.Emit(System.Reflection.Emit.OpCodes.Cgt);
569            il.Emit(System.Reflection.Emit.OpCodes.Brfalse, c1);
570            CompileInstructions(il, state, ds);
571            il.Emit(System.Reflection.Emit.OpCodes.Br, end);
572            il.MarkLabel(c1);
573            CompileInstructions(il, state, ds);
574            il.MarkLabel(end);
575            return;
576          }
577        case OpCodes.AND: {
578            Label falseBranch = il.DefineLabel();
579            Label end = il.DefineLabel();
580            CompileInstructions(il, state, ds);
581            for (int i = 1; i < nArgs; i++) {
582              il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_0); // > 0
583              il.Emit(System.Reflection.Emit.OpCodes.Cgt);
584              il.Emit(System.Reflection.Emit.OpCodes.Brfalse, falseBranch);
585              CompileInstructions(il, state, ds);
586            }
587            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_0); // > 0
588            il.Emit(System.Reflection.Emit.OpCodes.Cgt);
589            il.Emit(System.Reflection.Emit.OpCodes.Brfalse, falseBranch);
590            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // 1
591            il.Emit(System.Reflection.Emit.OpCodes.Br, end);
592            il.MarkLabel(falseBranch);
593            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // -1
594            il.Emit(System.Reflection.Emit.OpCodes.Neg);
595            il.MarkLabel(end);
596            return;
597          }
598        case OpCodes.OR: {
599            Label trueBranch = il.DefineLabel();
600            Label end = il.DefineLabel();
601            Label resultBranch = il.DefineLabel();
602            CompileInstructions(il, state, ds);
603            for (int i = 1; i < nArgs; i++) {
604              Label nextArgBranch = il.DefineLabel();
605              // complex definition because of special properties of NaN 
606              il.Emit(System.Reflection.Emit.OpCodes.Dup);
607              il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_0); // <= 0       
608              il.Emit(System.Reflection.Emit.OpCodes.Ble, nextArgBranch);
609              il.Emit(System.Reflection.Emit.OpCodes.Br, resultBranch);
610              il.MarkLabel(nextArgBranch);
611              il.Emit(System.Reflection.Emit.OpCodes.Pop);
612              CompileInstructions(il, state, ds);
613            }
614            il.MarkLabel(resultBranch);
615            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_0); // > 0
616            il.Emit(System.Reflection.Emit.OpCodes.Cgt);
617            il.Emit(System.Reflection.Emit.OpCodes.Brtrue, trueBranch);
618            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // -1
619            il.Emit(System.Reflection.Emit.OpCodes.Neg);
620            il.Emit(System.Reflection.Emit.OpCodes.Br, end);
621            il.MarkLabel(trueBranch);
622            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // 1
623            il.MarkLabel(end);
624            return;
625          }
626        case OpCodes.NOT: {
627            CompileInstructions(il, state, ds);
628            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_0); // > 0
629            il.Emit(System.Reflection.Emit.OpCodes.Cgt);
630            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 2.0); // * 2
631            il.Emit(System.Reflection.Emit.OpCodes.Mul);
632            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // - 1
633            il.Emit(System.Reflection.Emit.OpCodes.Sub);
634            il.Emit(System.Reflection.Emit.OpCodes.Neg); // * -1
635            return;
636          }
637        case OpCodes.GT: {
638            CompileInstructions(il, state, ds);
639            CompileInstructions(il, state, ds);
640
641            il.Emit(System.Reflection.Emit.OpCodes.Cgt); // 1 (>) / 0 (otherwise)
642            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 2.0); // * 2
643            il.Emit(System.Reflection.Emit.OpCodes.Mul);
644            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // - 1
645            il.Emit(System.Reflection.Emit.OpCodes.Sub);
646            return;
647          }
648        case OpCodes.LT: {
649            CompileInstructions(il, state, ds);
650            CompileInstructions(il, state, ds);
651            il.Emit(System.Reflection.Emit.OpCodes.Clt);
652            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 2.0); // * 2
653            il.Emit(System.Reflection.Emit.OpCodes.Mul);
654            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // - 1
655            il.Emit(System.Reflection.Emit.OpCodes.Sub);
656            return;
657          }
658        case OpCodes.TimeLag: {
659            LaggedTreeNode laggedTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
660            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row -= lag
661            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, laggedTreeNode.Lag);
662            il.Emit(System.Reflection.Emit.OpCodes.Add);
663            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
664            var prevLaggedContext = state.InLaggedContext;
665            state.InLaggedContext = true;
666            CompileInstructions(il, state, ds);
667            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row += lag
668            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, laggedTreeNode.Lag);
669            il.Emit(System.Reflection.Emit.OpCodes.Sub);
670            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
671            state.InLaggedContext = prevLaggedContext;
672            return;
673          }
674        case OpCodes.Integral: {
675            int savedPc = state.ProgramCounter;
676            LaggedTreeNode laggedTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
677            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row -= lag
678            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, laggedTreeNode.Lag);
679            il.Emit(System.Reflection.Emit.OpCodes.Add);
680            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
681            var prevLaggedContext = state.InLaggedContext;
682            state.InLaggedContext = true;
683            CompileInstructions(il, state, ds);
684            for (int l = laggedTreeNode.Lag; l < 0; l++) {
685              il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row += lag
686              il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_1);
687              il.Emit(System.Reflection.Emit.OpCodes.Add);
688              il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
689              state.ProgramCounter = savedPc;
690              CompileInstructions(il, state, ds);
691              il.Emit(System.Reflection.Emit.OpCodes.Add);
692            }
693            state.InLaggedContext = prevLaggedContext;
694            return;
695          }
696
697        //mkommend: derivate calculation taken from:
698        //http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
699        //one sided smooth differentiatior, N = 4
700        // y' = 1/8h (f_i + 2f_i-1, -2 f_i-3 - f_i-4)
701        case OpCodes.Derivative: {
702            int savedPc = state.ProgramCounter;
703            CompileInstructions(il, state, ds);
704            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row --
705            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_M1);
706            il.Emit(System.Reflection.Emit.OpCodes.Add);
707            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
708            state.ProgramCounter = savedPc;
709            var prevLaggedContext = state.InLaggedContext;
710            state.InLaggedContext = true;
711            CompileInstructions(il, state, ds);
712            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 2.0); // f_0 + 2 * f_1
713            il.Emit(System.Reflection.Emit.OpCodes.Mul);
714            il.Emit(System.Reflection.Emit.OpCodes.Add);
715
716            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row -=2
717            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_2);
718            il.Emit(System.Reflection.Emit.OpCodes.Sub);
719            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
720            state.ProgramCounter = savedPc;
721            CompileInstructions(il, state, ds);
722            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 2.0); // f_0 + 2 * f_1 - 2 * f_3
723            il.Emit(System.Reflection.Emit.OpCodes.Mul);
724            il.Emit(System.Reflection.Emit.OpCodes.Sub);
725
726            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row --
727            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_M1);
728            il.Emit(System.Reflection.Emit.OpCodes.Add);
729            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
730            state.ProgramCounter = savedPc;
731            CompileInstructions(il, state, ds);
732            il.Emit(System.Reflection.Emit.OpCodes.Sub); // f_0 + 2 * f_1 - 2 * f_3 - f_4
733            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 8.0); // / 8
734            il.Emit(System.Reflection.Emit.OpCodes.Div);
735
736            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row +=4
737            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_4);
738            il.Emit(System.Reflection.Emit.OpCodes.Add);
739            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
740            state.InLaggedContext = prevLaggedContext;
741            return;
742          }
743        case OpCodes.Call: {
744            throw new NotSupportedException(
745              "Automatically defined functions are not supported by the SymbolicDataAnalysisTreeILEmittingInterpreter. Either turn of ADFs or change the interpeter.");
746          }
747        case OpCodes.Arg: {
748            throw new NotSupportedException(
749              "Automatically defined functions are not supported by the SymbolicDataAnalysisTreeILEmittingInterpreter. Either turn of ADFs or change the interpeter.");
750          }
751        case OpCodes.Variable: {
752            VariableTreeNode varNode = (VariableTreeNode)currentInstr.dynamicNode;
753            if (!state.InLaggedContext) {
754              il.Emit(System.Reflection.Emit.OpCodes.Ldarg_1); // load columns array
755              il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, (int)currentInstr.iArg0);
756              // load correct column of the current variable
757              il.Emit(System.Reflection.Emit.OpCodes.Ldelem_Ref);
758              il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // rowIndex
759              il.Emit(System.Reflection.Emit.OpCodes.Call, listGetValue);
760              il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, varNode.Weight); // load weight
761              il.Emit(System.Reflection.Emit.OpCodes.Mul);
762            } else {
763              var nanResult = il.DefineLabel();
764              var normalResult = il.DefineLabel();
765              var cachedValue = il.DefineLabel();
766              var multiplyValue = il.DefineLabel();
767              il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // rowIndex
768              il.Emit(System.Reflection.Emit.OpCodes.Dup);
769              il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_0);
770              il.Emit(System.Reflection.Emit.OpCodes.Blt, nanResult);
771              il.Emit(System.Reflection.Emit.OpCodes.Dup);
772              il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, ds.Rows);
773              il.Emit(System.Reflection.Emit.OpCodes.Bge, nanResult);
774              il.Emit(System.Reflection.Emit.OpCodes.Ldarg_3);
775              il.Emit(System.Reflection.Emit.OpCodes.Bge, cachedValue);
776              il.Emit(System.Reflection.Emit.OpCodes.Ldarg_1); // load columns array
777              il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, (int)currentInstr.iArg0);
778              il.Emit(System.Reflection.Emit.OpCodes.Ldelem_Ref);
779              il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // rowIndex
780              il.Emit(System.Reflection.Emit.OpCodes.Call, listGetValue);
781              il.Emit(System.Reflection.Emit.OpCodes.Br, multiplyValue);
782              il.MarkLabel(cachedValue);
783              il.Emit(System.Reflection.Emit.OpCodes.Ldarg_2); // load cached values array
784              il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, (int)currentInstr.iArg0);
785              il.Emit(System.Reflection.Emit.OpCodes.Ldelem_Ref);
786              il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // rowIndex
787              il.Emit(System.Reflection.Emit.OpCodes.Ldarg_3); // startRow
788              il.Emit(System.Reflection.Emit.OpCodes.Sub); // startRow
789              il.Emit(System.Reflection.Emit.OpCodes.Call, listGetValue);
790              il.MarkLabel(multiplyValue);
791              il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, varNode.Weight); // load weight
792              il.Emit(System.Reflection.Emit.OpCodes.Mul);
793              il.Emit(System.Reflection.Emit.OpCodes.Br, normalResult);
794              il.MarkLabel(nanResult);
795              il.Emit(System.Reflection.Emit.OpCodes.Pop); // rowIndex
796              il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, double.NaN);
797              il.MarkLabel(normalResult);
798            }
799            return;
800          }
801        case OpCodes.LagVariable: {
802            var nanResult = il.DefineLabel();
803            var normalResult = il.DefineLabel();
804            var cachedValue = il.DefineLabel();
805            var multiplyValue = il.DefineLabel();
806            LaggedVariableTreeNode varNode = (LaggedVariableTreeNode)currentInstr.dynamicNode;
807            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, varNode.Lag); // lag
808            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // rowIndex
809            il.Emit(System.Reflection.Emit.OpCodes.Add); // actualRowIndex = rowIndex + sampleOffset
810            il.Emit(System.Reflection.Emit.OpCodes.Dup);
811            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_0);
812            il.Emit(System.Reflection.Emit.OpCodes.Blt, nanResult);
813            il.Emit(System.Reflection.Emit.OpCodes.Dup);
814            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, ds.Rows);
815            il.Emit(System.Reflection.Emit.OpCodes.Bge, nanResult);
816            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_3); // startindex
817            il.Emit(System.Reflection.Emit.OpCodes.Bge, cachedValue);
818            // normal value
819            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_1); // load columns array
820            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, (int)currentInstr.iArg0); // load correct column of the current variable
821            il.Emit(System.Reflection.Emit.OpCodes.Ldelem_Ref);
822            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, varNode.Lag); // lag
823            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // rowIndex
824            il.Emit(System.Reflection.Emit.OpCodes.Add); // actualRowIndex = rowIndex + sampleOffset
825            il.Emit(System.Reflection.Emit.OpCodes.Call, listGetValue);
826            il.Emit(System.Reflection.Emit.OpCodes.Br, multiplyValue);
827            il.MarkLabel(cachedValue);
828            // cached value
829            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_2); // load cached values
830            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, (int)currentInstr.iArg0); // load correct column of the current variable
831            il.Emit(System.Reflection.Emit.OpCodes.Ldelem_Ref);
832            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, varNode.Lag); // lag
833            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // rowIndex
834            il.Emit(System.Reflection.Emit.OpCodes.Add); // actualRowIndex = rowIndex + sampleOffset
835            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_3); // startRow
836            il.Emit(System.Reflection.Emit.OpCodes.Sub); // startRow           
837            il.Emit(System.Reflection.Emit.OpCodes.Call, listGetValue);
838
839            il.MarkLabel(multiplyValue);
840            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, varNode.Weight); // load weight
841            il.Emit(System.Reflection.Emit.OpCodes.Mul);
842            il.Emit(System.Reflection.Emit.OpCodes.Br, normalResult);
843            il.MarkLabel(nanResult);
844            il.Emit(System.Reflection.Emit.OpCodes.Pop); // pop the row index
845            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, double.NaN);
846            il.MarkLabel(normalResult);
847            return;
848          }
849        case OpCodes.Constant: {
850            ConstantTreeNode constNode = (ConstantTreeNode)currentInstr.dynamicNode;
851            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, constNode.Value);
852            return;
853          }
854
855        //mkommend: this symbol uses the logistic function f(x) = 1 / (1 + e^(-alpha * x) )
856        //to determine the relative amounts of the true and false branch see http://en.wikipedia.org/wiki/Logistic_function
857        case OpCodes.VariableCondition: {
858            throw new NotSupportedException("Interpretation of symbol " + currentInstr.dynamicNode.Symbol.Name +
859                                            " is not supported by the SymbolicDataAnalysisTreeILEmittingInterpreter");
860          }
861        default:
862          throw new NotSupportedException("Interpretation of symbol " + currentInstr.dynamicNode.Symbol.Name +
863                                          " is not supported by the SymbolicDataAnalysisTreeILEmittingInterpreter");
864      }
865    }
866
867    private byte MapSymbolToOpCode(ISymbolicExpressionTreeNode treeNode) {
868      byte opCode;
869      if (!symbolToOpcode.TryGetValue(treeNode.Symbol.GetType(), out opCode))
870        throw new NotSupportedException("Symbol: " + treeNode.Symbol);
871      return opCode;
872    }
873
874
875    public static double AiryA(double x) {
876      if (double.IsNaN(x)) return double.NaN;
877      double ai, aip, bi, bip;
878      alglib.airy(x, out ai, out aip, out bi, out bip);
879      return ai;
880    }
881
882    public static double AiryB(double x) {
883      if (double.IsNaN(x)) return double.NaN;
884      double ai, aip, bi, bip;
885      alglib.airy(x, out ai, out aip, out bi, out bip);
886      return bi;
887    }
888    public static double Dawson(double x) {
889      if (double.IsNaN(x)) return double.NaN;
890      return alglib.dawsonintegral(x);
891    }
892
893    public static double Gamma(double x) {
894      if (double.IsNaN(x)) return double.NaN;
895      return alglib.gammafunction(x);
896    }
897
898    public static double Psi(double x) {
899      if (double.IsNaN(x)) return double.NaN;
900      else if (x <= 0 && (Math.Floor(x) - x).IsAlmost(0)) return double.NaN;
901      return alglib.psi(x);
902    }
903
904    public static double ExpIntegralEi(double x) {
905      if (double.IsNaN(x)) return double.NaN;
906      return alglib.exponentialintegralei(x);
907    }
908
909    public static double SinIntegral(double x) {
910      if (double.IsNaN(x)) return double.NaN;
911      double si, ci;
912      alglib.sinecosineintegrals(x, out si, out ci);
913      return si;
914    }
915
916    public static double CosIntegral(double x) {
917      if (double.IsNaN(x)) return double.NaN;
918      double si, ci;
919      alglib.sinecosineintegrals(x, out si, out ci);
920      return ci;
921    }
922
923    public static double HypSinIntegral(double x) {
924      if (double.IsNaN(x)) return double.NaN;
925      double shi, chi;
926      alglib.hyperbolicsinecosineintegrals(x, out shi, out chi);
927      return shi;
928    }
929
930    public static double HypCosIntegral(double x) {
931      if (double.IsNaN(x)) return double.NaN;
932      double shi, chi;
933      alglib.hyperbolicsinecosineintegrals(x, out shi, out chi);
934      return chi;
935    }
936
937    public static double FresnelCosIntegral(double x) {
938      if (double.IsNaN(x)) return double.NaN;
939      double c = 0, s = 0;
940      alglib.fresnelintegral(x, ref c, ref s);
941      return c;
942    }
943
944    public static double FresnelSinIntegral(double x) {
945      if (double.IsNaN(x)) return double.NaN;
946      double c = 0, s = 0;
947      alglib.fresnelintegral(x, ref c, ref s);
948      return s;
949    }
950
951    public static double Norm(double x) {
952      if (double.IsNaN(x)) return double.NaN;
953      return alglib.normaldistribution(x);
954    }
955
956    public static double Erf(double x) {
957      if (double.IsNaN(x)) return double.NaN;
958      return alglib.errorfunction(x);
959    }
960
961    public static double Bessel(double x) {
962      if (double.IsNaN(x)) return double.NaN;
963      return alglib.besseli0(x);
964    }
965
966  }
967}
Note: See TracBrowser for help on using the repository browser.