Free cookie consent management tool by TermsFeed Policy Generator

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

Last change on this file since 7887 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: 25.2 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 HeuristicLab.Common;
25using HeuristicLab.Core;
26using HeuristicLab.Data;
27using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
28using HeuristicLab.Parameters;
29using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
30
31namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
32  [StorableClass]
33  [Item("SymbolicDataAnalysisExpressionTreeInterpreter", "Interpreter for symbolic expression trees including automatically defined functions.")]
34  public sealed class SymbolicDataAnalysisExpressionTreeInterpreter : ParameterizedNamedItem, ISymbolicDataAnalysisExpressionTreeInterpreter {
35    private const string CheckExpressionsWithIntervalArithmeticParameterName = "CheckExpressionsWithIntervalArithmetic";
36    private const string EvaluatedSolutionsParameterName = "EvaluatedSolutions";
37    #region private classes
38    private class InterpreterState {
39      private double[] argumentStack;
40      private int argumentStackPointer;
41      private Instruction[] code;
42      private int pc;
43      public int ProgramCounter {
44        get { return pc; }
45        set { pc = value; }
46      }
47      internal InterpreterState(Instruction[] code, int argumentStackSize) {
48        this.code = code;
49        this.pc = 0;
50        if (argumentStackSize > 0) {
51          this.argumentStack = new double[argumentStackSize];
52        }
53        this.argumentStackPointer = 0;
54      }
55
56      internal void Reset() {
57        this.pc = 0;
58        this.argumentStackPointer = 0;
59      }
60
61      internal Instruction NextInstruction() {
62        return code[pc++];
63      }
64      private void Push(double val) {
65        argumentStack[argumentStackPointer++] = val;
66      }
67      private double Pop() {
68        return argumentStack[--argumentStackPointer];
69      }
70
71      internal void CreateStackFrame(double[] argValues) {
72        // push in reverse order to make indexing easier
73        for (int i = argValues.Length - 1; i >= 0; i--) {
74          argumentStack[argumentStackPointer++] = argValues[i];
75        }
76        Push(argValues.Length);
77      }
78
79      internal void RemoveStackFrame() {
80        int size = (int)Pop();
81        argumentStackPointer -= size;
82      }
83
84      internal double GetStackFrameValue(ushort index) {
85        // layout of stack:
86        // [0]   <- argumentStackPointer
87        // [StackFrameSize = N + 1]
88        // [Arg0] <- argumentStackPointer - 2 - 0
89        // [Arg1] <- argumentStackPointer - 2 - 1
90        // [...]
91        // [ArgN] <- argumentStackPointer - 2 - N
92        // <Begin of stack frame>
93        return argumentStack[argumentStackPointer - index - 2];
94      }
95    }
96    private class OpCodes {
97      public const byte Add = 1;
98      public const byte Sub = 2;
99      public const byte Mul = 3;
100      public const byte Div = 4;
101
102      public const byte Sin = 5;
103      public const byte Cos = 6;
104      public const byte Tan = 7;
105
106      public const byte Log = 8;
107      public const byte Exp = 9;
108
109      public const byte IfThenElse = 10;
110
111      public const byte GT = 11;
112      public const byte LT = 12;
113
114      public const byte AND = 13;
115      public const byte OR = 14;
116      public const byte NOT = 15;
117
118
119      public const byte Average = 16;
120
121      public const byte Call = 17;
122
123      public const byte Variable = 18;
124      public const byte LagVariable = 19;
125      public const byte Constant = 20;
126      public const byte Arg = 21;
127
128      public const byte Power = 22;
129      public const byte Root = 23;
130      public const byte TimeLag = 24;
131      public const byte Integral = 25;
132      public const byte Derivative = 26;
133
134      public const byte VariableCondition = 27;
135
136      public const byte Square = 28;
137      public const byte SquareRoot = 29;
138      public const byte Gamma = 30;
139      public const byte Psi = 31;
140      public const byte Dawson = 32;
141      public const byte ExponentialIntegralEi = 33;
142      public const byte CosineIntegral = 34;
143      public const byte SineIntegral = 35;
144      public const byte HyperbolicCosineIntegral = 36;
145      public const byte HyperbolicSineIntegral = 37;
146      public const byte FresnelCosineIntegral = 38;
147      public const byte FresnelSineIntegral = 39;
148      public const byte AiryA = 40;
149      public const byte AiryB = 41;
150      public const byte Norm = 42;
151      public const byte Erf = 43;
152      public const byte Bessel = 44;
153    }
154    #endregion
155
156    private Dictionary<Type, byte> symbolToOpcode = new Dictionary<Type, byte>() {
157      { typeof(Addition), OpCodes.Add },
158      { typeof(Subtraction), OpCodes.Sub },
159      { typeof(Multiplication), OpCodes.Mul },
160      { typeof(Division), OpCodes.Div },
161      { typeof(Sine), OpCodes.Sin },
162      { typeof(Cosine), OpCodes.Cos },
163      { typeof(Tangent), OpCodes.Tan },
164      { typeof(Logarithm), OpCodes.Log },
165      { typeof(Exponential), OpCodes.Exp },
166      { typeof(IfThenElse), OpCodes.IfThenElse },
167      { typeof(GreaterThan), OpCodes.GT },
168      { typeof(LessThan), OpCodes.LT },
169      { typeof(And), OpCodes.AND },
170      { typeof(Or), OpCodes.OR },
171      { typeof(Not), OpCodes.NOT},
172      { typeof(Average), OpCodes.Average},
173      { typeof(InvokeFunction), OpCodes.Call },
174      { typeof(Variable), OpCodes.Variable },
175      { typeof(LaggedVariable), OpCodes.LagVariable },
176      { typeof(Constant), OpCodes.Constant },
177      { typeof(Argument), OpCodes.Arg },
178      { typeof(Power),OpCodes.Power},
179      { typeof(Root),OpCodes.Root},
180      { typeof(TimeLag), OpCodes.TimeLag},
181      { typeof(Integral), OpCodes.Integral},
182      { typeof(Derivative), OpCodes.Derivative},
183      { typeof(VariableCondition),OpCodes.VariableCondition},
184      { typeof(Square),OpCodes.Square},
185      { typeof(SquareRoot),OpCodes.SquareRoot},
186      { typeof(Gamma), OpCodes.Gamma },
187      { typeof(Psi), OpCodes.Psi },
188      { typeof(Dawson), OpCodes.Dawson},
189      { typeof(ExponentialIntegralEi), OpCodes.ExponentialIntegralEi },
190      { typeof(CosineIntegral), OpCodes.CosineIntegral },
191      { typeof(SineIntegral), OpCodes.SineIntegral },
192      { typeof(HyperbolicCosineIntegral), OpCodes.HyperbolicCosineIntegral },
193      { typeof(HyperbolicSineIntegral), OpCodes.HyperbolicSineIntegral },
194      { typeof(FresnelCosineIntegral), OpCodes.FresnelCosineIntegral },
195      { typeof(FresnelSineIntegral), OpCodes.FresnelSineIntegral },
196      { typeof(AiryA), OpCodes.AiryA },
197      { typeof(AiryB), OpCodes.AiryB },
198      { typeof(Norm), OpCodes.Norm},
199      { typeof(Erf), OpCodes.Erf},
200      { typeof(Bessel), OpCodes.Bessel}     
201    };
202
203    public override bool CanChangeName {
204      get { return false; }
205    }
206    public override bool CanChangeDescription {
207      get { return false; }
208    }
209
210    #region parameter properties
211    public IValueParameter<BoolValue> CheckExpressionsWithIntervalArithmeticParameter {
212      get { return (IValueParameter<BoolValue>)Parameters[CheckExpressionsWithIntervalArithmeticParameterName]; }
213    }
214
215    public IValueParameter<IntValue> EvaluatedSolutionsParameter {
216      get { return (IValueParameter<IntValue>)Parameters[EvaluatedSolutionsParameterName]; }
217    }
218    #endregion
219
220    #region properties
221    public BoolValue CheckExpressionsWithIntervalArithmetic {
222      get { return CheckExpressionsWithIntervalArithmeticParameter.Value; }
223      set { CheckExpressionsWithIntervalArithmeticParameter.Value = value; }
224    }
225
226    public IntValue EvaluatedSolutions {
227      get { return EvaluatedSolutionsParameter.Value; }
228      set { EvaluatedSolutionsParameter.Value = value; }
229    }
230    #endregion
231
232    [StorableConstructor]
233    private SymbolicDataAnalysisExpressionTreeInterpreter(bool deserializing) : base(deserializing) { }
234    private SymbolicDataAnalysisExpressionTreeInterpreter(SymbolicDataAnalysisExpressionTreeInterpreter original, Cloner cloner) : base(original, cloner) { }
235    public override IDeepCloneable Clone(Cloner cloner) {
236      return new SymbolicDataAnalysisExpressionTreeInterpreter(this, cloner);
237    }
238
239    public SymbolicDataAnalysisExpressionTreeInterpreter()
240      : base("SymbolicDataAnalysisExpressionTreeInterpreter", "Interpreter for symbolic expression trees including automatically defined functions.") {
241      Parameters.Add(new ValueParameter<BoolValue>(CheckExpressionsWithIntervalArithmeticParameterName, "Switch that determines if the interpreter checks the validity of expressions with interval arithmetic before evaluating the expression.", new BoolValue(false)));
242      Parameters.Add(new ValueParameter<IntValue>(EvaluatedSolutionsParameterName, "A counter for the total number of solutions the interpreter has evaluated", new IntValue(0)));
243    }
244
245    [StorableHook(HookType.AfterDeserialization)]
246    private void AfterDeserialization() {
247      if (!Parameters.ContainsKey(EvaluatedSolutionsParameterName))
248        Parameters.Add(new ValueParameter<IntValue>(EvaluatedSolutionsParameterName, "A counter for the total number of solutions the interpreter has evaluated", new IntValue(0)));
249    }
250
251    #region IStatefulItem
252    public void InitializeState() {
253      EvaluatedSolutions.Value = 0;
254    }
255
256    public void ClearState() {
257    }
258    #endregion
259
260    public IEnumerable<double> GetSymbolicExpressionTreeValues(ISymbolicExpressionTree tree, Dataset dataset, IEnumerable<int> rows) {
261      if (CheckExpressionsWithIntervalArithmetic.Value)
262        throw new NotSupportedException("Interval arithmetic is not yet supported in the symbolic data analysis interpreter.");
263      EvaluatedSolutions.Value++; // increment the evaluated solutions counter
264      var compiler = new SymbolicExpressionTreeCompiler();
265      Instruction[] code = compiler.Compile(tree, MapSymbolToOpCode);
266      int necessaryArgStackSize = 0;
267      for (int i = 0; i < code.Length; i++) {
268        Instruction instr = code[i];
269        if (instr.opCode == OpCodes.Variable) {
270          var variableTreeNode = instr.dynamicNode as VariableTreeNode;
271          instr.iArg0 = dataset.GetReadOnlyDoubleValues(variableTreeNode.VariableName);
272          code[i] = instr;
273        } else if (instr.opCode == OpCodes.LagVariable) {
274          var laggedVariableTreeNode = instr.dynamicNode as LaggedVariableTreeNode;
275          instr.iArg0 = dataset.GetReadOnlyDoubleValues(laggedVariableTreeNode.VariableName);
276          code[i] = instr;
277        } else if (instr.opCode == OpCodes.VariableCondition) {
278          var variableConditionTreeNode = instr.dynamicNode as VariableConditionTreeNode;
279          instr.iArg0 = dataset.GetReadOnlyDoubleValues(variableConditionTreeNode.VariableName);
280        } else if (instr.opCode == OpCodes.Call) {
281          necessaryArgStackSize += instr.nArguments + 1;
282        }
283      }
284      var state = new InterpreterState(code, necessaryArgStackSize);
285
286      foreach (var rowEnum in rows) {
287        int row = rowEnum;
288        state.Reset();
289        yield return Evaluate(dataset, ref row, state);
290      }
291    }
292
293    private double Evaluate(Dataset dataset, ref int row, InterpreterState state) {
294      Instruction currentInstr = state.NextInstruction();
295      switch (currentInstr.opCode) {
296        case OpCodes.Add: {
297            double s = Evaluate(dataset, ref row, state);
298            for (int i = 1; i < currentInstr.nArguments; i++) {
299              s += Evaluate(dataset, ref row, state);
300            }
301            return s;
302          }
303        case OpCodes.Sub: {
304            double s = Evaluate(dataset, ref row, state);
305            for (int i = 1; i < currentInstr.nArguments; i++) {
306              s -= Evaluate(dataset, ref row, state);
307            }
308            if (currentInstr.nArguments == 1) s = -s;
309            return s;
310          }
311        case OpCodes.Mul: {
312            double p = Evaluate(dataset, ref row, state);
313            for (int i = 1; i < currentInstr.nArguments; i++) {
314              p *= Evaluate(dataset, ref row, state);
315            }
316            return p;
317          }
318        case OpCodes.Div: {
319            double p = Evaluate(dataset, ref row, state);
320            for (int i = 1; i < currentInstr.nArguments; i++) {
321              p /= Evaluate(dataset, ref row, state);
322            }
323            if (currentInstr.nArguments == 1) p = 1.0 / p;
324            return p;
325          }
326        case OpCodes.Average: {
327            double sum = Evaluate(dataset, ref row, state);
328            for (int i = 1; i < currentInstr.nArguments; i++) {
329              sum += Evaluate(dataset, ref row, state);
330            }
331            return sum / currentInstr.nArguments;
332          }
333        case OpCodes.Cos: {
334            return Math.Cos(Evaluate(dataset, ref row, state));
335          }
336        case OpCodes.Sin: {
337            return Math.Sin(Evaluate(dataset, ref row, state));
338          }
339        case OpCodes.Tan: {
340            return Math.Tan(Evaluate(dataset, ref row, state));
341          }
342        case OpCodes.Square: {
343            return Math.Pow(Evaluate(dataset, ref row, state), 2);
344          }
345        case OpCodes.Power: {
346            double x = Evaluate(dataset, ref row, state);
347            double y = Math.Round(Evaluate(dataset, ref row, state));
348            return Math.Pow(x, y);
349          }
350        case OpCodes.SquareRoot: {
351            return Math.Sqrt(Evaluate(dataset, ref row, state));
352          }
353        case OpCodes.Root: {
354            double x = Evaluate(dataset, ref row, state);
355            double y = Math.Round(Evaluate(dataset, ref row, state));
356            return Math.Pow(x, 1 / y);
357          }
358        case OpCodes.Exp: {
359            return Math.Exp(Evaluate(dataset, ref row, state));
360          }
361        case OpCodes.Log: {
362            return Math.Log(Evaluate(dataset, ref row, state));
363          }
364        case OpCodes.Gamma: {
365            var x = Evaluate(dataset, ref row, state);
366            if (double.IsNaN(x)) return double.NaN;
367            else return alglib.gammafunction(x);
368          }
369        case OpCodes.Psi: {
370            var x = Evaluate(dataset, ref row, state);
371            if (double.IsNaN(x)) return double.NaN;
372            else if (x.IsAlmost(0.0)) return double.NaN;
373            else if ((Math.Floor(x) - x).IsAlmost(0)) return double.NaN;
374            return alglib.psi(x);
375          }
376        case OpCodes.Dawson: {
377            var x = Evaluate(dataset, ref row, state);
378            if (double.IsNaN(x)) return double.NaN;
379            return alglib.dawsonintegral(x);
380          }
381        case OpCodes.ExponentialIntegralEi: {
382            var x = Evaluate(dataset, ref row, state);
383            if (double.IsNaN(x)) return double.NaN;
384            return alglib.exponentialintegralei(x);
385          }
386        case OpCodes.SineIntegral: {
387            double si, ci;
388            var x = Evaluate(dataset, ref row, state);
389            if (double.IsNaN(x)) return double.NaN;
390            else {
391              alglib.sinecosineintegrals(x, out si, out ci);
392              return si;
393            }
394          }
395        case OpCodes.CosineIntegral: {
396            double si, ci;
397            var x = Evaluate(dataset, ref row, state);
398            if (double.IsNaN(x)) return double.NaN;
399            else {
400              alglib.sinecosineintegrals(x, out si, out ci);
401              return ci;
402            }
403          }
404        case OpCodes.HyperbolicSineIntegral: {
405            double shi, chi;
406            var x = Evaluate(dataset, ref row, state);
407            if (double.IsNaN(x)) return double.NaN;
408            else {
409              alglib.hyperbolicsinecosineintegrals(x, out shi, out chi);
410              return shi;
411            }
412          }
413        case OpCodes.HyperbolicCosineIntegral: {
414            double shi, chi;
415            var x = Evaluate(dataset, ref row, state);
416            if (double.IsNaN(x)) return double.NaN;
417            else {
418              alglib.hyperbolicsinecosineintegrals(x, out shi, out chi);
419              return chi;
420            }
421          }
422        case OpCodes.FresnelCosineIntegral: {
423            double c = 0, s = 0;
424            var x = Evaluate(dataset, ref row, state);
425            if (double.IsNaN(x)) return double.NaN;
426            else {
427              alglib.fresnelintegral(x, ref c, ref s);
428              return c;
429            }
430          }
431        case OpCodes.FresnelSineIntegral: {
432            double c = 0, s = 0;
433            var x = Evaluate(dataset, ref row, state);
434            if (double.IsNaN(x)) return double.NaN;
435            else {
436              alglib.fresnelintegral(x, ref c, ref s);
437              return s;
438            }
439          }
440        case OpCodes.AiryA: {
441            double ai, aip, bi, bip;
442            var x = Evaluate(dataset, ref row, state);
443            if (double.IsNaN(x)) return double.NaN;
444            else {
445              alglib.airy(x, out ai, out aip, out bi, out bip);
446              return ai;
447            }
448          }
449        case OpCodes.AiryB: {
450            double ai, aip, bi, bip;
451            var x = Evaluate(dataset, ref row, state);
452            if (double.IsNaN(x)) return double.NaN;
453            else {
454              alglib.airy(x, out ai, out aip, out bi, out bip);
455              return bi;
456            }
457          }
458        case OpCodes.Norm: {
459            var x = Evaluate(dataset, ref row, state);
460            if (double.IsNaN(x)) return double.NaN;
461            else return alglib.normaldistribution(x);
462          }
463        case OpCodes.Erf: {
464            var x = Evaluate(dataset, ref row, state);
465            if (double.IsNaN(x)) return double.NaN;
466            else return alglib.errorfunction(x);
467          }
468        case OpCodes.Bessel: {
469            var x = Evaluate(dataset, ref row, state);
470            if (double.IsNaN(x)) return double.NaN;
471            else return alglib.besseli0(x);
472          }
473        case OpCodes.IfThenElse: {
474            double condition = Evaluate(dataset, ref row, state);
475            double result;
476            if (condition > 0.0) {
477              result = Evaluate(dataset, ref row, state); SkipInstructions(state);
478            } else {
479              SkipInstructions(state); result = Evaluate(dataset, ref row, state);
480            }
481            return result;
482          }
483        case OpCodes.AND: {
484            double result = Evaluate(dataset, ref row, state);
485            for (int i = 1; i < currentInstr.nArguments; i++) {
486              if (result > 0.0) result = Evaluate(dataset, ref row, state);
487              else {
488                SkipInstructions(state);
489              }
490            }
491            return result > 0.0 ? 1.0 : -1.0;
492          }
493        case OpCodes.OR: {
494            double result = Evaluate(dataset, ref row, state);
495            for (int i = 1; i < currentInstr.nArguments; i++) {
496              if (result <= 0.0) result = Evaluate(dataset, ref row, state);
497              else {
498                SkipInstructions(state);
499              }
500            }
501            return result > 0.0 ? 1.0 : -1.0;
502          }
503        case OpCodes.NOT: {
504            return Evaluate(dataset, ref row, state) > 0.0 ? -1.0 : 1.0;
505          }
506        case OpCodes.GT: {
507            double x = Evaluate(dataset, ref row, state);
508            double y = Evaluate(dataset, ref row, state);
509            if (x > y) return 1.0;
510            else return -1.0;
511          }
512        case OpCodes.LT: {
513            double x = Evaluate(dataset, ref row, state);
514            double y = Evaluate(dataset, ref row, state);
515            if (x < y) return 1.0;
516            else return -1.0;
517          }
518        case OpCodes.TimeLag: {
519            var timeLagTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
520            row += timeLagTreeNode.Lag;
521            double result = Evaluate(dataset, ref row, state);
522            row -= timeLagTreeNode.Lag;
523            return result;
524          }
525        case OpCodes.Integral: {
526            int savedPc = state.ProgramCounter;
527            var timeLagTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
528            double sum = 0.0;
529            for (int i = 0; i < Math.Abs(timeLagTreeNode.Lag); i++) {
530              row += Math.Sign(timeLagTreeNode.Lag);
531              sum += Evaluate(dataset, ref row, state);
532              state.ProgramCounter = savedPc;
533            }
534            row -= timeLagTreeNode.Lag;
535            sum += Evaluate(dataset, ref row, state);
536            return sum;
537          }
538
539        //mkommend: derivate calculation taken from:
540        //http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
541        //one sided smooth differentiatior, N = 4
542        // y' = 1/8h (f_i + 2f_i-1, -2 f_i-3 - f_i-4)
543        case OpCodes.Derivative: {
544            int savedPc = state.ProgramCounter;
545            double f_0 = Evaluate(dataset, ref row, state); row--;
546            state.ProgramCounter = savedPc;
547            double f_1 = Evaluate(dataset, ref row, state); row -= 2;
548            state.ProgramCounter = savedPc;
549            double f_3 = Evaluate(dataset, ref row, state); row--;
550            state.ProgramCounter = savedPc;
551            double f_4 = Evaluate(dataset, ref row, state);
552            row += 4;
553
554            return (f_0 + 2 * f_1 - 2 * f_3 - f_4) / 8; // h = 1
555          }
556        case OpCodes.Call: {
557            // evaluate sub-trees
558            double[] argValues = new double[currentInstr.nArguments];
559            for (int i = 0; i < currentInstr.nArguments; i++) {
560              argValues[i] = Evaluate(dataset, ref row, state);
561            }
562            // push on argument values on stack
563            state.CreateStackFrame(argValues);
564
565            // save the pc
566            int savedPc = state.ProgramCounter;
567            // set pc to start of function 
568            state.ProgramCounter = (ushort)currentInstr.iArg0;
569            // evaluate the function
570            double v = Evaluate(dataset, ref row, state);
571
572            // delete the stack frame
573            state.RemoveStackFrame();
574
575            // restore the pc => evaluation will continue at point after my subtrees 
576            state.ProgramCounter = savedPc;
577            return v;
578          }
579        case OpCodes.Arg: {
580            return state.GetStackFrameValue((ushort)currentInstr.iArg0);
581          }
582        case OpCodes.Variable: {
583            if (row < 0 || row >= dataset.Rows)
584              return double.NaN;
585            var variableTreeNode = (VariableTreeNode)currentInstr.dynamicNode;
586            return ((IList<double>)currentInstr.iArg0)[row] * variableTreeNode.Weight;
587          }
588        case OpCodes.LagVariable: {
589            var laggedVariableTreeNode = (LaggedVariableTreeNode)currentInstr.dynamicNode;
590            int actualRow = row + laggedVariableTreeNode.Lag;
591            if (actualRow < 0 || actualRow >= dataset.Rows)
592              return double.NaN;
593            return ((IList<double>)currentInstr.iArg0)[actualRow] * laggedVariableTreeNode.Weight;
594          }
595        case OpCodes.Constant: {
596            var constTreeNode = currentInstr.dynamicNode as ConstantTreeNode;
597            return constTreeNode.Value;
598          }
599
600        //mkommend: this symbol uses the logistic function f(x) = 1 / (1 + e^(-alpha * x) )
601        //to determine the relative amounts of the true and false branch see http://en.wikipedia.org/wiki/Logistic_function
602        case OpCodes.VariableCondition: {
603            if (row < 0 || row >= dataset.Rows)
604              return double.NaN;
605            var variableConditionTreeNode = (VariableConditionTreeNode)currentInstr.dynamicNode;
606            double variableValue = ((IList<double>)currentInstr.iArg0)[row];
607            double x = variableValue - variableConditionTreeNode.Threshold;
608            double p = 1 / (1 + Math.Exp(-variableConditionTreeNode.Slope * x));
609
610            double trueBranch = Evaluate(dataset, ref row, state);
611            double falseBranch = Evaluate(dataset, ref row, state);
612
613            return trueBranch * p + falseBranch * (1 - p);
614          }
615        default: throw new NotSupportedException();
616      }
617    }
618
619    private byte MapSymbolToOpCode(ISymbolicExpressionTreeNode treeNode) {
620      if (symbolToOpcode.ContainsKey(treeNode.Symbol.GetType()))
621        return symbolToOpcode[treeNode.Symbol.GetType()];
622      else
623        throw new NotSupportedException("Symbol: " + treeNode.Symbol);
624    }
625
626    // skips a whole branch
627    private void SkipInstructions(InterpreterState state) {
628      int i = 1;
629      while (i > 0) {
630        i += state.NextInstruction().nArguments;
631        i--;
632      }
633    }
634  }
635}
Note: See TracBrowser for help on using the repository browser.