Free cookie consent management tool by TermsFeed Policy Generator

source: branches/ScatterSearch (trunk integration)/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/SymbolicDataAnalysisExpressionTreeILEmittingInterpreter.cs @ 7994

Last change on this file since 7994 was 7708, checked in by gkronber, 13 years ago

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