Index: /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs	(revision 15364)
+++ /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs	(revision 15364)
@@ -0,0 +1,202 @@
+﻿#region License Information
+/* HeuristicLab
+ * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
+ *
+ * This file is part of HeuristicLab.
+ *
+ * HeuristicLab is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * HeuristicLab is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
+ */
+#endregion
+
+using System;
+using System.Collections.Concurrent;
+using System.Collections.Generic;
+using System.Linq;
+using System.Threading;
+using System.Threading.Tasks;
+using HeuristicLab.Analysis;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Data;
+using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
+using HeuristicLab.Optimization;
+using HeuristicLab.Parameters;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+using HeuristicLab.Problems.DataAnalysis;
+using HeuristicLab.Problems.DataAnalysis.Symbolic;
+using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;
+
+namespace HeuristicLab.Algorithms.DataAnalysis.Experimental {
+  // UNFINISHED
+  [Item("Generalized Additive Modelling", "GAM")]
+  [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 102)]
+  [StorableClass]
+  public sealed class GAM : FixedDataAnalysisAlgorithm<IRegressionProblem> {
+    [StorableConstructor]
+    private GAM(bool deserializing) : base(deserializing) { }
+    [StorableHook(HookType.AfterDeserialization)]
+    private void AfterDeserialization() {
+    }
+
+    private GAM(GAM original, Cloner cloner)
+      : base(original, cloner) {
+    }
+    public override IDeepCloneable Clone(Cloner cloner) {
+      return new GAM(this, cloner);
+    }
+
+    public GAM()
+      : base() {
+      Problem = new RegressionProblem();
+      Parameters.Add(new ValueParameter<DoubleValue>("Lambda", "Regularization for smoothing splines", new DoubleValue(1.0)));
+      Parameters.Add(new ValueParameter<IntValue>("Max iterations", "", new IntValue(100)));
+    }
+
+
+    protected override void Run(CancellationToken cancellationToken) {
+
+      double lambda = ((IValueParameter<DoubleValue>)Parameters["Lambda"]).Value.Value;
+      int maxIters = ((IValueParameter<IntValue>)Parameters["Max iterations"]).Value.Value;
+
+      // calculates a GAM model using a linear representation + independent non-linear functions of each variable
+      // using backfitting algorithm (see The Elements of Statistical Learning page 298)
+
+      var problemData = Problem.ProblemData;
+      var y = problemData.TargetVariableTrainingValues.ToArray();
+      var avgY = y.Average();
+      var inputVars = Problem.ProblemData.AllowedInputVariables.ToArray();
+      IRegressionModel[] f = new IRegressionModel[inputVars.Length * 2]; // linear(x) + nonlinear(x)
+      for(int i=0;i<f.Length;i++) {
+        f[i] = new ConstantModel(0.0, problemData.TargetVariable);
+      }
+
+      var rmseTable = new DataTable("RMSE");
+      var rmseRow = new DataRow("RMSE (train)");
+      var rmseRowTest = new DataRow("RMSE (test)");
+      rmseTable.Rows.Add(rmseRow);
+      rmseTable.Rows.Add(rmseRowTest);
+
+      Results.Add(new Result("RMSE", rmseTable));
+      rmseRow.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TrainingIndices).StandardDeviation()); // -1 index to use all predictors
+      rmseRowTest.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TestIndices).StandardDeviation()); 
+
+      // until convergence
+      int iters = 0;
+      var t = new double[y.Length];
+      while (iters++ < maxIters) {
+        int j = 0;
+        foreach (var inputVar in inputVars) {
+          var res = CalculateResiduals(problemData, f, j, avgY, problemData.TrainingIndices);
+          f[j] = RegressLR(problemData, inputVar, res);
+          j++;
+        }
+        foreach (var inputVar in inputVars) {
+          var res = CalculateResiduals(problemData, f, j, avgY, problemData.TrainingIndices);
+          f[j] = RegressSpline(problemData, inputVar, res, lambda);
+          j++;
+        }
+
+        rmseRow.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TrainingIndices).StandardDeviation()); // -1 index to use all predictors
+        rmseRowTest.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TestIndices).StandardDeviation()); 
+      }
+
+      var model = new RegressionEnsembleModel(f.Concat(new[] { new ConstantModel(avgY, problemData.TargetVariable) }));
+      model.AverageModelEstimates = false;
+      Results.Add(new Result("Ensemble solution", model.CreateRegressionSolution((IRegressionProblemData)problemData.Clone())));
+    }
+
+    private double[] CalculateResiduals(IRegressionProblemData problemData, IRegressionModel[] f, int j, double avgY, IEnumerable<int> rows) {
+      var y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, rows);
+      double[] t = y.Select(yi => yi - avgY).ToArray();
+      // collect other predictions
+      for (int k = 0; k < f.Length; k++) {
+        if (k != j) {
+          var pred = f[k].GetEstimatedValues(problemData.Dataset, rows).ToArray();
+          // determine target for this smoother
+          for (int i = 0; i < t.Length; i++) {
+            t[i] -= pred[i];
+          }
+        }
+      }
+      return t;
+    }
+
+    private IRegressionModel RegressLR(IRegressionProblemData problemData, string inputVar, double[] target) {
+      // Umständlich!
+      var ds = new Dataset(new string[] { inputVar, "y" }, new[] { problemData.Dataset.GetDoubleValues(inputVar, problemData.TrainingIndices).ToList(), target.ToList() });
+      var pd = new RegressionProblemData(ds, new string[] { inputVar }, "y");
+      double rmsError, cvRmsError;
+      return LinearRegression.CreateLinearRegressionSolution(pd, out rmsError, out cvRmsError).Model;
+    }
+
+    private IRegressionModel RegressSpline(IRegressionProblemData problemData, string inputVar, double[] target, double lambda) {
+      // Umständlich!
+      return Splines.CalculateSmoothingSplineReinsch(
+        problemData.Dataset.GetDoubleValues(inputVar, problemData.TrainingIndices).ToArray(),
+        (double[])target.Clone(),
+        new string[] { inputVar },
+        s: lambda,
+        targetVar: problemData.TargetVariable);
+    }
+
+  }
+
+
+  // UNFINISHED
+  public class RBFModel : NamedItem, IRegressionModel {
+    private alglib.rbfmodel model;
+
+    public string TargetVariable { get; set; }
+
+    public IEnumerable<string> VariablesUsedForPrediction { get; private set; }
+    private ITransformation<double>[] scaling;
+
+    public event EventHandler TargetVariableChanged;
+
+    public RBFModel(RBFModel orig, Cloner cloner) : base(orig, cloner) {
+      this.TargetVariable = orig.TargetVariable;
+      this.VariablesUsedForPrediction = orig.VariablesUsedForPrediction.ToArray();
+      this.model = (alglib.rbfmodel)orig.model.make_copy();
+      this.scaling = orig.scaling.Select(s => cloner.Clone(s)).ToArray();
+    }
+    public RBFModel(alglib.rbfmodel model, string targetVar, string[] inputs, IEnumerable<ITransformation<double>> scaling) : base("RBFModel", "RBFModel") {
+      this.model = model;
+      this.TargetVariable = targetVar;
+      this.VariablesUsedForPrediction = inputs;
+      this.scaling = scaling.ToArray();
+    }
+
+    public override IDeepCloneable Clone(Cloner cloner) {
+      return new RBFModel(this, cloner);
+    }
+
+    public IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
+      return new RegressionSolution(this, (IRegressionProblemData)problemData.Clone());
+    }
+
+    public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
+      double[] x = new double[VariablesUsedForPrediction.Count()];
+      double[] y;
+      foreach (var r in rows) {
+        int c = 0;
+        foreach (var v in VariablesUsedForPrediction) {
+          x[c] = scaling[c].Apply(dataset.GetDoubleValue(v, r).ToEnumerable()).First(); // OUCH!
+          c++;
+        }
+        alglib.rbfcalc(model, x, out y);
+        yield return y[0];
+      }
+    }
+  }
+}
Index: /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/HeuristicLab.Algorithms.DataAnalysis.Experimental.csproj
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/HeuristicLab.Algorithms.DataAnalysis.Experimental.csproj	(revision 15363)
+++ /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/HeuristicLab.Algorithms.DataAnalysis.Experimental.csproj	(revision 15364)
@@ -134,4 +134,5 @@
   <ItemGroup>
     <Compile Include="ForwardSelection.cs" />
+    <Compile Include="GAM.cs" />
     <Compile Include="RBF.cs" />
     <Compile Include="Splines.cs" />
Index: /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs	(revision 15363)
+++ /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs	(revision 15364)
@@ -37,4 +37,6 @@
 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,7 +68,10 @@
           "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,5 +94,5 @@
         switch (type) {
           case "Monotone":
-            alglib.spline1dbuildmonotone(x, y, out c); 
+            alglib.spline1dbuildmonotone(x, y, out c);
             AddAlglibSplineResult(c, inputVars);
             break;
@@ -132,8 +137,243 @@
               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,8 +386,123 @@
       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
