Changeset 15364

Timestamp:

09/15/17 18:20:14 (7 years ago)

Author:

gkronber

Message:

#2789 experimenting with GAMs and Splines

Location:

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental

Files:

• GAM.cs (added)
• HeuristicLab.Algorithms.DataAnalysis.Experimental.csproj (modified) (1 diff)
• Splines.cs (modified) (5 diffs)

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/HeuristicLab.Algorithms.DataAnalysis.Experimental.csproj

@@ -134 +134 @@
   <ItemGroup>
     <Compile Include="ForwardSelection.cs" />
+    <Compile Include="GAM.cs" />
     <Compile Include="RBF.cs" />
     <Compile Include="Splines.cs" />

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs

@@ -37 +37 @@
 using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;
 using MathNet.Numerics.Interpolation;
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.LinearAlgebra.Double.Solvers;
 
 namespace HeuristicLab.Algorithms.DataAnalysis.Experimental {
@@ -66 +68 @@
           "LogLinear (Math.NET)",
           "Common (Math.NET)",
+          "Smoothing Spline (Mine)",
+          "Smoothing Spline (Reinsch)",
       }.Select(s => new StringValue(s)));
 
       Parameters.Add(new ConstrainedValueParameter<StringValue>("Type", "The type of spline (as supported by alglib)", validTypes, validTypes.First()));
+      Parameters.Add(new ValueParameter<DoubleValue>("Lambda", "Regularization parameter for smoothing splines (0..+inf)", new DoubleValue(100)));
       Problem = new RegressionProblem();
     }
@@ -89 +94 @@
         switch (type) {
           case "Monotone":
-            alglib.spline1dbuildmonotone(x, y, out c); 
+            alglib.spline1dbuildmonotone(x, y, out c);
             AddAlglibSplineResult(c, inputVars);
             break;
@@ -132 +137 @@
               break;
             }
+          case "Smoothing Spline (Mine)": {
+              CalculateSmoothingSpline(x, y, inputVars);
+              break;
+            }
+          case "Smoothing Spline (Reinsch)": {
+              CalculateSmoothingSplineReinsch(x, y, inputVars);
+              break;
+            }
           default: throw new NotSupportedException();
         }
 
       }
