source: branches/symbreg-factors-2650/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/SymbolicDataAnalysisExpressionTreeILEmittingInterpreter.cs @ 14237

Last change on this file since 14237 was 14237, checked in by gkronber, 3 years ago

#2650: work in progress..

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