Free cookie consent management tool by TermsFeed Policy Generator

source: branches/2972_PDPRowSelect/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/SymbolicDataAnalysisExpressionTreeILEmittingInterpreter.cs @ 16724

Last change on this file since 16724 was 16360, checked in by gkronber, 6 years ago

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