source: branches/HeuristicLab.Problems.GeneticProgramming.BloodGlucosePrediction/Interpreter.cs @ 14310

using System;
using System.Collections;
using System.Collections.Generic;
using System.Collections.ObjectModel;
using System.Collections.Specialized;
using System.Drawing.Design;
using System.Linq;
using HeuristicLab.Common;
using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
using HeuristicLab.Problems.DataAnalysis;
using HeuristicLab.Problems.DataAnalysis.Symbolic;

namespace HeuristicLab.Problems.GeneticProgramming.GlucosePrediction {
  public static class Interpreter {
    private class Data {
      public double[] realGluc;
      public double[] realIns;
      public double[] realCh;
      public Dictionary<ISymbolicExpressionTreeNode, double[]> precalculatedValues;
    }

    public static IEnumerable<double> Apply(ISymbolicExpressionTreeNode model, IDataset dataset, IEnumerable<int> rows) {
      double[] targetGluc = dataset.GetDoubleValues("Glucose_target", rows).ToArray(); // only for skipping rows for which we should not produce an output

      var data = new Data {
        realGluc = dataset.GetDoubleValues("Glucose_Interpol", rows).ToArray(),
        realIns = dataset.GetDoubleValues("Insuline", rows).ToArray(),
        realCh = dataset.GetDoubleValues("CH", rows).ToArray(),
        precalculatedValues = CreatePrecalculatedValues(model, dataset)
      };
      var predictions = new double[targetGluc.Length];
      var rowsEnumerator = rows.GetEnumerator();
      for (int k = 0; k < predictions.Length; k++, rowsEnumerator.MoveNext()) {
        if (double.IsNaN(targetGluc[k])) {
          predictions[k] = double.NaN;
        } else {
          var rawPred = InterpretRec(model, data, rowsEnumerator.Current);
          predictions[k] = rawPred; 
        }
      }
      return predictions;
    }

    private static Dictionary<ISymbolicExpressionTreeNode, double[]> CreatePrecalculatedValues(ISymbolicExpressionTreeNode root, IDataset dataset) {
      var dict = new Dictionary<ISymbolicExpressionTreeNode, double[]>();
      // here we integrate ins or ch inputs over the whole day to generate smoothed ins/ch values with the same number of rows
      // the integrated values are reset to zero whenever a new evluation period starts
      foreach (var node in root.IterateNodesPrefix()) {
        var curvedInsNode = node as CurvedInsVariableTreeNode;
        var curvedChNode = node as CurvedChVariableTreeNode;
        if (curvedInsNode != null) {
          dict.Add(curvedInsNode, Integrate(curvedInsNode, dataset));
        } else if (curvedChNode != null) {
          dict.Add(curvedChNode, Integrate(curvedChNode, dataset));
        }
      }
      return dict;
    }

    private static double[] Integrate(CurvedInsVariableTreeNode node, IDataset dataset) {
      // d Q1 / dt = ins(t) - alpha * Q1(t)
      // d Q2 / dt = alpha * (Q1(t) - Q2(t))
      // S = Q1(t) + Q2(t)
      var alpha = node.Alpha;

      var ins = dataset.GetReadOnlyDoubleValues("Insuline");
      var time = dataset.GetReadOnlyDoubleValues("HourMin").ToArray();

      // TODO reset for new time intervals

      double q1, q2, q1_prev, q2_prev;
      // starting values: zeros
      q1 = q2 = q1_prev = q2_prev = 0;
      double[] s = new double[dataset.Rows];

      for (int t = 1; t < dataset.Rows; t++) {
        if (IsStartOfNewPeriod(time, t)) {
          q1 = q2 = q1_prev = q2_prev = 0;
        }
        q1 = q1_prev + ins[t] - alpha * q1_prev;
        q2 = q2_prev + alpha * (q1_prev - q2_prev);
        s[t] = q1 + q2;
        q1_prev = q1;
        q2_prev = q2;

      }
      return s;
    }

    private static bool IsStartOfNewPeriod(double[] time, int t) {
      return t == 0 ||
             (time[t].IsAlmost(2005) && !time[t - 1].IsAlmost(2000));
    }


