Changeset 14310 for branches/HeuristicLab.Problems.GeneticProgramming.BloodGlucosePrediction/Interpreter.cs

Timestamp:

09/27/16 16:30:54 (8 years ago)

Author:

gkronber

Message:

changed interpreter and grammar to smooth ch and ins uptake using a compartement model

File:

• branches/HeuristicLab.Problems.GeneticProgramming.BloodGlucosePrediction/Interpreter.cs (modified) (5 diffs)

branches/HeuristicLab.Problems.GeneticProgramming.BloodGlucosePrediction/Interpreter.cs

@@ -2 +2 @@
 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;
@@ -15 +17 @@
       public double[] realIns;
       public double[] realCh;
-      public double[] predGluc;
+      public Dictionary<ISymbolicExpressionTreeNode, double[]> precalculatedValues;
     }
 
@@ -25 +27 @@
         realIns = dataset.GetDoubleValues("Insuline", rows).ToArray(),
         realCh = dataset.GetDoubleValues("CH", rows).ToArray(),
+        precalculatedValues = CreatePrecalculatedValues(model, dataset)
       };
-      data.predGluc = new double[data.realGluc.Length];
-      Array.Copy(data.realGluc, data.predGluc, data.predGluc.Length);
-      for (int k = 0; k < data.predGluc.Length; k++) {
+      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])) {
-          data.predGluc[k] = double.NaN;
+          predictions[k] = double.NaN;
         } else {
-          var rawPred = InterpretRec(model, data, k);
-          data.predGluc[k] = Math.Max(0, Math.Min(400, rawPred)); // limit output values of the model to 0 ... 400
-        }
-      }
-      return data.predGluc;
+          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;
     }
 
@@ -102 +187 @@
         }
       } else if (node.Symbol is PredictedGlucoseVariableSymbol) {
-        var n = (PredictedGlucoseVariableTreeNode)node;
-        return n.Weight * data.predGluc[k + n.RowOffset];
+        throw new NotSupportedException();
       } else if (node.Symbol is RealGlucoseVariableSymbol) {
         var n = (RealGlucoseVariableTreeNode)node;
-        return n.Weight * data.realGluc[k + n.RowOffset];
+        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) {
-        var n = (CurvedChVariableTreeNode)node;
-        double prevVal;
-        int prevValDistance;
-        GetPrevDataAndDistance(data.realCh, k, out prevVal, out prevValDistance, maxDistance: 48);
-        return n.Weight * prevVal * Beta(prevValDistance / 48.0, n.Alpha, n.Beta);
+        return data.precalculatedValues[node][k];
       } else if (node.Symbol is RealInsulineVariableSymbol) {
-        var n = (RealInsulineVariableTreeNode)node;
-        return n.Weight * data.realIns[k + n.RowOffset];
+        throw new NotSupportedException();
       } else if (node.Symbol is CurvedInsVariableSymbol) {
-        var n = (CurvedInsVariableTreeNode)node;
-        double maxVal;
-        int maxValDistance;
-        var sum = GetSumOfValues(48, k, data.realIns);
-
-        GetMaxValueAndDistance(data.realIns, k, out maxVal, out maxValDistance, maxDistance: 48);
-        return n.Weight * (sum - maxVal) * maxVal * Beta(maxValDistance / 48.0, n.Alpha, n.Beta);
+        return data.precalculatedValues[node][k];
       } else if (node.Symbol is Constant) {
         var n = (ConstantTreeNode)node;
@@ -135 +209 @@
       return 1.0 / alglib.beta(alpha, beta) * Math.Pow(x, alpha - 1) * Math.Pow(1 - x, beta - 1);
     }
-
-    private static void GetPrevDataAndDistance(double[] vals, int k, out double val, out int dist, int maxDistance = 48, double threshold = 0.0) {
-      // look backward from the current idx k and find the first value above the threshold
-      for (int i = k; i >= 0 && i >= (k - maxDistance); i--) {
-        if (vals[i] > threshold) {
-          val = vals[i];
-          dist = k - i;
-          return;
-        }
-      }
-      val = 0;
-      dist = maxDistance;
-    }
-
-    private static double GetSumOfValues(int windowSize, int k, double[] vals) {
-      var sum = 0.0;
-      for (int i = k; i >= 0 && i >= k - windowSize; i--)
-        sum += vals[i];
-      return sum;
-    }
-
-    private static void GetMaxValueAndDistance(double[] vals, int k, out double maxVal, out int dist, int maxDistance = 48) {
-      // look backward from the current idx k and find the max value and it's distance
-      maxVal = vals[k];
-      dist = 0;
-      for (int i = k; i >= 0 && i >= (k - maxDistance); i--) {
-        if (vals[i] > maxVal) {
-          maxVal = vals[i];
-          dist = k - i;
-        }
-      }
-    }
   }
 }