Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/SymbolicDataAnalysisExpressionTreeILEmittingInterpreter.cs @ 16762

Last change on this file since 16762 was 16670, checked in by gkronber, 5 years ago

#2866: bugfixes for r16668
File size: 38.4 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2019 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 HEAL.Attic;
33
34namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
35  [StorableType("426718E3-2A57-4CA4-98A1-65EDD0B0BDBF")]
36  [Item("SymbolicDataAnalysisExpressionTreeILEmittingInterpreter", "Interpreter for symbolic expression trees.")]
37  public sealed class SymbolicDataAnalysisExpressionTreeILEmittingInterpreter : ParameterizedNamedItem, ISymbolicDataAnalysisExpressionTreeInterpreter {
38    private static readonly Type thisType = typeof(SymbolicDataAnalysisExpressionTreeILEmittingInterpreter);
39    internal delegate double CompiledFunction(int sampleIndex, IList<double>[] columns);
40
41    #region method infos
42    private static MethodInfo listGetValue = typeof(IList<double>).GetProperty("Item", new Type[] { typeof(int) }).GetGetMethod();
43
44    private static MethodInfo cos = typeof(Math).GetMethod("Cos", new Type[] { typeof(double) });
45    private static MethodInfo sin = typeof(Math).GetMethod("Sin", new Type[] { typeof(double) });
46    private static MethodInfo tan = typeof(Math).GetMethod("Tan", new Type[] { typeof(double) });
47    private static MethodInfo tanh = typeof(Math).GetMethod("Tanh", new Type[] { typeof(double) });
48    private static MethodInfo exp = typeof(Math).GetMethod("Exp", new Type[] { typeof(double) });
49    private static MethodInfo log = typeof(Math).GetMethod("Log", new Type[] { typeof(double) });
50    private static MethodInfo power = typeof(Math).GetMethod("Pow", new Type[] { typeof(double), typeof(double) });
51    private static MethodInfo round = typeof(Math).GetMethod("Round", new Type[] { typeof(double) });
52    private static MethodInfo sqrt = typeof(Math).GetMethod("Sqrt", new Type[] { typeof(double) });
53    private static MethodInfo abs = typeof(Math).GetMethod("Abs", new Type[] { typeof(double) });
54
55    private static MethodInfo airyA = thisType.GetMethod("AiryA", new Type[] { typeof(double) });
56    private static MethodInfo airyB = thisType.GetMethod("AiryB", new Type[] { typeof(double) });
57    private static MethodInfo gamma = thisType.GetMethod("Gamma", new Type[] { typeof(double) });
58    private static MethodInfo psi = thisType.GetMethod("Psi", new Type[] { typeof(double) });
59    private static MethodInfo dawson = thisType.GetMethod("Dawson", new Type[] { typeof(double) });
60    private static MethodInfo expIntegralEi = thisType.GetMethod("ExpIntegralEi", new Type[] { typeof(double) });
61    private static MethodInfo sinIntegral = thisType.GetMethod("SinIntegral", new Type[] { typeof(double) });
62    private static MethodInfo cosIntegral = thisType.GetMethod("CosIntegral", new Type[] { typeof(double) });
63    private static MethodInfo hypSinIntegral = thisType.GetMethod("HypSinIntegral", new Type[] { typeof(double) });
64    private static MethodInfo hypCosIntegral = thisType.GetMethod("HypCosIntegral", new Type[] { typeof(double) });
65    private static MethodInfo fresnelCosIntegral = thisType.GetMethod("FresnelCosIntegral", new Type[] { typeof(double) });
66    private static MethodInfo fresnelSinIntegral = thisType.GetMethod("FresnelSinIntegral", new Type[] { typeof(double) });
67    private static MethodInfo norm = thisType.GetMethod("Norm", new Type[] { typeof(double) });
68    private static MethodInfo erf = thisType.GetMethod("Erf", new Type[] { typeof(double) });
69    private static MethodInfo bessel = thisType.GetMethod("Bessel", new Type[] { typeof(double) });
70    private static MethodInfo string_eq = typeof(string).GetMethod("Equals", new Type[] { typeof(string) });
71    #endregion
72
73    private const string CheckExpressionsWithIntervalArithmeticParameterName = "CheckExpressionsWithIntervalArithmetic";
74    private const string CheckExpressionsWithIntervalArithmeticParameterDescription = "Switch that determines if the interpreter checks the validity of expressions with interval arithmetic before evaluating the expression.";
75    private const string EvaluatedSolutionsParameterName = "EvaluatedSolutions";
76
77    public override bool CanChangeName {
78      get { return false; }
79    }
80
81    public override bool CanChangeDescription {
82      get { return false; }
83    }
84
85    #region parameter properties
86    public IFixedValueParameter<BoolValue> CheckExpressionsWithIntervalArithmeticParameter {
87      get { return (IFixedValueParameter<BoolValue>)Parameters[CheckExpressionsWithIntervalArithmeticParameterName]; }
88    }
89
90    public IFixedValueParameter<IntValue> EvaluatedSolutionsParameter {
91      get { return (IFixedValueParameter<IntValue>)Parameters[EvaluatedSolutionsParameterName]; }
92    }
93    #endregion
94
95    #region properties
96    public bool CheckExpressionsWithIntervalArithmetic {
97      get { return CheckExpressionsWithIntervalArithmeticParameter.Value.Value; }
98      set { CheckExpressionsWithIntervalArithmeticParameter.Value.Value = value; }
99    }
100    public int EvaluatedSolutions {
101      get { return EvaluatedSolutionsParameter.Value.Value; }
102      set { EvaluatedSolutionsParameter.Value.Value = value; }
103    }
104    #endregion
105
106    [StorableConstructor]
107    private SymbolicDataAnalysisExpressionTreeILEmittingInterpreter(StorableConstructorFlag _) : base(_) { }
108
109    private SymbolicDataAnalysisExpressionTreeILEmittingInterpreter(SymbolicDataAnalysisExpressionTreeILEmittingInterpreter original, Cloner cloner) : base(original, cloner) { }
110    public override IDeepCloneable Clone(Cloner cloner) {
111      return new SymbolicDataAnalysisExpressionTreeILEmittingInterpreter(this, cloner);
112    }
113
114    public SymbolicDataAnalysisExpressionTreeILEmittingInterpreter()
115      : base("SymbolicDataAnalysisExpressionTreeILEmittingInterpreter", "Interpreter for symbolic expression trees.") {
116      Parameters.Add(new FixedValueParameter<BoolValue>(CheckExpressionsWithIntervalArithmeticParameterName,
117        "Switch that determines if the interpreter checks the validity of expressions with interval arithmetic before evaluating the expression.", new BoolValue(false)));
118      Parameters.Add(new FixedValueParameter<IntValue>(EvaluatedSolutionsParameterName, "A counter for the total number of solutions the interpreter has evaluated", new IntValue(0)));
119    }
120
121    public SymbolicDataAnalysisExpressionTreeILEmittingInterpreter(string name, string description)
122      : base(name, description) {
123      Parameters.Add(new FixedValueParameter<BoolValue>(CheckExpressionsWithIntervalArithmeticParameterName,
124        "Switch that determines if the interpreter checks the validity of expressions with interval arithmetic before evaluating the expression.", new BoolValue(false)));
125      Parameters.Add(new FixedValueParameter<IntValue>(EvaluatedSolutionsParameterName, "A counter for the total number of solutions the interpreter has evaluated", new IntValue(0)));
126    }
127
128    [StorableHook(HookType.AfterDeserialization)]
129    private void AfterDeserialization() {
130      var evaluatedSolutions = new IntValue(0);
131      var checkExpressionsWithIntervalArithmetic = new BoolValue(false);
132      if (Parameters.ContainsKey(EvaluatedSolutionsParameterName)) {
133        var evaluatedSolutionsParameter = (IValueParameter<IntValue>)Parameters[EvaluatedSolutionsParameterName];
134        evaluatedSolutions = evaluatedSolutionsParameter.Value;
135        Parameters.Remove(EvaluatedSolutionsParameterName);
136      }
137      Parameters.Add(new FixedValueParameter<IntValue>(EvaluatedSolutionsParameterName, "A counter for the total number of solutions the interpreter has evaluated", evaluatedSolutions));
138      if (Parameters.ContainsKey(CheckExpressionsWithIntervalArithmeticParameterName)) {
139        var checkExpressionsWithIntervalArithmeticParameter = (IValueParameter<BoolValue>)Parameters[CheckExpressionsWithIntervalArithmeticParameterName];
140        Parameters.Remove(CheckExpressionsWithIntervalArithmeticParameterName);
141        checkExpressionsWithIntervalArithmetic = checkExpressionsWithIntervalArithmeticParameter.Value;
142      }
143      Parameters.Add(new FixedValueParameter<BoolValue>(CheckExpressionsWithIntervalArithmeticParameterName, CheckExpressionsWithIntervalArithmeticParameterDescription, checkExpressionsWithIntervalArithmetic));
144    }
145
146    #region IStatefulItem
147    public void InitializeState() {
148      EvaluatedSolutions = 0;
149    }
150
151    public void ClearState() {
152    }
153    #endregion
154
155    private readonly object syncRoot = new object();
156    public IEnumerable<double> GetSymbolicExpressionTreeValues(ISymbolicExpressionTree tree, IDataset dataset, IEnumerable<int> rows) {
157      if (CheckExpressionsWithIntervalArithmetic)
158        throw new NotSupportedException("Interval arithmetic is not yet supported in the symbolic data analysis interpreter.");
159
160      lock (syncRoot) {
161        EvaluatedSolutions++; // increment the evaluated solutions counter
162      }
163      var state = PrepareInterpreterState(tree, dataset);
164
165      Type[] methodArgs = { typeof(int), typeof(IList<double>[]) };
166      DynamicMethod testFun = new DynamicMethod("TestFun", typeof(double), methodArgs, typeof(SymbolicDataAnalysisExpressionTreeILEmittingInterpreter).Module);
167
168      ILGenerator il = testFun.GetILGenerator();
169      CompileInstructions(il, state, dataset);
170      il.Emit(System.Reflection.Emit.OpCodes.Conv_R8);
171      il.Emit(System.Reflection.Emit.OpCodes.Ret);
172      var function = (CompiledFunction)testFun.CreateDelegate(typeof(CompiledFunction));
173
174      IList<double>[] columns = dataset.DoubleVariables.Select(v => dataset.GetReadOnlyDoubleValues(v)).ToArray();
175
176      foreach (var row in rows) {
177        yield return function(row, columns);
178      }
179    }
180
181    private InterpreterState PrepareInterpreterState(ISymbolicExpressionTree tree, IDataset dataset) {
182      Instruction[] code = SymbolicExpressionTreeCompiler.Compile(tree, OpCodes.MapSymbolToOpCode);
183      Dictionary<string, int> doubleVariableNames = dataset.DoubleVariables.Select((x, i) => new { x, i }).ToDictionary(e => e.x, e => e.i);
184      int necessaryArgStackSize = 0;
185      foreach (Instruction instr in code) {
186        if (instr.opCode == OpCodes.Variable) {
187          var variableTreeNode = (VariableTreeNode)instr.dynamicNode;
188          instr.data = doubleVariableNames[variableTreeNode.VariableName];
189        } else if (instr.opCode == OpCodes.LagVariable) {
190          var laggedVariableTreeNode = (LaggedVariableTreeNode)instr.dynamicNode;
191          instr.data = doubleVariableNames[laggedVariableTreeNode.VariableName];
192        } else if (instr.opCode == OpCodes.VariableCondition) {
193          var variableConditionTreeNode = (VariableConditionTreeNode)instr.dynamicNode;
194          instr.data = doubleVariableNames[variableConditionTreeNode.VariableName];
195        } else if (instr.opCode == OpCodes.Call) {
196          necessaryArgStackSize += instr.nArguments + 1;
197        }
198      }
199      return new InterpreterState(code, necessaryArgStackSize);
200    }
201
202    private void CompileInstructions(ILGenerator il, InterpreterState state, IDataset ds) {
203      Instruction currentInstr = state.NextInstruction();
204      int nArgs = currentInstr.nArguments;
205
206      switch (currentInstr.opCode) {
207        case OpCodes.Add: {
208            if (nArgs > 0) {
209              CompileInstructions(il, state, ds);
210            }
211            for (int i = 1; i < nArgs; i++) {
212              CompileInstructions(il, state, ds);
213              il.Emit(System.Reflection.Emit.OpCodes.Add);
214            }
215            return;
216          }
217        case OpCodes.Sub: {
218            if (nArgs == 1) {
219              CompileInstructions(il, state, ds);
220              il.Emit(System.Reflection.Emit.OpCodes.Neg);
221              return;
222            }
223            if (nArgs > 0) {
224              CompileInstructions(il, state, ds);
225            }
226            for (int i = 1; i < nArgs; i++) {
227              CompileInstructions(il, state, ds);
228              il.Emit(System.Reflection.Emit.OpCodes.Sub);
229            }
230            return;
231          }
232        case OpCodes.Mul: {
233            if (nArgs > 0) {
234              CompileInstructions(il, state, ds);
235            }
236            for (int i = 1; i < nArgs; i++) {
237              CompileInstructions(il, state, ds);
238              il.Emit(System.Reflection.Emit.OpCodes.Mul);
239            }
240            return;
241          }
242        case OpCodes.Div: {
243            if (nArgs == 1) {
244              il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0);
245              CompileInstructions(il, state, ds);
246              il.Emit(System.Reflection.Emit.OpCodes.Div);
247              return;
248            }
249            if (nArgs > 0) {
250              CompileInstructions(il, state, ds);
251            }
252            for (int i = 1; i < nArgs; i++) {
253              CompileInstructions(il, state, ds);
254              il.Emit(System.Reflection.Emit.OpCodes.Div);
255            }
256            return;
257          }
258        case OpCodes.Average: {
259            CompileInstructions(il, state, ds);
260            for (int i = 1; i < nArgs; i++) {
261              CompileInstructions(il, state, ds);
262              il.Emit(System.Reflection.Emit.OpCodes.Add);
263            }
264            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, (double)nArgs);
265            il.Emit(System.Reflection.Emit.OpCodes.Div);
266            return;
267          }
268        case OpCodes.Absolute: {
269            CompileInstructions(il, state, ds);
270            il.Emit(System.Reflection.Emit.OpCodes.Call, abs);
271            return;
272          }
273        case OpCodes.Cos: {
274            CompileInstructions(il, state, ds);
275            il.Emit(System.Reflection.Emit.OpCodes.Call, cos);
276            return;
277          }
278        case OpCodes.Sin: {
279            CompileInstructions(il, state, ds);
280            il.Emit(System.Reflection.Emit.OpCodes.Call, sin);
281            return;
282          }
283        case OpCodes.Tan: {
284            CompileInstructions(il, state, ds);
285            il.Emit(System.Reflection.Emit.OpCodes.Call, tan);
286            return;
287          }
288        case OpCodes.Tanh: {
289            CompileInstructions(il, state, ds);
290            il.Emit(System.Reflection.Emit.OpCodes.Call, tanh);
291            return;
292          }
293        case OpCodes.Power: {
294            CompileInstructions(il, state, ds);
295            CompileInstructions(il, state, ds);
296            il.Emit(System.Reflection.Emit.OpCodes.Call, round);
297            il.Emit(System.Reflection.Emit.OpCodes.Call, power);
298            return;
299          }
300        case OpCodes.Root: {
301            CompileInstructions(il, state, ds);
302            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // 1 / round(...)
303            CompileInstructions(il, state, ds);
304            il.Emit(System.Reflection.Emit.OpCodes.Call, round);
305            il.Emit(System.Reflection.Emit.OpCodes.Div);
306            il.Emit(System.Reflection.Emit.OpCodes.Call, power);
307            return;
308          }
309        case OpCodes.Exp: {
310            CompileInstructions(il, state, ds);
311            il.Emit(System.Reflection.Emit.OpCodes.Call, exp);
312            return;
313          }
314        case OpCodes.Log: {
315            CompileInstructions(il, state, ds);
316            il.Emit(System.Reflection.Emit.OpCodes.Call, log);
317            return;
318          }
319        case OpCodes.Square: {
320            CompileInstructions(il, state, ds);
321            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 2.0);
322            il.Emit(System.Reflection.Emit.OpCodes.Call, power);
323            return;
324          }
325        case OpCodes.Cube: {
326            CompileInstructions(il, state, ds);
327            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 3.0);
328            il.Emit(System.Reflection.Emit.OpCodes.Call, power);
329            return;
330          }
331        case OpCodes.SquareRoot: {
332            CompileInstructions(il, state, ds);
333            il.Emit(System.Reflection.Emit.OpCodes.Call, sqrt);
334            return;
335          }
336        case OpCodes.CubeRoot: {
337            CompileInstructions(il, state, ds);
338            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0 / 3.0);
339            il.Emit(System.Reflection.Emit.OpCodes.Call, power);
340            return;
341          }
342        case OpCodes.AiryA: {
343            CompileInstructions(il, state, ds);
344            il.Emit(System.Reflection.Emit.OpCodes.Call, airyA);
345            return;
346          }
347        case OpCodes.AiryB: {
348            CompileInstructions(il, state, ds);
349            il.Emit(System.Reflection.Emit.OpCodes.Call, airyB);
350            return;
351          }
352        case OpCodes.Bessel: {
353            CompileInstructions(il, state, ds);
354            il.Emit(System.Reflection.Emit.OpCodes.Call, bessel);
355            return;
356          }
357        case OpCodes.CosineIntegral: {
358            CompileInstructions(il, state, ds);
359            il.Emit(System.Reflection.Emit.OpCodes.Call, cosIntegral);
360            return;
361          }
362        case OpCodes.Dawson: {
363            CompileInstructions(il, state, ds);
364            il.Emit(System.Reflection.Emit.OpCodes.Call, dawson);
365            return;
366          }
367        case OpCodes.Erf: {
368            CompileInstructions(il, state, ds);
369            il.Emit(System.Reflection.Emit.OpCodes.Call, erf);
370            return;
371          }
372        case OpCodes.ExponentialIntegralEi: {
373            CompileInstructions(il, state, ds);
374            il.Emit(System.Reflection.Emit.OpCodes.Call, expIntegralEi);
375            return;
376          }
377        case OpCodes.FresnelCosineIntegral: {
378            CompileInstructions(il, state, ds);
379            il.Emit(System.Reflection.Emit.OpCodes.Call, fresnelCosIntegral);
380            return;
381          }
382        case OpCodes.FresnelSineIntegral: {
383            CompileInstructions(il, state, ds);
384            il.Emit(System.Reflection.Emit.OpCodes.Call, fresnelSinIntegral);
385            return;
386          }
387        case OpCodes.Gamma: {
388            CompileInstructions(il, state, ds);
389            il.Emit(System.Reflection.Emit.OpCodes.Call, gamma);
390            return;
391          }
392        case OpCodes.HyperbolicCosineIntegral: {
393            CompileInstructions(il, state, ds);
394            il.Emit(System.Reflection.Emit.OpCodes.Call, hypCosIntegral);
395            return;
396          }
397        case OpCodes.HyperbolicSineIntegral: {
398            CompileInstructions(il, state, ds);
399            il.Emit(System.Reflection.Emit.OpCodes.Call, hypSinIntegral);
400            return;
401          }
402        case OpCodes.Norm: {
403            CompileInstructions(il, state, ds);
404            il.Emit(System.Reflection.Emit.OpCodes.Call, norm);
405            return;
406          }
407        case OpCodes.Psi: {
408            CompileInstructions(il, state, ds);
409            il.Emit(System.Reflection.Emit.OpCodes.Call, psi);
410            return;
411          }
412        case OpCodes.SineIntegral: {
413            CompileInstructions(il, state, ds);
414            il.Emit(System.Reflection.Emit.OpCodes.Call, sinIntegral);
415            return;
416          }
417        case OpCodes.AnalyticQuotient: {
418            CompileInstructions(il, state, ds); // x1
419            CompileInstructions(il, state, ds); // x2
420
421            il.Emit(System.Reflection.Emit.OpCodes.Dup);
422            il.Emit(System.Reflection.Emit.OpCodes.Mul); // x2*x2
423            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0);
424            il.Emit(System.Reflection.Emit.OpCodes.Mul); // 1+x2*x2
425            il.Emit(System.Reflection.Emit.OpCodes.Call, sqrt);
426            il.Emit(System.Reflection.Emit.OpCodes.Div);
427            return;
428          }
429        case OpCodes.IfThenElse: {
430            Label end = il.DefineLabel();
431            Label c1 = il.DefineLabel();
432            CompileInstructions(il, state, ds);
433            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 0.0); // > 0
434            il.Emit(System.Reflection.Emit.OpCodes.Cgt);
435            il.Emit(System.Reflection.Emit.OpCodes.Brfalse, c1);
436            CompileInstructions(il, state, ds);
437            il.Emit(System.Reflection.Emit.OpCodes.Br, end);
438            il.MarkLabel(c1);
439            CompileInstructions(il, state, ds);
440            il.MarkLabel(end);
441            return;
442          }
443        case OpCodes.AND: {
444            Label falseBranch = il.DefineLabel();
445            Label end = il.DefineLabel();
446            CompileInstructions(il, state, ds);
447            for (int i = 1; i < nArgs; i++) {
448              il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 0.0); // > 0
449              il.Emit(System.Reflection.Emit.OpCodes.Cgt);
450              il.Emit(System.Reflection.Emit.OpCodes.Brfalse, falseBranch);
451              CompileInstructions(il, state, ds);
452            }
453            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 0.0); // > 0
454            il.Emit(System.Reflection.Emit.OpCodes.Cgt);
455            il.Emit(System.Reflection.Emit.OpCodes.Brfalse, falseBranch);
456            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // 1
457            il.Emit(System.Reflection.Emit.OpCodes.Br, end);
458            il.MarkLabel(falseBranch);
459            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // -1
460            il.Emit(System.Reflection.Emit.OpCodes.Neg);
461            il.MarkLabel(end);
462            return;
463          }
464        case OpCodes.OR: {
465            Label trueBranch = il.DefineLabel();
466            Label end = il.DefineLabel();
467            Label resultBranch = il.DefineLabel();
468            CompileInstructions(il, state, ds);
469            for (int i = 1; i < nArgs; i++) {
470              Label nextArgBranch = il.DefineLabel();
471              // complex definition because of special properties of NaN 
472              il.Emit(System.Reflection.Emit.OpCodes.Dup);
473              il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 0.0); // <= 0       
474              il.Emit(System.Reflection.Emit.OpCodes.Ble, nextArgBranch);
475              il.Emit(System.Reflection.Emit.OpCodes.Br, resultBranch);
476              il.MarkLabel(nextArgBranch);
477              il.Emit(System.Reflection.Emit.OpCodes.Pop);
478              CompileInstructions(il, state, ds);
479            }
480            il.MarkLabel(resultBranch);
481            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 0.0); // > 0
482            il.Emit(System.Reflection.Emit.OpCodes.Cgt);
483            il.Emit(System.Reflection.Emit.OpCodes.Brtrue, trueBranch);
484            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // -1
485            il.Emit(System.Reflection.Emit.OpCodes.Neg);
486            il.Emit(System.Reflection.Emit.OpCodes.Br, end);
487            il.MarkLabel(trueBranch);
488            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // 1
489            il.MarkLabel(end);
490            return;
491          }
492        case OpCodes.NOT: {
493            CompileInstructions(il, state, ds);
494            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 0.0); // > 0
495            il.Emit(System.Reflection.Emit.OpCodes.Cgt);
496            il.Emit(System.Reflection.Emit.OpCodes.Conv_R8); // convert to float64
497            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 2.0); // * 2
498            il.Emit(System.Reflection.Emit.OpCodes.Mul);
499            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // - 1
500            il.Emit(System.Reflection.Emit.OpCodes.Sub);
501            il.Emit(System.Reflection.Emit.OpCodes.Neg); // * -1
502            return;
503          }
504        case OpCodes.XOR: {
505            CompileInstructions(il, state, ds);
506            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 0.0);
507            il.Emit(System.Reflection.Emit.OpCodes.Cgt);// > 0
508
509            for (int i = 1; i < nArgs; i++) {
510              CompileInstructions(il, state, ds);
511              il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 0.0);
512              il.Emit(System.Reflection.Emit.OpCodes.Cgt);// > 0
513              il.Emit(System.Reflection.Emit.OpCodes.Xor);
514            }
515            il.Emit(System.Reflection.Emit.OpCodes.Conv_R8); // convert to float64
516
517            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 2.0); // * 2
518            il.Emit(System.Reflection.Emit.OpCodes.Mul);
519            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // - 1
520            il.Emit(System.Reflection.Emit.OpCodes.Sub);
521            return;
522          }
523        case OpCodes.GT: {
524            CompileInstructions(il, state, ds);
525            CompileInstructions(il, state, ds);
526
527            il.Emit(System.Reflection.Emit.OpCodes.Cgt); // 1 (>) / 0 (otherwise)
528            il.Emit(System.Reflection.Emit.OpCodes.Conv_R8); // convert to float64
529            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 2.0); // * 2
530            il.Emit(System.Reflection.Emit.OpCodes.Mul);
531            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // - 1
532            il.Emit(System.Reflection.Emit.OpCodes.Sub);
533            return;
534          }
535        case OpCodes.LT: {
536            CompileInstructions(il, state, ds);
537            CompileInstructions(il, state, ds);
538            il.Emit(System.Reflection.Emit.OpCodes.Clt);
539            il.Emit(System.Reflection.Emit.OpCodes.Conv_R8); // convert to float64
540            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 2.0); // * 2
541            il.Emit(System.Reflection.Emit.OpCodes.Mul);
542            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 1.0); // - 1
543            il.Emit(System.Reflection.Emit.OpCodes.Sub);
544            return;
545          }
546        case OpCodes.TimeLag: {
547            LaggedTreeNode laggedTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
548            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row -= lag
549            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, laggedTreeNode.Lag);
550            il.Emit(System.Reflection.Emit.OpCodes.Add);
551            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
552            var prevLaggedContext = state.InLaggedContext;
553            state.InLaggedContext = true;
554            CompileInstructions(il, state, ds);
555            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row += lag
556            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, laggedTreeNode.Lag);
557            il.Emit(System.Reflection.Emit.OpCodes.Sub);
558            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
559            state.InLaggedContext = prevLaggedContext;
560            return;
561          }
562        case OpCodes.Integral: {
563            int savedPc = state.ProgramCounter;
564            LaggedTreeNode laggedTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
565            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row -= lag
566            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, laggedTreeNode.Lag);
567            il.Emit(System.Reflection.Emit.OpCodes.Add);
568            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
569            var prevLaggedContext = state.InLaggedContext;
570            state.InLaggedContext = true;
571            CompileInstructions(il, state, ds);
572            for (int l = laggedTreeNode.Lag; l < 0; l++) {
573              il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row += lag
574              il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_1);
575              il.Emit(System.Reflection.Emit.OpCodes.Add);
576              il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
577              state.ProgramCounter = savedPc;
578              CompileInstructions(il, state, ds);
579              il.Emit(System.Reflection.Emit.OpCodes.Add);
580            }
581            state.InLaggedContext = prevLaggedContext;
582            return;
583          }
584
585        //mkommend: derivate calculation taken from:
586        //http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
587        //one sided smooth differentiatior, N = 4
588        // y' = 1/8h (f_i + 2f_i-1, -2 f_i-3 - f_i-4)
589        case OpCodes.Derivative: {
590            int savedPc = state.ProgramCounter;
591            CompileInstructions(il, state, ds);
592            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row --
593            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_M1);
594            il.Emit(System.Reflection.Emit.OpCodes.Add);
595            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
596            state.ProgramCounter = savedPc;
597            var prevLaggedContext = state.InLaggedContext;
598            state.InLaggedContext = true;
599            CompileInstructions(il, state, ds);
600            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 2.0); // f_0 + 2 * f_1
601            il.Emit(System.Reflection.Emit.OpCodes.Mul);
602            il.Emit(System.Reflection.Emit.OpCodes.Add);
603
604            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row -=2
605            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_2);
606            il.Emit(System.Reflection.Emit.OpCodes.Sub);
607            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
608            state.ProgramCounter = savedPc;
609            CompileInstructions(il, state, ds);
610            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 2.0); // f_0 + 2 * f_1 - 2 * f_3
611            il.Emit(System.Reflection.Emit.OpCodes.Mul);
612            il.Emit(System.Reflection.Emit.OpCodes.Sub);
613
614            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row --
615            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_M1);
616            il.Emit(System.Reflection.Emit.OpCodes.Add);
617            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
618            state.ProgramCounter = savedPc;
619            CompileInstructions(il, state, ds);
620            il.Emit(System.Reflection.Emit.OpCodes.Sub); // f_0 + 2 * f_1 - 2 * f_3 - f_4
621            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, 8.0); // / 8
622            il.Emit(System.Reflection.Emit.OpCodes.Div);
623
624            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // row +=4
625            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_4);
626            il.Emit(System.Reflection.Emit.OpCodes.Add);
627            il.Emit(System.Reflection.Emit.OpCodes.Starg, 0);
628            state.InLaggedContext = prevLaggedContext;
629            return;
630          }
631        case OpCodes.Call: {
632            throw new NotSupportedException(
633              "Automatically defined functions are not supported by the SymbolicDataAnalysisTreeILEmittingInterpreter. Either turn of ADFs or change the interpeter.");
634          }
635        case OpCodes.Arg: {
636            throw new NotSupportedException(
637              "Automatically defined functions are not supported by the SymbolicDataAnalysisTreeILEmittingInterpreter. Either turn of ADFs or change the interpeter.");
638          }
639        case OpCodes.Variable: {
640            VariableTreeNode varNode = (VariableTreeNode)currentInstr.dynamicNode;
641            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_1); // load columns array
642            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, (int)currentInstr.data);
643            // load correct column of the current variable
644            il.Emit(System.Reflection.Emit.OpCodes.Ldelem_Ref);
645            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // rowIndex
646            if (!state.InLaggedContext) {
647              il.Emit(System.Reflection.Emit.OpCodes.Call, listGetValue);
648              il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, varNode.Weight); // load weight
649              il.Emit(System.Reflection.Emit.OpCodes.Mul);
650            } else {
651              var nanResult = il.DefineLabel();
652              var normalResult = il.DefineLabel();
653              il.Emit(System.Reflection.Emit.OpCodes.Dup);
654              il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_0);
655              il.Emit(System.Reflection.Emit.OpCodes.Blt, nanResult);
656              il.Emit(System.Reflection.Emit.OpCodes.Dup);
657              il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, ds.Rows);
658              il.Emit(System.Reflection.Emit.OpCodes.Bge, nanResult);
659              il.Emit(System.Reflection.Emit.OpCodes.Call, listGetValue);
660              il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, varNode.Weight); // load weight
661              il.Emit(System.Reflection.Emit.OpCodes.Mul);
662              il.Emit(System.Reflection.Emit.OpCodes.Br, normalResult);
663              il.MarkLabel(nanResult);
664              il.Emit(System.Reflection.Emit.OpCodes.Pop); // rowIndex
665              il.Emit(System.Reflection.Emit.OpCodes.Pop); // column reference
666              il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, double.NaN);
667              il.MarkLabel(normalResult);
668            }
669            return;
670          }
671        case OpCodes.LagVariable: {
672            var nanResult = il.DefineLabel();
673            var normalResult = il.DefineLabel();
674            LaggedVariableTreeNode varNode = (LaggedVariableTreeNode)currentInstr.dynamicNode;
675            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_1); // load columns array
676            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, (int)currentInstr.data);
677            // load correct column of the current variable
678            il.Emit(System.Reflection.Emit.OpCodes.Ldelem_Ref);
679            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, varNode.Lag); // lag
680            il.Emit(System.Reflection.Emit.OpCodes.Ldarg_0); // rowIndex
681            il.Emit(System.Reflection.Emit.OpCodes.Add); // actualRowIndex = rowIndex + sampleOffset
682            il.Emit(System.Reflection.Emit.OpCodes.Dup);
683            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4_0);
684            il.Emit(System.Reflection.Emit.OpCodes.Blt, nanResult);
685            il.Emit(System.Reflection.Emit.OpCodes.Dup);
686            il.Emit(System.Reflection.Emit.OpCodes.Ldc_I4, ds.Rows);
687            il.Emit(System.Reflection.Emit.OpCodes.Bge, nanResult);
688            il.Emit(System.Reflection.Emit.OpCodes.Call, listGetValue);
689            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, varNode.Weight); // load weight
690            il.Emit(System.Reflection.Emit.OpCodes.Mul);
691            il.Emit(System.Reflection.Emit.OpCodes.Br, normalResult);
692            il.MarkLabel(nanResult);
693            il.Emit(System.Reflection.Emit.OpCodes.Pop); // sample index
694            il.Emit(System.Reflection.Emit.OpCodes.Pop); // column reference
695            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, double.NaN);
696            il.MarkLabel(normalResult);
697            return;
698          }
699        case OpCodes.Constant: {
700            ConstantTreeNode constNode = (ConstantTreeNode)currentInstr.dynamicNode;
701            il.Emit(System.Reflection.Emit.OpCodes.Ldc_R8, constNode.Value);
702            return;
703          }
704
705        //mkommend: this symbol uses the logistic function f(x) = 1 / (1 + e^(-alpha * x) )
706        //to determine the relative amounts of the true and false branch see http://en.wikipedia.org/wiki/Logistic_function
707        case OpCodes.VariableCondition: {
708            throw new NotSupportedException("Interpretation of symbol " + currentInstr.dynamicNode.Symbol.Name +
709                                            " is not supported by the SymbolicDataAnalysisTreeILEmittingInterpreter");
710          }
711        default:
712          throw new NotSupportedException("Interpretation of symbol " + currentInstr.dynamicNode.Symbol.Name +
713                                          " is not supported by the SymbolicDataAnalysisTreeILEmittingInterpreter");
714      }
715    }
716
717    public static double AiryA(double x) {
718      if (double.IsNaN(x)) return double.NaN;
719      double ai, aip, bi, bip;
720      alglib.airy(x, out ai, out aip, out bi, out bip);
721      return ai;
722    }
723
724    public static double AiryB(double x) {
725      if (double.IsNaN(x)) return double.NaN;
726      double ai, aip, bi, bip;
727      alglib.airy(x, out ai, out aip, out bi, out bip);
728      return bi;
729    }
730    public static double Dawson(double x) {
731      if (double.IsNaN(x)) return double.NaN;
732      return alglib.dawsonintegral(x);
733    }
734
735    public static double Gamma(double x) {
736      if (double.IsNaN(x)) return double.NaN;
737      return alglib.gammafunction(x);
738    }
739
740    public static double Psi(double x) {
741      if (double.IsNaN(x)) return double.NaN;
742      else if (x <= 0 && (Math.Floor(x) - x).IsAlmost(0)) return double.NaN;
743      return alglib.psi(x);
744    }
745
746    public static double ExpIntegralEi(double x) {
747      if (double.IsNaN(x)) return double.NaN;
748      return alglib.exponentialintegralei(x);
749    }
750
751    public static double SinIntegral(double x) {
752      if (double.IsNaN(x)) return double.NaN;
753      double si, ci;
754      alglib.sinecosineintegrals(x, out si, out ci);
755      return si;
756    }
757
758    public static double CosIntegral(double x) {
759      if (double.IsNaN(x)) return double.NaN;
760      double si, ci;
761      alglib.sinecosineintegrals(x, out si, out ci);
762      return ci;
763    }
764
765    public static double HypSinIntegral(double x) {
766      if (double.IsNaN(x)) return double.NaN;
767      double shi, chi;
768      alglib.hyperbolicsinecosineintegrals(x, out shi, out chi);
769      return shi;
770    }
771
772    public static double HypCosIntegral(double x) {
773      if (double.IsNaN(x)) return double.NaN;
774      double shi, chi;
775      alglib.hyperbolicsinecosineintegrals(x, out shi, out chi);
776      return chi;
777    }
778
779    public static double FresnelCosIntegral(double x) {
780      if (double.IsNaN(x)) return double.NaN;
781      double c = 0, s = 0;
782      alglib.fresnelintegral(x, ref c, ref s);
783      return c;
784    }
785
786    public static double FresnelSinIntegral(double x) {
787      if (double.IsNaN(x)) return double.NaN;
788      double c = 0, s = 0;
789      alglib.fresnelintegral(x, ref c, ref s);
790      return s;
791    }
792
793    public static double Norm(double x) {
794      if (double.IsNaN(x)) return double.NaN;
795      return alglib.normaldistribution(x);
796    }
797
798    public static double Erf(double x) {
799      if (double.IsNaN(x)) return double.NaN;
800      return alglib.errorfunction(x);
801    }
802
803    public static double Bessel(double x) {
804      if (double.IsNaN(x)) return double.NaN;
805      return alglib.besseli0(x);
806    }
807  }
808}
Note: See TracBrowser for help on using the repository browser.