    private static double[] Integrate(CurvedChVariableTreeNode node, IDataset dataset) {
      // d Q1 / dt = ins(t) - alpha * Q1(t)
      // d Q2 / dt = alpha * (Q1(t) - Q2(t))
      // S = Q1(t) + Q2(t)
      var alpha = node.Alpha;

      var ch = dataset.GetReadOnlyDoubleValues("CH");
      var time = dataset.GetReadOnlyDoubleValues("HourMin").ToArray();

      // TODO reset for new time intervals

      double q1, q2, q1_prev, q2_prev;
      // starting values: zeros
      q1 = q2 = q1_prev = q2_prev = 0;
      double[] s = new double[dataset.Rows];

      for (int t = 1; t < dataset.Rows; t++) {
        if (IsStartOfNewPeriod(time, t)) {
          q1 = q2 = q1_prev = q2_prev = 0;
        }
        q1 = q1_prev + ch[t] - alpha * q1_prev;
        q2 = q2_prev + alpha * (q1_prev - q2_prev);
        s[t] = q1 + q2;
        q1_prev = q1;
        q2_prev = q2;

      }
      return s;
    }

    private static double InterpretRec(ISymbolicExpressionTreeNode node, Data data, int k) {
      if (node.Symbol is SimpleSymbol) {
        switch (node.Symbol.Name) {
          case "+":
          case "+Ins":
          case "+Ch": {
              return InterpretRec(node.GetSubtree(0), data, k) + InterpretRec(node.GetSubtree(1), data, k);
            }
          case "-":
          case "-Ins":
          case "-Ch": {
              return InterpretRec(node.GetSubtree(0), data, k) - InterpretRec(node.GetSubtree(1), data, k);
            }
          case "*":
          case "*Ins":
          case "*Ch": {
              return InterpretRec(node.GetSubtree(0), data, k) * InterpretRec(node.GetSubtree(1), data, k);
            }
          case "/Ch":
          case "/Ins":
          case "/": {
              return InterpretRec(node.GetSubtree(0), data, k) / InterpretRec(node.GetSubtree(1), data, k);
            }
          case "Exp":
          case "ExpIns":
          case "ExpCh": {
              return Math.Exp(InterpretRec(node.GetSubtree(0), data, k));
            }
          case "Sin":
          case "SinIns":
          case "SinCh": {
              return Math.Sin(InterpretRec(node.GetSubtree(0), data, k));
            }
          case "CosCh":
          case "CosIns":
          case "Cos": {
              return Math.Cos(InterpretRec(node.GetSubtree(0), data, k));
            }
          case "LogCh":
          case "LogIns":
          case "Log": {
              return Math.Log(InterpretRec(node.GetSubtree(0), data, k));
            }
          case "Func": {
              // <exprgluc> + <exprch> - <exprins>
              return InterpretRec(node.GetSubtree(0), data, k)
                     + InterpretRec(node.GetSubtree(1), data, k)
                     - InterpretRec(node.GetSubtree(2), data, k);
            }
          case "ExprGluc": {
              return InterpretRec(node.GetSubtree(0), data, k);
            }
          case "ExprCh": {
              return InterpretRec(node.GetSubtree(0), data, k);
            }
          case "ExprIns": {
              return InterpretRec(node.GetSubtree(0), data, k);
            }
          default: {
              throw new InvalidProgramException("Found an unknown symbol " + node.Symbol);
            }
        }
      } else if (node.Symbol is PredictedGlucoseVariableSymbol) {
        throw new NotSupportedException();
      } else if (node.Symbol is RealGlucoseVariableSymbol) {
        var n = (RealGlucoseVariableTreeNode)node;
        if (k + n.RowOffset < 0 || k + n.RowOffset >= data.realGluc.Length) return double.NaN;
        return data.realGluc[k + n.RowOffset];
      } else if (node.Symbol is CurvedChVariableSymbol) {
        return data.precalculatedValues[node][k];
      } else if (node.Symbol is RealInsulineVariableSymbol) {
        throw new NotSupportedException();
      } else if (node.Symbol is CurvedInsVariableSymbol) {
        return data.precalculatedValues[node][k];
      } else if (node.Symbol is Constant) {
        var n = (ConstantTreeNode)node;
        return n.Value;
      } else {
        throw new InvalidProgramException("found unknown symbol " + node.Symbol);
      }
    }

    private static double Beta(double x, double alpha, double beta) {
      return 1.0 / alglib.beta(alpha, beta) * Math.Pow(x, alpha - 1) * Math.Pow(1 - x, beta - 1);
    }
  }
}