+    }
+
+    // array with non-zero lower index
+    private class MyArray<T> {
+      private T[] arr;
+      private int lowerBound;
+
+      public T this[int key] {
+        get {
+          return arr[key - lowerBound];
+        }
+        set {
+          arr[key - lowerBound] = value;
+        }
+      }
+
+      public MyArray(int lowerBound, int numElements) {
+        this.lowerBound = lowerBound;
+        arr = new T[numElements];
+      }
+      public MyArray(int lowerBound, T[] source) : this(lowerBound, source.Length) {
+        Array.Copy(source, arr, source.Length);
+      }
+
+      public T[] ToArray() {
+        var res = new T[arr.Length];
+        Array.Copy(arr, res, res.Length);
+        return res;
+      }
+    }
+
+    public void CalculateSmoothingSplineReinsch(double[] xOrig, double[] yOrig, string[] inputVars) {
+      double s = ((IValueParameter<DoubleValue>)Parameters["Lambda"]).Value.Value; // controls extent of smoothing
+
+      var model = CalculateSmoothingSplineReinsch(xOrig, yOrig, inputVars, s, Problem.ProblemData.TargetVariable);
+
+      Results.Add(new Result("Solution", new RegressionSolution(model, (IRegressionProblemData)Problem.ProblemData.Clone())));
+
+
+    }
+
+    public static IRegressionModel CalculateSmoothingSplineReinsch(double[] xOrig, double[] yOrig, string[] inputVars, double s, string targetVar) {
+      var minX = xOrig.Min();
+      var maxX = xOrig.Max();
+      var range = maxX - minX;
+
+      SortAndBin(xOrig, yOrig, out xOrig, out yOrig, scaling: true);
+
+      // See Smoothing by Spline Functions, Reinsch, 1967
+      // move x and y into an array indexed with 1..n to match indexing in Reinsch paper
+      int n = xOrig.Length;
+      var x = new MyArray<double>(1, xOrig);
+      var y = new MyArray<double>(1, yOrig);
+      var inv_dy = new MyArray<double>(1, Enumerable.Repeat(1.0, n).ToArray()); // equal weights
+
+      int n1 = 1;
+      int n2 = n - 1;
+
+      // results
+      var a = new MyArray<double>(n1, n);
+      var b = new MyArray<double>(n1, n);
+      var c = new MyArray<double>(n1, n);
+      var d = new MyArray<double>(n1, n);
+
+      // smooth
+      {
+        int i, m1, m2; double e, f, f2, g = 0.0, h, p;
+        MyArray<double>
+          r = new MyArray<double>(n1 - 1, n + 2),
+          r1 = new MyArray<double>(n1 - 1, n + 2),
+          r2 = new MyArray<double>(n1 - 1, n + 2),
+          t = new MyArray<double>(n1 - 1, n + 2),
+          t1 = new MyArray<double>(n1 - 1, n + 2),
+          u = new MyArray<double>(n1 - 1, n + 2),
+          v = new MyArray<double>(n1 - 1, n + 2);
+        m1 = n1 - 1;
+        m2 = n2 + 1;
+        r[m1] = r[n1] = r1[n2] = r2[n2] = r2[m2] =
+          u[m1] = u[n1] = u[n2] = u[m2] = p = 0.0;
+        m1 = n1 + 1;
+        m2 = n2 - 1;
+        h = x[m1] - x[n1]; f = (y[m1] - y[n1]) / h;
+        for (i = m1; i <= m2; i++) {
+          g = h; h = x[i + 1] - x[i];
+          e = f; f = (y[i + 1] - y[i]) / h;
+          a[i] = f - e; t[i] = 2 * (g + h) / 3; t1[i] = h / 3;
+          r2[i] = inv_dy[i - 1] / g; r[i] = inv_dy[i + 1] / h;
+          r1[i] = -inv_dy[i] / g - inv_dy[i] / h;
+        }
+        for (i = m1; i <= m2; i++) {
+          b[i] = r[i] * r[i] + r1[i] * r1[i] + r2[i] * r2[i];
+          c[i] = r[i] * r1[i + 1] + r1[i] * r2[i + 1];
+          d[i] = r[i] * r2[i + 2];
+        }
+        f2 = -s;
+        next:
+        for (i = m1; i <= m2; i++) {
+          r1[i - 1] = f * r[i - 1]; r2[i - 2] = g * r[i - 2];
+          r[i] = 1 / (p * b[i] + t[i] - f * r1[i - 1] - g * r2[i - 2]);
+          u[i] = a[i] - r1[i - 1] * u[i - 1] - r2[i - 2] * u[i - 2];
+          f = p * c[i] + t1[i] - h * r1[i - 1]; g = h; h = d[i] * p;
+        }
+        for (i = m2; i >= m1; i--) {
+          u[i] = r[i] * u[i] - r1[i] * u[i + 1] - r2[i] * u[i + 2];
+        }
+        e = h = 0;
+        for (i = n1; i <= m2; i++) {
+          g = h; h = (u[i + 1] - u[i]) / (x[i + 1] - x[i]);
+          v[i] = (h - g) * inv_dy[i] * inv_dy[i]; e = e + v[i] * (h - g);
+        }
+        g = v[n2] = -h * inv_dy[n2] * inv_dy[n2]; e = e - g * h;
+        g = f2; f2 = e * p * p;
+        if (f2 >= s || f2 <= g) goto fin;
+        f = 0; h = (v[m1] - v[n1]) / (x[m1] - x[n1]);
+        for (i = m1; i <= m2; i++) {
+          g = h; h = (v[i + 1] - v[i]) / (x[i + 1] - x[i]);
+          g = h - g - r1[i - 1] * r[i - 1] - r2[i - 2] * r[i - 2];
+          f = f + g * r[i] * g; r[i] = g;
+        }
+        h = e - p * f; if (h <= 0) goto fin;
+        p = p + (s - f2) / ((Math.Sqrt(s / e) + p) * h); goto next;
+
+        fin:
+        for (i = n1; i <= n2; i++) {
+          a[i] = y[i] - p * v[i];
+          c[i] = u[i];
+        }
+        for (i = n1; i <= m2; i++) {
+          h = x[i + 1] - x[i];
+          d[i] = (c[i + 1] - c[i]) / (3 * h);
+          b[i] = (a[i + 1] - a[i]) / h - (h * d[i] + c[i]) * h;
+        }
+      }
+
+      return new ReinschSmoothingSplineModel(a.ToArray(), b.ToArray(), c.ToArray(), d.ToArray(), xOrig, targetVar, inputVars, minX, 1.0 / range);
+    }
+
+    private void CalculateSmoothingSpline(double[] x, double[] y, string[] inputVars) {
+      // see Smoothing and Non-Parametric Regression, Germán Rodríguez, 2001 2.3.1
+      SortAndBin(x, y, out x, out y, scaling: false);
+      int n = x.Length;
+
+      SparseMatrix delta = new SparseMatrix(n - 2, n);
+      // double[,] delta = new double[n - 2, n];
+      //double[,] W = new double[n - 2, n - 2];
+      SparseMatrix W = new SparseMatrix(n - 2, n - 2);
+      Matrix WInvD = new DenseMatrix(n - 2, n);
+
+      // double[,] W_inv_D = new double[n - 2, n];
+      // double[,] K = new double[n, n];
+
+      // go over successive knots to determine distances and fill Delta and W
+      for (int i = 0; i < n - 2; i++) {
+        double h = x[i + 1] - x[i];
+        double h_next = x[i + 2] - x[i + 1];
+        delta[i, i] = 1.0 / h;
+        delta[i, i + 1] = -1.0 / h - 1.0 / h_next;
+        delta[i, i + 2] = 1.0 / h_next;
+        W[i, i] = (h + h_next) / 3;
+        if (i > 0) {
+          W[i - 1, i] =
+           W[i, i - 1] = h / 6.0;
+        }
+      }
+
+      var solvResult = W.TrySolveIterative(delta, WInvD, new MlkBiCgStab());
+
+      // alglib.ablas.rmatrixgemm(n - 2, n, n - 2, 1.0, W, 0, 0, 0, delta, 0, 0, 0, 1.0, W_inv_D, 0, 0); // W^-1(M = n-2, K = n-2) D(K = n-2, N=n)   
+      // alglib.ablas.rmatrixgemm(n, n, n - 2, 1.0, delta, 0, 0, 1, W_inv_D, 0, 0, 0, 1.0, K, 0, 0); // D(M=n-2, K=n)^T * W^-1D (K=n, N=n-2)
+
+      var K = delta.TransposeThisAndMultiply(WInvD);
+
+      double lambda = ((IValueParameter<DoubleValue>)Parameters["Lambda"]).Value.Value;
+
+      for (int i = 0; i < n; i++) {
+        for (int j = 0; j < n; j++) {
+          K[i, j] *= lambda;
+          if (i == j) K[i, j] += 1;
+        }
+      }
+
+      // solve for y
+      // double[] s;
+      // int solverInfo;
+      // alglib.densesolverreport solverRep;
+      // alglib.rmatrixsolve(K, n, y, out solverInfo, out solverRep, out s);
+
+      var s = K.Solve(new DenseVector(y)).ToArray();
+
+      Results.Add(new Result("Solution", new RegressionSolution(new SmoothingSplineModel(s, x, Problem.ProblemData.TargetVariable, inputVars),
+ (IRegressionProblemData)Problem.ProblemData.Clone())));
+    }
+
+    private static void SortAndBin(double[] x, double[] y, out double[] x2, out double[] y2, bool scaling = false) {
+      // sort y by x
+      Array.Sort(x, y);
+
+      var xl = new List<double>();
+      var yl = new List<double>();
+
+      int n = x.Length;
+      var range = x[n - 1] - x[0];
+      var binSize = range / n / 20;
+      {
+        // binning
+        int i = 0;
+        while (i < n) {
+          int k = 0;
+          int j = i;
+          double sumX = 0.0;
+          double sumY = 0.0;
+          while (j < n && x[j] - x[i] <= binSize) {
+            k++;
+            sumX += x[j];
+            sumY += y[j];
+            j++;
+          }
+          var avgX = sumX / k;
+          if (scaling) avgX = (avgX - x[0]) / range;
+          xl.Add(avgX);     // scale 
+          yl.Add(sumY / k);
+          i += k;
+        }
+      }
+
+      x2 = xl.ToArray();
+      y2 = yl.ToArray();
     }
 
