Free cookie consent management tool by TermsFeed Policy Generator

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

Last change on this file since 13248 was 13248, checked in by mkommend, 8 years ago

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