Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
03/04/13 14:43:44 (11 years ago)
Author:
bburlacu
Message:

#2021: Initial implementation of the SymbolicExpressionTreeLinearInterpreter.

Location:
branches/HeuristicLab.DataAnalysis.Symbolic.LinearInterpreter
Files:
1 added
2 copied

Legend:

Unmodified
Added
Removed

  • branches/HeuristicLab.DataAnalysis.Symbolic.LinearInterpreter/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/SymbolicDataAnalysisExpressionTreeLinearInterpreter.cs

    r9257 r9271  
    3131namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
    3232  [StorableClass]
    33   [Item("SymbolicDataAnalysisExpressionTreeInterpreter", "Interpreter for symbolic expression trees including automatically defined functions.")]
    34   public class SymbolicDataAnalysisExpressionTreeInterpreter : ParameterizedNamedItem, ISymbolicDataAnalysisExpressionTreeInterpreter {
     33  [Item("SymbolicDataAnalysisExpressionTreeLinearInterpreter", "Interpreter for symbolic expression trees including automatically defined functions.")]
     34  public class SymbolicDataAnalysisExpressionTreeLinearInterpreter : ParameterizedNamedItem, ISymbolicDataAnalysisExpressionTreeInterpreter {
    3535    private const string CheckExpressionsWithIntervalArithmeticParameterName = "CheckExpressionsWithIntervalArithmetic";
    3636    private const string EvaluatedSolutionsParameterName = "EvaluatedSolutions";
    … …  
    6262
    6363    [StorableConstructor]
    64     protected SymbolicDataAnalysisExpressionTreeInterpreter(bool deserializing) : base(deserializing) { }
    65     protected SymbolicDataAnalysisExpressionTreeInterpreter(SymbolicDataAnalysisExpressionTreeInterpreter original, Cloner cloner) : base(original, cloner) { }
     64    protected SymbolicDataAnalysisExpressionTreeLinearInterpreter(bool deserializing) : base(deserializing) { }
     65    protected SymbolicDataAnalysisExpressionTreeLinearInterpreter(SymbolicDataAnalysisExpressionTreeLinearInterpreter original, Cloner cloner) : base(original, cloner) { }
    6666    public override IDeepCloneable Clone(Cloner cloner) {
    67       return new SymbolicDataAnalysisExpressionTreeInterpreter(this, cloner);
    68     }
    69 
    70     public SymbolicDataAnalysisExpressionTreeInterpreter()
    71       : base("SymbolicDataAnalysisExpressionTreeInterpreter", "Interpreter for symbolic expression trees including automatically defined functions.") {
     67      return new SymbolicDataAnalysisExpressionTreeLinearInterpreter(this, cloner);
     68    }
     69
     70    public SymbolicDataAnalysisExpressionTreeLinearInterpreter()
     71      : base("SymbolicDataAnalysisExpressionTreeLinearInterpreter", "Interpreter for symbolic expression trees including automatically defined functions.") {
    7272      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)));
    7373      Parameters.Add(new ValueParameter<IntValue>(EvaluatedSolutionsParameterName, "A counter for the total number of solutions the interpreter has evaluated", new IntValue(0)));
    7474    }
    7575
    76     protected SymbolicDataAnalysisExpressionTreeInterpreter(string name, string description)
     76    protected SymbolicDataAnalysisExpressionTreeLinearInterpreter(string name, string description)
    7777      : base(name, description) {
    7878      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)));
    … …  
    102102        EvaluatedSolutions.Value++; // increment the evaluated solutions counter
    103103      }
    104       var state = PrepareInterpreterState(tree, dataset);
     104
     105      Instruction[] code = SymbolicExpressionTreeCompiler.CompilePostfix(tree, OpCodes.MapSymbolToOpCode);
     106      foreach (Instruction instr in code) {
     107        switch (instr.opCode) {
     108          case OpCodes.Variable: {
     109              var variableTreeNode = (VariableTreeNode)instr.dynamicNode;
     110              instr.iArg0 = dataset.GetReadOnlyDoubleValues(variableTreeNode.VariableName);
     111            }
     112            break;
     113          case OpCodes.LagVariable: {
     114              var laggedVariableTreeNode = (LaggedVariableTreeNode)instr.dynamicNode;
     115              instr.iArg0 = dataset.GetReadOnlyDoubleValues(laggedVariableTreeNode.VariableName);
     116            }
     117            break;
     118          case OpCodes.VariableCondition: {
     119              var variableConditionTreeNode = (VariableConditionTreeNode)instr.dynamicNode;
     120              instr.iArg0 = dataset.GetReadOnlyDoubleValues(variableConditionTreeNode.VariableName);
     121            }
     122            break;
     123        }
     124      }
    105125
    106126      foreach (var rowEnum in rows) {
    107127        int row = rowEnum;
    108         yield return Evaluate(dataset, ref row, state);
    109         state.Reset();
     128        yield return Evaluate(dataset, ref row, code);
    110129      }
    111130    }
    112131
    113     private static InterpreterState PrepareInterpreterState(ISymbolicExpressionTree tree, Dataset dataset) {
    114       Instruction[] code = SymbolicExpressionTreeCompiler.Compile(tree, OpCodes.MapSymbolToOpCode);
    115       int necessaryArgStackSize = 0;
    116       foreach (Instruction instr in code) {
    117         if (instr.opCode == OpCodes.Variable) {
    118           var variableTreeNode = (VariableTreeNode)instr.dynamicNode;
    119           instr.iArg0 = dataset.GetReadOnlyDoubleValues(variableTreeNode.VariableName);
    120         } else if (instr.opCode == OpCodes.LagVariable) {
    121           var laggedVariableTreeNode = (LaggedVariableTreeNode)instr.dynamicNode;
    122           instr.iArg0 = dataset.GetReadOnlyDoubleValues(laggedVariableTreeNode.VariableName);
    123         } else if (instr.opCode == OpCodes.VariableCondition) {
    124           var variableConditionTreeNode = (VariableConditionTreeNode)instr.dynamicNode;
    125           instr.iArg0 = dataset.GetReadOnlyDoubleValues(variableConditionTreeNode.VariableName);
    126         } else if (instr.opCode == OpCodes.Call) {
    127           necessaryArgStackSize += instr.nArguments + 1;
     132    protected virtual double Evaluate(Dataset dataset, ref int row, Instruction[] code) {
     133      int count = 0;
     134      for (int j = 0; j != code.Length; ++j) {
     135        Instruction currentInstr = code[j];
     136
     137        switch (currentInstr.opCode) {
     138          case OpCodes.Add: {
     139              double s = code[j - currentInstr.nArguments].value;
     140              for (int i = j - currentInstr.nArguments + 1; i < j; i++) {
     141                s += code[i].value;
     142              }
     143              currentInstr.value = s;
     144            }
     145            break;
     146          case OpCodes.Sub: {
     147              double s = code[j - currentInstr.nArguments].value;
     148              for (int i = j - currentInstr.nArguments + 1; i < j; i++) {
     149                s -= code[i].value;
     150              }
     151              if (currentInstr.nArguments == 1) s = -s;
     152              currentInstr.value = s;
     153            }
     154            break;
     155          case OpCodes.Mul: {
     156              double p = code[j - currentInstr.nArguments].value;
     157              for (int i = j - currentInstr.nArguments + 1; i < j; i++) {
     158                p *= code[i].value;
     159              }
     160              currentInstr.value = p;
     161            }
     162            break;
     163          case OpCodes.Div: {
     164              double p = code[j - currentInstr.nArguments].value;
     165              for (int i = j - currentInstr.nArguments + 1; i < j; i++) {
     166                p /= code[i].value;
     167              }
     168              if (currentInstr.nArguments == 1) p = 1.0 / p;
     169              currentInstr.value = p;
     170            }
     171            break;
     172          case OpCodes.Average: {
     173              double sum = code[j - currentInstr.nArguments].value;
     174              for (int i = j - currentInstr.nArguments + 1; i < j; i++) {
     175                sum += code[i].value;
     176              }
     177              currentInstr.value = sum / currentInstr.nArguments;
     178            }
     179            break;
     180          case OpCodes.Cos: {
     181              currentInstr.value = Math.Cos(code[j - 1].value);
     182            }
     183            break;
     184          case OpCodes.Sin: {
     185              currentInstr.value = Math.Sin(code[j - 1].value);
     186              break;
     187            }
     188          case OpCodes.Tan: {
     189              currentInstr.value = Math.Tan(code[j - 1].value);
     190            }
     191            break;
     192          case OpCodes.Square: {
     193              currentInstr.value = Math.Pow(code[j - 1].value, 2);
     194            }
     195            break;
     196          case OpCodes.Power: {
     197              double x = code[j - 2].value;
     198              double y = Math.Round(code[j - 1].value);
     199              currentInstr.value = Math.Pow(x, y);
     200            }
     201            break;
     202          case OpCodes.SquareRoot: {
     203              currentInstr.value = Math.Sqrt(code[j - 1].value);
     204            }
     205            break;
     206          case OpCodes.Root: {
     207              double x = code[j - 2].value;
     208              double y = Math.Round(code[j - 1].value);
     209              currentInstr.value = Math.Pow(x, 1 / y);
     210            }
     211            break;
     212          case OpCodes.Exp: {
     213              currentInstr.value = Math.Exp(code[j - 1].value);
     214            }
     215            break;
     216          case OpCodes.Log: {
     217              currentInstr.value = Math.Log(code[j - 1].value);
     218            }
     219            break;
     220          case OpCodes.Gamma: {
     221              currentInstr.value = double.IsNaN(code[j - 1].value) ? double.NaN : alglib.gammafunction(code[j - 1].value);
     222            }
     223            break;
     224          case OpCodes.Psi: {
     225              var x = code[j - 1].value;
     226              if (double.IsNaN(x)) currentInstr.value = double.NaN;
     227              else if (x <= 0 && (Math.Floor(x) - x).IsAlmost(0)) currentInstr.value = double.NaN;
     228              else currentInstr.value = alglib.psi(x);
     229            }
     230            break;
     231          case OpCodes.Dawson: {
     232              var x = code[j - 1].value;
     233              currentInstr.value = double.IsNaN(x) ? double.NaN : alglib.dawsonintegral(x);
     234            }
     235            break;
     236          case OpCodes.ExponentialIntegralEi: {
     237              var x = code[j - 1].value;
     238              currentInstr.value = double.IsNaN(x) ? double.NaN : alglib.exponentialintegralei(x);
     239            }
     240            break;
     241          case OpCodes.SineIntegral: {
     242              double si, ci;
     243              var x = code[j - 1].value;
     244              if (double.IsNaN(x)) currentInstr.value = double.NaN;
     245              else {
     246                alglib.sinecosineintegrals(x, out si, out ci);
     247                currentInstr.value = si;
     248              }
     249            }
     250            break;
     251          case OpCodes.CosineIntegral: {
     252              double si, ci;
     253              var x = code[j - 1].value;
     254              if (double.IsNaN(x)) currentInstr.value = double.NaN;
     255              else {
     256                alglib.sinecosineintegrals(x, out si, out ci);
     257                currentInstr.value = ci;
     258              }
     259            }
     260            break;
     261          case OpCodes.HyperbolicSineIntegral: {
     262              double shi, chi;
     263              var x = code[j - 1].value;
     264              if (double.IsNaN(x)) currentInstr.value = double.NaN;
     265              else {
     266                alglib.hyperbolicsinecosineintegrals(x, out shi, out chi);
     267                currentInstr.value = shi;
     268              }
     269            }
     270            break;
     271          case OpCodes.HyperbolicCosineIntegral: {
     272              double shi, chi;
     273              var x = code[j - 1].value;
     274              if (double.IsNaN(x)) currentInstr.value = double.NaN;
     275              else {
     276                alglib.hyperbolicsinecosineintegrals(x, out shi, out chi);
     277                currentInstr.value = chi;
     278              }
     279            }
     280            break;
     281          case OpCodes.FresnelCosineIntegral: {
     282              double c = 0, s = 0;
     283              var x = code[j - 1].value;
     284              if (double.IsNaN(x)) currentInstr.value = double.NaN;
     285              else {
     286                alglib.fresnelintegral(x, ref c, ref s);
     287                currentInstr.value = c;
     288              }
     289            }
     290            break;
     291          case OpCodes.FresnelSineIntegral: {
     292              double c = 0, s = 0;
     293              var x = code[j - 1].value;
     294              if (double.IsNaN(x)) currentInstr.value = double.NaN;
     295              else {
     296                alglib.fresnelintegral(x, ref c, ref s);
     297                currentInstr.value = s;
     298              }
     299            }
     300            break;
     301          case OpCodes.AiryA: {
     302              double ai, aip, bi, bip;
     303              var x = code[j - 1].value;
     304              if (double.IsNaN(x)) currentInstr.value = double.NaN;
     305              else {
     306                alglib.airy(x, out ai, out aip, out bi, out bip);
     307                currentInstr.value = ai;
     308              }
     309            }
     310            break;
     311          case OpCodes.AiryB: {
     312              double ai, aip, bi, bip;
     313              var x = code[j - 1].value;
     314              if (double.IsNaN(x)) currentInstr.value = double.NaN;
     315              else {
     316                alglib.airy(x, out ai, out aip, out bi, out bip);
     317                currentInstr.value = bi;
     318              }
     319            }
     320            break;
     321          case OpCodes.Norm: {
     322              var x = code[j - 1].value;
     323              currentInstr.value = double.IsNaN(x) ? double.NaN : alglib.normaldistribution(x);
     324            }
     325            break;
     326          case OpCodes.Erf: {
     327              var x = code[j - 1].value;
     328              currentInstr.value = double.IsNaN(x) ? double.NaN : alglib.errorfunction(x);
     329            }
     330            break;
     331          case OpCodes.Bessel: {
     332              var x = code[j - 1].value;
     333              currentInstr.value = double.IsNaN(x) ? double.NaN : alglib.besseli0(x);
     334            }
     335            break;
     336          case OpCodes.IfThenElse: {
     337              double condition = code[j - currentInstr.nArguments].value;
     338              double result = condition > 0.0 ? code[j - 2].value : code[j - 1].value;
     339              currentInstr.value = result;
     340            }
     341            break;
     342          case OpCodes.AND: {
     343              double result = code[j - currentInstr.nArguments].value;
     344              for (int i = j - currentInstr.nArguments + 1; i < j; i++) {
     345                if (result > 0.0) result = code[i].value;
     346              }
     347              currentInstr.value = result > 0.0 ? 1.0 : -1.0;
     348            }
     349            break;
     350          case OpCodes.OR: {
     351              double result = code[j - currentInstr.nArguments].value;
     352              for (int i = 1; i < currentInstr.nArguments; i++) {
     353                if (result <= 0.0) result = code[i].value; ;
     354              }
     355              currentInstr.value = result > 0.0 ? 1.0 : -1.0;
     356            }
     357            break;
     358          case OpCodes.NOT: {
     359              currentInstr.value = code[j - 1].value > 0.0 ? -1.0 : 1.0;
     360            }
     361            break;
     362          case OpCodes.GT: {
     363              double x = code[j - 2].value;
     364              double y = code[j - 1].value;
     365              currentInstr.value = x > y ? 1.0 : -1.0;
     366            }
     367            break;
     368          case OpCodes.LT: {
     369              double x = code[j - 2].value;
     370              double y = code[j - 1].value;
     371              currentInstr.value = x < y ? 1.0 : -1.0;
     372            }
     373            break;
     374          case OpCodes.TimeLag: {
     375              throw new NotSupportedException();
     376            }
     377          case OpCodes.Integral: {
     378              throw new NotSupportedException();
     379            }
     380
     381          //mkommend: derivate calculation taken from:
     382          //http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
     383          //one sided smooth differentiatior, N = 4
     384          // y' = 1/8h (f_i + 2f_i-1, -2 f_i-3 - f_i-4)
     385          case OpCodes.Derivative: {
     386              throw new NotSupportedException();
     387            }
     388          case OpCodes.Call: {
     389              throw new NotSupportedException();
     390            }
     391          case OpCodes.Arg: {
     392              throw new NotSupportedException();
     393            }
     394          case OpCodes.Variable: {
     395              if (row < 0 || row >= dataset.Rows) currentInstr.value = double.NaN;
     396              else {
     397                var variableTreeNode = (VariableTreeNode)currentInstr.dynamicNode;
     398                currentInstr.value = ((IList<double>)currentInstr.iArg0)[row] * variableTreeNode.Weight;
     399              }
     400            }
     401            break;
     402          case OpCodes.LagVariable: {
     403              var laggedVariableTreeNode = (LaggedVariableTreeNode)currentInstr.dynamicNode;
     404              int actualRow = row + laggedVariableTreeNode.Lag;
     405              if (actualRow < 0 || actualRow >= dataset.Rows) currentInstr.value = double.NaN;
     406              else {
     407                currentInstr.value = ((IList<double>)currentInstr.iArg0)[actualRow] * laggedVariableTreeNode.Weight;
     408              }
     409            }
     410            break;
     411          case OpCodes.Constant: {
     412              var constTreeNode = (ConstantTreeNode)currentInstr.dynamicNode;
     413              currentInstr.value = constTreeNode.Value;
     414            }
     415            break;
     416
     417          //mkommend: this symbol uses the logistic function f(x) = 1 / (1 + e^(-alpha * x) )
     418          //to determine the relative amounts of the true and false branch see http://en.wikipedia.org/wiki/Logistic_function
     419          case OpCodes.VariableCondition: {
     420              if (row < 0 || row >= dataset.Rows) currentInstr.value = double.NaN;
     421              else {
     422                var variableConditionTreeNode = (VariableConditionTreeNode)currentInstr.dynamicNode;
     423                double variableValue = ((IList<double>)currentInstr.iArg0)[row];
     424                double x = variableValue - variableConditionTreeNode.Threshold;
     425                double p = 1 / (1 + Math.Exp(-variableConditionTreeNode.Slope * x));
     426
     427                double trueBranch = code[j - currentInstr.nArguments].value;
     428                double falseBranch = code[j - currentInstr.nArguments + 1].value;
     429
     430                currentInstr.value = trueBranch * p + falseBranch * (1 - p);
     431              }
     432            }
     433            break;
     434          default:
     435            throw new NotSupportedException();
    128436        }
     437
     438        // book-keeping in order to keep all the values correct
     439        int narg = currentInstr.nArguments;
     440
     441        // we count the values before the next function is encountered
     442        if (narg == 0) { ++count; continue; }
     443
     444        if (count > narg) {
     445          for (int i = 1; i <= count - narg; ++i)
     446            code[j - i].value = code[j - i - narg].value;
     447        }
     448
     449        count -= (narg - 1); // because after being evaluated the current instruction becomes a 'value'
    129450      }
    130       return new InterpreterState(code, necessaryArgStackSize);
    131     }
    132 
    133 
    134     protected virtual double Evaluate(Dataset dataset, ref int row, InterpreterState state) {
    135       Instruction currentInstr = state.NextInstruction();
    136       switch (currentInstr.opCode) {
    137         case OpCodes.Add: {
    138             double s = Evaluate(dataset, ref row, state);
    139             for (int i = 1; i < currentInstr.nArguments; i++) {
    140               s += Evaluate(dataset, ref row, state);
    141             }
    142             return s;
    143           }
    144         case OpCodes.Sub: {
    145             double s = Evaluate(dataset, ref row, state);
    146             for (int i = 1; i < currentInstr.nArguments; i++) {
    147               s -= Evaluate(dataset, ref row, state);
    148             }
    149             if (currentInstr.nArguments == 1) s = -s;
    150             return s;
    151           }
    152         case OpCodes.Mul: {
    153             double p = Evaluate(dataset, ref row, state);
    154             for (int i = 1; i < currentInstr.nArguments; i++) {
    155               p *= Evaluate(dataset, ref row, state);
    156             }
    157             return p;
    158           }
    159         case OpCodes.Div: {
    160             double p = Evaluate(dataset, ref row, state);
    161             for (int i = 1; i < currentInstr.nArguments; i++) {
    162               p /= Evaluate(dataset, ref row, state);
    163             }
    164             if (currentInstr.nArguments == 1) p = 1.0 / p;
    165             return p;
    166           }
    167         case OpCodes.Average: {
    168             double sum = Evaluate(dataset, ref row, state);
    169             for (int i = 1; i < currentInstr.nArguments; i++) {
    170               sum += Evaluate(dataset, ref row, state);
    171             }
    172             return sum / currentInstr.nArguments;
    173           }
    174         case OpCodes.Cos: {
    175             return Math.Cos(Evaluate(dataset, ref row, state));
    176           }
    177         case OpCodes.Sin: {
    178             return Math.Sin(Evaluate(dataset, ref row, state));
    179           }
    180         case OpCodes.Tan: {
    181             return Math.Tan(Evaluate(dataset, ref row, state));
    182           }
    183         case OpCodes.Square: {
    184             return Math.Pow(Evaluate(dataset, ref row, state), 2);
    185           }
    186         case OpCodes.Power: {
    187             double x = Evaluate(dataset, ref row, state);
    188             double y = Math.Round(Evaluate(dataset, ref row, state));
    189             return Math.Pow(x, y);
    190           }
    191         case OpCodes.SquareRoot: {
    192             return Math.Sqrt(Evaluate(dataset, ref row, state));
    193           }
    194         case OpCodes.Root: {
    195             double x = Evaluate(dataset, ref row, state);
    196             double y = Math.Round(Evaluate(dataset, ref row, state));
    197             return Math.Pow(x, 1 / y);
    198           }
    199         case OpCodes.Exp: {
    200             return Math.Exp(Evaluate(dataset, ref row, state));
    201           }
    202         case OpCodes.Log: {
    203             return Math.Log(Evaluate(dataset, ref row, state));
    204           }
    205         case OpCodes.Gamma: {
    206             var x = Evaluate(dataset, ref row, state);
    207             if (double.IsNaN(x)) return double.NaN;
    208             else return alglib.gammafunction(x);
    209           }
    210         case OpCodes.Psi: {
    211             var x = Evaluate(dataset, ref row, state);
    212             if (double.IsNaN(x)) return double.NaN;
    213             else if (x <= 0 && (Math.Floor(x) - x).IsAlmost(0)) return double.NaN;
    214             return alglib.psi(x);
    215           }
    216         case OpCodes.Dawson: {
    217             var x = Evaluate(dataset, ref row, state);
    218             if (double.IsNaN(x)) return double.NaN;
    219             return alglib.dawsonintegral(x);
    220           }
    221         case OpCodes.ExponentialIntegralEi: {
    222             var x = Evaluate(dataset, ref row, state);
    223             if (double.IsNaN(x)) return double.NaN;
    224             return alglib.exponentialintegralei(x);
    225           }
    226         case OpCodes.SineIntegral: {
    227             double si, ci;
    228             var x = Evaluate(dataset, ref row, state);
    229             if (double.IsNaN(x)) return double.NaN;
    230             else {
    231               alglib.sinecosineintegrals(x, out si, out ci);
    232               return si;
    233             }
    234           }
    235         case OpCodes.CosineIntegral: {
    236             double si, ci;
    237             var x = Evaluate(dataset, ref row, state);
    238             if (double.IsNaN(x)) return double.NaN;
    239             else {
    240               alglib.sinecosineintegrals(x, out si, out ci);
    241               return ci;
    242             }
    243           }
    244         case OpCodes.HyperbolicSineIntegral: {
    245             double shi, chi;
    246             var x = Evaluate(dataset, ref row, state);
    247             if (double.IsNaN(x)) return double.NaN;
    248             else {
    249               alglib.hyperbolicsinecosineintegrals(x, out shi, out chi);
    250               return shi;
    251             }
    252           }
    253         case OpCodes.HyperbolicCosineIntegral: {
    254             double shi, chi;
    255             var x = Evaluate(dataset, ref row, state);
    256             if (double.IsNaN(x)) return double.NaN;
    257             else {
    258               alglib.hyperbolicsinecosineintegrals(x, out shi, out chi);
    259               return chi;
    260             }
    261           }
    262         case OpCodes.FresnelCosineIntegral: {
    263             double c = 0, s = 0;
    264             var x = Evaluate(dataset, ref row, state);
    265             if (double.IsNaN(x)) return double.NaN;
    266             else {
    267               alglib.fresnelintegral(x, ref c, ref s);
    268               return c;
    269             }
    270           }
    271         case OpCodes.FresnelSineIntegral: {
    272             double c = 0, s = 0;
    273             var x = Evaluate(dataset, ref row, state);
    274             if (double.IsNaN(x)) return double.NaN;
    275             else {
    276               alglib.fresnelintegral(x, ref c, ref s);
    277               return s;
    278             }
    279           }
    280         case OpCodes.AiryA: {
    281             double ai, aip, bi, bip;
    282             var x = Evaluate(dataset, ref row, state);
    283             if (double.IsNaN(x)) return double.NaN;
    284             else {
    285               alglib.airy(x, out ai, out aip, out bi, out bip);
    286               return ai;
    287             }
    288           }
    289         case OpCodes.AiryB: {
    290             double ai, aip, bi, bip;
    291             var x = Evaluate(dataset, ref row, state);
    292             if (double.IsNaN(x)) return double.NaN;
    293             else {
    294               alglib.airy(x, out ai, out aip, out bi, out bip);
    295               return bi;
    296             }
    297           }
    298         case OpCodes.Norm: {
    299             var x = Evaluate(dataset, ref row, state);
    300             if (double.IsNaN(x)) return double.NaN;
    301             else return alglib.normaldistribution(x);
    302           }
    303         case OpCodes.Erf: {
    304             var x = Evaluate(dataset, ref row, state);
    305             if (double.IsNaN(x)) return double.NaN;
    306             else return alglib.errorfunction(x);
    307           }
    308         case OpCodes.Bessel: {
    309             var x = Evaluate(dataset, ref row, state);
    310             if (double.IsNaN(x)) return double.NaN;
    311             else return alglib.besseli0(x);
    312           }
    313         case OpCodes.IfThenElse: {
    314             double condition = Evaluate(dataset, ref row, state);
    315             double result;
    316             if (condition > 0.0) {
    317               result = Evaluate(dataset, ref row, state); state.SkipInstructions();
    318             } else {
    319               state.SkipInstructions(); result = Evaluate(dataset, ref row, state);
    320             }
    321             return result;
    322           }
    323         case OpCodes.AND: {
    324             double result = Evaluate(dataset, ref row, state);
    325             for (int i = 1; i < currentInstr.nArguments; i++) {
    326               if (result > 0.0) result = Evaluate(dataset, ref row, state);
    327               else {
    328                 state.SkipInstructions();
    329               }
    330             }
    331             return result > 0.0 ? 1.0 : -1.0;
    332           }
    333         case OpCodes.OR: {
    334             double result = Evaluate(dataset, ref row, state);
    335             for (int i = 1; i < currentInstr.nArguments; i++) {
    336               if (result <= 0.0) result = Evaluate(dataset, ref row, state);
    337               else {
    338                 state.SkipInstructions();
    339               }
    340             }
    341             return result > 0.0 ? 1.0 : -1.0;
    342           }
    343         case OpCodes.NOT: {
    344             return Evaluate(dataset, ref row, state) > 0.0 ? -1.0 : 1.0;
    345           }
    346         case OpCodes.GT: {
    347             double x = Evaluate(dataset, ref row, state);
    348             double y = Evaluate(dataset, ref row, state);
    349             if (x > y) return 1.0;
    350             else return -1.0;
    351           }
    352         case OpCodes.LT: {
    353             double x = Evaluate(dataset, ref row, state);
    354             double y = Evaluate(dataset, ref row, state);
    355             if (x < y) return 1.0;
    356             else return -1.0;
    357           }
    358         case OpCodes.TimeLag: {
    359             var timeLagTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
    360             row += timeLagTreeNode.Lag;
    361             double result = Evaluate(dataset, ref row, state);
    362             row -= timeLagTreeNode.Lag;
    363             return result;
    364           }
    365         case OpCodes.Integral: {
    366             int savedPc = state.ProgramCounter;
    367             var timeLagTreeNode = (LaggedTreeNode)currentInstr.dynamicNode;
    368             double sum = 0.0;
    369             for (int i = 0; i < Math.Abs(timeLagTreeNode.Lag); i++) {
    370               row += Math.Sign(timeLagTreeNode.Lag);
    371               sum += Evaluate(dataset, ref row, state);
    372               state.ProgramCounter = savedPc;
    373             }
    374             row -= timeLagTreeNode.Lag;
    375             sum += Evaluate(dataset, ref row, state);
    376             return sum;
    377           }
    378 
    379         //mkommend: derivate calculation taken from:
    380         //http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
    381         //one sided smooth differentiatior, N = 4
    382         // y' = 1/8h (f_i + 2f_i-1, -2 f_i-3 - f_i-4)
    383         case OpCodes.Derivative: {
    384             int savedPc = state.ProgramCounter;
    385             double f_0 = Evaluate(dataset, ref row, state); row--;
    386             state.ProgramCounter = savedPc;
    387             double f_1 = Evaluate(dataset, ref row, state); row -= 2;
    388             state.ProgramCounter = savedPc;
    389             double f_3 = Evaluate(dataset, ref row, state); row--;
    390             state.ProgramCounter = savedPc;
    391             double f_4 = Evaluate(dataset, ref row, state);
    392             row += 4;
    393 
    394             return (f_0 + 2 * f_1 - 2 * f_3 - f_4) / 8; // h = 1
    395           }
    396         case OpCodes.Call: {
    397             // evaluate sub-trees
    398             double[] argValues = new double[currentInstr.nArguments];
    399             for (int i = 0; i < currentInstr.nArguments; i++) {
    400               argValues[i] = Evaluate(dataset, ref row, state);
    401             }
    402             // push on argument values on stack
    403             state.CreateStackFrame(argValues);
    404 
    405             // save the pc
    406             int savedPc = state.ProgramCounter;
    407             // set pc to start of function 
    408             state.ProgramCounter = (ushort)currentInstr.iArg0;
    409             // evaluate the function
    410             double v = Evaluate(dataset, ref row, state);
    411 
    412             // delete the stack frame
    413             state.RemoveStackFrame();
    414 
    415             // restore the pc => evaluation will continue at point after my subtrees 
    416             state.ProgramCounter = savedPc;
    417             return v;
    418           }
    419         case OpCodes.Arg: {
    420             return state.GetStackFrameValue((ushort)currentInstr.iArg0);
    421           }
    422         case OpCodes.Variable: {
    423             if (row < 0 || row >= dataset.Rows) return double.NaN;
    424             var variableTreeNode = (VariableTreeNode)currentInstr.dynamicNode;
    425             return ((IList<double>)currentInstr.iArg0)[row] * variableTreeNode.Weight;
    426           }
    427         case OpCodes.LagVariable: {
    428             var laggedVariableTreeNode = (LaggedVariableTreeNode)currentInstr.dynamicNode;
    429             int actualRow = row + laggedVariableTreeNode.Lag;
    430             if (actualRow < 0 || actualRow >= dataset.Rows) return double.NaN;
    431             return ((IList<double>)currentInstr.iArg0)[actualRow] * laggedVariableTreeNode.Weight;
    432           }
    433         case OpCodes.Constant: {
    434             var constTreeNode = (ConstantTreeNode)currentInstr.dynamicNode;
    435             return constTreeNode.Value;
    436           }
    437 
    438         //mkommend: this symbol uses the logistic function f(x) = 1 / (1 + e^(-alpha * x) )
    439         //to determine the relative amounts of the true and false branch see http://en.wikipedia.org/wiki/Logistic_function
    440         case OpCodes.VariableCondition: {
    441             if (row < 0 || row >= dataset.Rows) return double.NaN;
    442             var variableConditionTreeNode = (VariableConditionTreeNode)currentInstr.dynamicNode;
    443             double variableValue = ((IList<double>)currentInstr.iArg0)[row];
    444             double x = variableValue - variableConditionTreeNode.Threshold;
    445             double p = 1 / (1 + Math.Exp(-variableConditionTreeNode.Slope * x));
    446 
    447             double trueBranch = Evaluate(dataset, ref row, state);
    448             double falseBranch = Evaluate(dataset, ref row, state);
    449 
    450             return trueBranch * p + falseBranch * (1 - p);
    451           }
    452         default: throw new NotSupportedException();
    453       }
     451
     452      return code[code.Length - 1].value;
    454453    }
    455454  }
Note: See TracChangeset for help on using the changeset viewer.