@@ -146 +386 @@
       Results.Add(new Result("Solution", new RegressionSolution(new MathNetSplineModel(c, Problem.ProblemData.TargetVariable, inputVars),
         (IRegressionProblemData)Problem.ProblemData.Clone())));
-
     }
   }
 
+
+
+  // UNFINISHED
+  public class ReinschSmoothingSplineModel : NamedItem, IRegressionModel {
+    private double[] a;
+    private double[] b;
+    private double[] c;
+    private double[] d;
+    private double[] x;
+    private double offset;
+    private double scale;
+
+    public string TargetVariable { get; set; }
+
+    public IEnumerable<string> VariablesUsedForPrediction { get; private set; }
+
+    public event EventHandler TargetVariableChanged;
+
+    public ReinschSmoothingSplineModel(ReinschSmoothingSplineModel orig, Cloner cloner) : base(orig, cloner) {
+      this.TargetVariable = orig.TargetVariable;
+      this.VariablesUsedForPrediction = orig.VariablesUsedForPrediction.ToArray();
+      this.a = orig.a;
+      this.b = orig.b;
+      this.c = orig.c;
+      this.d = orig.d;
+      this.x = orig.x;
+      this.scale = orig.scale;
+      this.offset = orig.offset;
+    }
+    public ReinschSmoothingSplineModel(double[] a, double[] b, double[] c, double[] d, double[] x, string targetVar, string[] inputs, double offset, double scale) : base("SplineModel", "SplineModel") {
+      this.a = a;
+      this.b = b;
+      this.c = c;
+      this.d = d;
+      this.x = x;
+      this.TargetVariable = targetVar;
+      this.VariablesUsedForPrediction = inputs;
+      this.scale = scale;
+      this.offset = offset;
+    }
+
+    public override IDeepCloneable Clone(Cloner cloner) {
+      return new ReinschSmoothingSplineModel(this, cloner);
+    }
+
+    public IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
+      return new RegressionSolution(this, (IRegressionProblemData)problemData.Clone());
+    }
+
+    public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
+      int n = x.Length;
+      foreach (var xx in dataset.GetDoubleValues(VariablesUsedForPrediction.First(), rows).Select(xi => (xi - offset) * scale)) {
+        if (xx < x[0]) yield return a[0];
+        else if (xx >= x[n - 1]) yield return a[n - 2];
+        else {
+          // binary search
+          int lower = 0;
+          int upper = n - 1;
+          while (true) {
+            if (upper < lower) throw new InvalidProgramException();
+            int i = lower + (upper - lower) / 2;
+            if (x[i] <= xx && xx < x[i + 1]) {
+              double h = xx - x[i];
+              yield return a[i] + h * (b[i] + h * (c[i] + h * d[i]));
+              break;
+            } else if (xx < x[i]) {
+              upper = i - 1;
+            } else {
+              lower = i + 1;
+            }
+          }
+        }
+      }
+    }
+  }
+
+  // UNFINISHED
+  public class SmoothingSplineModel : NamedItem, IRegressionModel {
+    private double[] s;
+    private IInterpolation interpolation;
+
+    public string TargetVariable { get; set; }
+
+    public IEnumerable<string> VariablesUsedForPrediction { get; private set; }
+
+    public event EventHandler TargetVariableChanged;
+
+    public SmoothingSplineModel(SmoothingSplineModel orig, Cloner cloner) : base(orig, cloner) {
+      this.TargetVariable = orig.TargetVariable;
+      this.VariablesUsedForPrediction = orig.VariablesUsedForPrediction.ToArray();
+      this.s = orig.s; // TODO
+      this.interpolation = orig.interpolation;
+    }
+    public SmoothingSplineModel(double[] s, double[] x, string targetVar, string[] inputs) : base("SplineModel", "SplineModel") {
+      this.s = s;
+      this.TargetVariable = targetVar;
+      this.VariablesUsedForPrediction = inputs;
+      this.interpolation = MathNet.Numerics.Interpolate.CubicSpline(x, s);
+    }
+
+    public override IDeepCloneable Clone(Cloner cloner) {
+      return new SmoothingSplineModel(this, cloner);
+    }
+
+    public IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
+      return new RegressionSolution(this, (IRegressionProblemData)problemData.Clone());
+    }
+
+    public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
+      foreach (var x in dataset.GetDoubleValues(VariablesUsedForPrediction.First(), rows)) {
+
+        yield return interpolation.Interpolate(x);
+
+      }
+    }
+  }
 
   // UNFINISHED