Index: /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/CubicSplineGCV.cs
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/CubicSplineGCV.cs	(revision 15467)
+++ /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/CubicSplineGCV.cs	(revision 15468)
@@ -178,5 +178,5 @@
     }
     public static ReinschSmoothingSplineModel CalculateCubicSpline(double[] x, double[] y,
-      string targetVariable, string[] inputVars,
+      string targetVariable, string[] inputVars, 
       out CubGcvReport report) {
       var w = Enumerable.Repeat(1.0, x.Length).ToArray();
Index: /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs	(revision 15467)
+++ /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs	(revision 15468)
@@ -21,20 +21,15 @@
 
 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 {
@@ -228,10 +223,10 @@
           product = product.Zip(problemData.Dataset.GetDoubleValues(inputVars[i], problemData.TrainingIndices), (pi, vi) => pi * vi).ToArray();
         }
-        // CubicSplineGCV.CubGcvReport report;
-        // return CubicSplineGCV.CalculateCubicSpline(
-        //   product,
-        //   (double[])target.Clone(), 
-        //   problemData.TargetVariable, inputVars, out report
-        //   );
+        CubicSplineGCV.CubGcvReport report;
+        return CubicSplineGCV.CalculateCubicSpline(
+          product,
+          (double[])target.Clone(), 
+          problemData.TargetVariable, inputVars, out report
+          );
 
         double optTolerance; double cvRMSE;
@@ -274,5 +269,5 @@
         // return Splines.CalculatePenalizedRegressionSpline(product, (double[])target.Clone(), lambda, problemData.TargetVariable, inputVars, out avgTrainRMSE, out cvRMSE, out bestPredictions);
         SBART.SBART_Report rep;
-        var model = SBART.CalculateSBART(product, (double[])target.Clone(), problemData.TargetVariable, inputVars, (float)lambda, out rep);
+        var model = SBART.CalculateSBART(product, (double[])target.Clone(), problemData.TargetVariable, inputVars, out rep);
         Console.WriteLine("{0} {1:N5} {2:N5} {3:N5} {4:N5}", string.Join(",", inputVars), rep.gcv, rep.leverage.Sum(), product.StandardDeviation(), target.StandardDeviation());
         return model;
Index: /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/SBART.cs
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/SBART.cs	(revision 15467)
+++ /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/SBART.cs	(revision 15468)
@@ -112,4 +112,60 @@
     [DllImport("sbart_x64.dll", CallingConvention = CallingConvention.Cdecl, EntryPoint = "setreg")]
     public static extern void setreg_x64(float[] x, float[] y, float[] w, ref int n, float[] xw, ref int nx, ref float min, ref float range, float[] knot, ref int nk);
+    /*
+     * calculates value at  x  of  jderiv-th derivative of spline from b-repr.
+    c  the spline is taken to be continuous from the right, EXCEPT at the
+    c  rightmost knot, where it is taken to be continuous from the left.
+    c
+    c******  i n p u t ******
+    c  t, bcoef, n, k......forms the b-representation of the spline  f  to
+    c        be evaluated. specifically,
+    c  t.....knot sequence, of length  n+k, assumed nondecreasing.
+    c  bcoef.....b-coefficient sequence, of length  n .
+    c  n.....length of  bcoef  and dimension of spline(k,t),
+    c        a s s u m e d  positive .
+    c  k.....order of the spline .
+    c
+    c  w a r n i n g . . .   the restriction  k .le. kmax (=20)  is imposed
+    c        arbitrarily by the dimension statement for  aj, dl, dr  below,
+    c        but is  n o w h e r e  c h e c k e d  for.
+    c
+    c  x.....the point at which to evaluate .
+    c  jderiv.....integer giving the order of the derivative to be evaluated
+    c        a s s u m e d  to be zero or positive.
+    c
+    c******  o u t p u t  ******
+    c  bvalue.....the value of the (jderiv)-th derivative of  f  at  x .
+    c
+    c******  m e t h o d  ******
+    c     The nontrivial knot interval  (t(i),t(i+1))  containing  x  is lo-
+    c  cated with the aid of  interv . The  k  b-coeffs of  f  relevant for
+    c  this interval are then obtained from  bcoef (or taken to be zero if
+    c  not explicitly available) and are then differenced  jderiv  times to
+    c  obtain the b-coeffs of  (d**jderiv)f  relevant for that interval.
+    c  Precisely, with  j = jderiv, we have from x.(12) of the text that
+    c
+    c     (d**j)f  =  sum ( bcoef(.,j)*b(.,k-j,t) )
+    c
+    c  where
+    c                   / bcoef(.),                     ,  j .eq. 0
+    c                   /
+    c    bcoef(.,j)  =  / bcoef(.,j-1) - bcoef(.-1,j-1)
+    c                   / ----------------------------- ,  j .gt. 0
+    c                   /    (t(.+k-j) - t(.))/(k-j)
+    c
+    c     Then, we use repeatedly the fact that
+    c
+    c    sum ( a(.)*b(.,m,t)(x) )  =  sum ( a(.,x)*b(.,m-1,t)(x) )
+    c  with
+    c                 (x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1)
+    c    a(.,x)  =    ---------------------------------------
+    c                 (x - t(.))      + (t(.+m-1) - x)
+    c
+    c  to write  (d**j)f(x)  eventually as a linear combination of b-splines
+    c  of order  1 , and the coefficient for  b(i,1,t)(x)  must then be the
+    c  desired number  (d**j)f(x). (see x.(17)-(19) of text).
+     */
+    [DllImport("sbart_x64.dll", CallingConvention = CallingConvention.Cdecl, EntryPoint = "bvalue")]
+    public static extern float bvalue(float[] t, float[] bcoeff, ref int n, ref int k, ref float x, ref int jderiv);
 
     public class SBART_Report {
@@ -119,6 +175,114 @@
     }
 
+
+    public static IRegressionModel CalculateSBART(double[] x, double[] y, double[] w, int nKnots, string targetVariable, string[] inputVars, out SBART_Report rep) {
+      // use kMeans to find knot points
+      double[,] xy = new double[x.Length, 1];
+      for (int i = 0; i < x.Length; i++) xy[i, 0] = x[i];
+      double[,] c;
+      int[] xyc;
+      int info;
+      alglib.kmeansgenerate(xy, x.Length, 1, nKnots, 10, out info, out c, out xyc);
+      double[] knots = new double[nKnots];
+      for (int i = 0; i < knots.Length; i++) knots[i] = c[0, i];
+      return CalculateSBART(x, y, w, knots, targetVariable, inputVars, out rep);
+    }
+
+    public static IRegressionModel CalculateSBART(double[] x, double[] y, double[] w, double[] knots, string targetVariable, string[] inputVars, out SBART_Report rep) {
+      float[] xs = x.Select(xi=>(float)xi).ToArray();
+      float[] ys = y.Select(xi => (float)xi).ToArray();
+      float[] ws = w.Select(xi => (float)xi).ToArray();
+      float[] k = knots.Select(xi => (float)xi).ToArray();
+
+      int n = xs.Length;
+      if (n < 4) throw new ArgumentException("n < 4");
+      if (knots.Length > n + 2) throw new ArgumentException("more than n+2 knots");
+      float[] xw = new float[n];
+      int nx = -99;
+      float min = 0.0f;
+      float range = 0.0f;
+      int nk = -99;
+      float[] regKnots = new float[n + 6];
+
+      // sort xs together with ys and ws
+      // combine rows with duplicate x values
+      // transform x to range [0 .. 1]
+      // create a set of knots (using a heuristic for the number of knots)
+      // knots are located at data points. denser regions of x contain more knots.
+      SBART.setreg_x64(xs, ys, ws,
+        ref n, xw, ref nx, ref min, ref range, regKnots, ref nk);
+
+      // in this case we want to use the knots supplied by the caller.
+      // the knot values produced by setreg are overwritten with scaled knots supplied by caller.
+      // knots must be ordered as well.
+      int i = 0;
+      // left boundary
+      regKnots[i++] = 0.0f;
+      regKnots[i++] = 0.0f;
+      regKnots[i++] = 0.0f;
+      regKnots[i++] = 0.0f;
+      foreach (var knot in knots.OrderBy(ki=>ki)) {
+        regKnots[i++] = ((float)knot - min) / range;
+      }
+      // right boundary
+      regKnots[i++] = 1.0f;
+      regKnots[i++] = 1.0f;
+      regKnots[i++] = 1.0f;
+      regKnots[i++] = 1.0f;
+      nk = knots.Length + 4;
+
+      float criterion = -99.0f; // GCV
+      int icrit = 1; // calculate GCV
+      float smoothingParameter = -99.0f;
+      int smoothingParameterIndicator = 0;
+      float lowerSmoothingParameter = 0.5f;
+      float upperSmoothingParameter = 1.0f;
+      float tol = 0.01f;
+      int isetup = 0; // not setup?
+
+      // results
+      float[] coeff = new float[nk];
+      float[] leverage = new float[nx];
+      float[] y_smoothed = new float[nx];
+      int ier = -99;
+
+
+      // working arrays for sbart
+      float[] xwy = new float[nk];
+      float[] hs0 = new float[nk];
+      float[] hs1 = new float[nk];
+      float[] hs2 = new float[nk];
+      float[] hs3 = new float[nk];
+      float[] sg0 = new float[nk];
+      float[] sg1 = new float[nk];
+      float[] sg2 = new float[nk];
+      float[] sg3 = new float[nk];
+      int ld4 = 4;
+      float[,] adb = new float[ld4, nk];
+
+      float[,] p1ip = new float[nk, ld4];
+      int ldnk = nk + 4;
+      float[,] p2ip = new float[nk, nx];
+
+      SBART.sbart_x64(xs.Take(nx).ToArray(), ys.Take(nx).ToArray(), ws.Take(nx).ToArray(), ref nx,
+        regKnots, ref nk,
+        coeff, y_smoothed, leverage,
+        ref criterion, ref icrit,
+        ref smoothingParameter, ref smoothingParameterIndicator, ref lowerSmoothingParameter, ref upperSmoothingParameter,
+        ref tol, ref isetup,
+        xwy, hs0, hs1, hs2, hs3, sg0, sg1, sg2, sg3, adb, p1ip, p2ip, ref ld4, ref ldnk, ref ier);
+
+      if (ier > 0) throw new ArgumentException(ier.ToString());
+
+      rep = new SBART_Report();
+      rep.gcv = criterion;
+      rep.smoothingParameter = smoothingParameter;
+      rep.leverage = leverage.Select(li => (double)li).ToArray();
+
+      return new BartRegressionModel(regKnots.Take(nk+4).ToArray(), coeff, targetVariable, inputVars, min, range);
+    }
+
     public static IRegressionModel CalculateSBART(double[] x, double[] y,
-      string targetVariable, string[] inputVars, float smoothingParameter,
+      string targetVariable, string[] inputVars,
       out SBART_Report report) {
       var w = Enumerable.Repeat(1.0, x.Length).ToArray();
@@ -133,9 +297,9 @@
       int icrit = 1; // 0..don't calc CV,  1 .. GCV, 2 CV
 
-      // float smoothingParameter = -99.0f;
+      float smoothingParameter = -99.0f;
       int smoothingParameterIndicator = 0; 
       float lowerSmoothingParameter = 0f;
       float upperSmoothingParameter = 1.0f;
-      float tol = 0.01f;
+      float tol = 0.02f;
       int isetup = 0; // not setup?
 
@@ -149,4 +313,9 @@
         float[] ys = y.Select(yi => (float)yi).ToArray();
         float[] ws = w.Select(wi => (float)wi).ToArray();
+        
+        // sort xs together with ys and ws
+        // combine rows with duplicate x values
+        // create a set of knots (using a heuristic for the number of knots)
+        // knots are located at data points. denser regions of x contain more knots.
         SBART.setreg_x64(xs, ys, ws,
           ref n, xw, ref nx, ref min, ref range, knots, ref nk);
@@ -163,4 +332,5 @@
         */
 
+        /*
         // use uniform grid of knots
         nk = 20;
@@ -176,5 +346,5 @@
         knots[nk + 2] = xs[nx - 1];
         knots[nk + 3] = xs[nx - 1];
-
+        */
         if (nx < 4) {
           report = new SBART_Report();
@@ -220,5 +390,5 @@
         report.leverage = leverage.Select(li => (double)li).ToArray();
 
-        return new BartRegressionModel(xs.Take(nx).ToArray(), y_smoothed.Take(nx).ToArray(), targetVariable, inputVars, min, range);
+        return new BartRegressionModel(knots.Take(nk+4).ToArray(), coeff, targetVariable, inputVars, min, range);
 
       } else {
@@ -229,6 +399,6 @@
 
     public class BartRegressionModel : NamedItem, IRegressionModel {
-      private float[] x;
-      private float[] y;
+      private float[] knots;
+      private float[] bcoeff;
       private double min;
       private double range;
@@ -243,6 +413,6 @@
 
       public BartRegressionModel(BartRegressionModel orig, Cloner cloner) {
-        this.x = orig.x;
-        this.y = orig.y;
+        this.knots = orig.knots;
+        this.bcoeff = orig.bcoeff;
         this.min = orig.min;
         this.range = orig.range;
@@ -250,9 +420,9 @@
         this.variablesUsedForPrediction = orig.variablesUsedForPrediction;
       }
-      public BartRegressionModel(float[] x, float[] y_smoothed, string targetVariable, string[] inputVars, double min, double range) {
+      public BartRegressionModel(float[] knots, float[] bcoeff, string targetVariable, string[] inputVars, double min, double range) {
         this.variablesUsedForPrediction = inputVars;
         this.TargetVariable = targetVariable;
-        this.x = x;
-        this.y = y_smoothed;
+        this.knots = knots;
+        this.bcoeff = bcoeff;
         this.range = range;
         this.min = min;
@@ -266,32 +436,69 @@
 
       public double GetEstimatedValue(double xx) {
-        xx = (xx - min) / range;
-        // piecewise linear approximation
-        int n = x.Length;
-        if (xx <= x[0]) {
-          double h = xx - x[0];
-          return h * (y[1] - y[0]) / (x[1] - x[0]) + y[0];
-        } else if (xx >= x[n-1]) {
-          double h = xx - x[n-1];
-          return h * (y[n-1] - y[n-2]) / (x[n-1] - x[n-2]) + y[n-1];
-        } 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];
-              double k = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
-              return h * k + y[i];
-            } else if (xx < x[i]) {
-              upper = i - 1;
-            } else {
-              lower = i + 1;
-            }
-          }
-        }
-        return 0.0;
+        float x = (float)((xx - min) / range);
+        int n = bcoeff.Length;
+        int k = 4;
+        int zero = 0;
+        int one = 1;
+        int two = 2;
+
+
+        // linear extrapolation
+        if (x < 0) {
+          float x0 = 0.0f;
+          var y0 = bvalue(knots, bcoeff, ref n, ref k, ref x0, ref zero);
+          var y0d = bvalue(knots, bcoeff, ref n, ref k, ref x0, ref one);
+          return y0 + x * y0d;
+        }
+        if (x > 1) {
+          float x1 = 1.0f;
+          var y1 = bvalue(knots, bcoeff, ref n, ref k, ref x1, ref zero);
+          var y1d = bvalue(knots, bcoeff, ref n, ref k, ref x1, ref one);
+          return y1 + (x-1) * y1d;
+        }
+
+        lock (this) {
+          return bvalue(knots, bcoeff, ref n, ref k, ref x, ref zero);
+        }
+
+        // piecewise constant approximation
+        // if (xx <= x[0]) return bcoeff[0];
+        // if (xx >= x[n - 1]) return bcoeff[n - 1];
+        // for(int i=1;i<n-2;i++) {
+        //   var h1 = xx - x[i];
+        //   var h2 = xx - x[i + 1];
+        //   if(h1 > 0 && h2 <= 0) {
+        //     if (h1 < h2) return bcoeff[i]; else return bcoeff[i + 1];
+        //   }
+        // }
+        // return 0.0;
+
+        // // piecewise linear approximation
+        // int n = x.Length;
+        // if (xx <= x[0]) {
+        //   double h = xx - x[0];
+        //   return h * (y[1] - y[0]) / (x[1] - x[0]) + y[0];
+        // } else if (xx >= x[n-1]) {
+        //   double h = xx - x[n-1];
+        //   return h * (y[n-1] - y[n-2]) / (x[n-1] - x[n-2]) + y[n-1];
+        // } 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];
+        //       double k = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
+        //       return h * k + y[i];
+        //     } else if (xx < x[i]) {
+        //       upper = i - 1;
+        //     } else {
+        //       lower = i + 1;
+        //     }
+        //   }
+        // }
+        // return 0.0;
       }
 
Index: /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs	(revision 15467)
+++ /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs	(revision 15468)
@@ -148,11 +148,4 @@
       var validTypes = new ItemSet<StringValue>(
         new[] {
-        "Monotone", "Akima", "Catmull-Rom", "Cubic", "Linear",
-          "Cubic - Natural (Math.NET)",
-          "Polynomial (Math.NET)",
-          "Rational (Math.NET)",
-          "LogLinear (Math.NET)",
-          "Common (Math.NET)",
-          "Smoothing Spline (Mine)",
           "Smoothing Spline (Reinsch)",
           "Smoothing Spline (Reinsch with automatic tolerance determination using LOOCV)",
@@ -185,52 +178,4 @@
         var type = ((StringValue)Parameters["Type"].ActualValue).Value;
         switch (type) {
-          case "Monotone":
-            alglib.spline1dbuildmonotone(x, y, out c);
-            AddAlglibSplineResult(c, inputVars);
-            break;
-          case "Akima":
-            alglib.spline1dbuildakima(x, y, out c); AddAlglibSplineResult(c, inputVars);
-            break;
-            ;
-          case "Catmull-Rom":
-            alglib.spline1dbuildcatmullrom(x, y, out c); AddAlglibSplineResult(c, inputVars);
-            break;
-
-          case "Cubic":
-            alglib.spline1dbuildcubic(x, y, out c); AddAlglibSplineResult(c, inputVars);
-            break;
-
-          case "Linear":
-            alglib.spline1dbuildlinear(x, y, out c); AddAlglibSplineResult(c, inputVars);
-            break;
-          case "Cubic - Natural (Math.NET)": {
-              var spline = MathNet.Numerics.Interpolation.CubicSpline.InterpolateNatural(x, y);
-              AddMathNetSplineResult(spline, inputVars);
-              break;
-            }
-          case "Common (Math.NET)": {
-              var spline = MathNet.Numerics.Interpolate.Common(x, y);
-              AddMathNetSplineResult(spline, inputVars);
-              break;
-            }
-          case "LogLinear (Math.NET)": {
-              var spline = MathNet.Numerics.Interpolate.LogLinear(x, y);
-              AddMathNetSplineResult(spline, inputVars);
-              break;
-            }
-          case "Polynomial (Math.NET)": {
-              var spline = MathNet.Numerics.Interpolate.Polynomial(x, y);
-              AddMathNetSplineResult(spline, inputVars);
-              break;
-            }
-          case "Rational (Math.NET)": {
-              var spline = MathNet.Numerics.Interpolate.RationalWithoutPoles(x, y);
-              AddMathNetSplineResult(spline, inputVars);
-              break;
-            }
-          case "Smoothing Spline (Mine)": {
-              CalculateSmoothingSpline(x, y, inputVars);
-              break;
-            }
           case "Smoothing Spline (Reinsch)": {
               CalculateSmoothingSplineReinsch(x, y, inputVars);
@@ -264,8 +209,8 @@
           case "Finnbar O'Sullivan (sbart)": {
               SBART.SBART_Report report;
-              var lambda = ((IValueParameter<DoubleValue>)Parameters["Lambda"]).Value.Value; // controls extent of smoothing
+           
               var model =
                 SBART.CalculateSBART(x, y,
-                Problem.ProblemData.TargetVariable, inputVars, (float)lambda, out report);
+                Problem.ProblemData.TargetVariable, inputVars, out report);
               var targetVar = Problem.ProblemData.TargetVariable;
               var problemData = (IRegressionProblemData)Problem.ProblemData.Clone();
@@ -275,8 +220,4 @@
               break;
 
-            }
-          case "B-Spline Smoothing": {
-              BSplineSmoothing(x, y, inputVars);
-              break;
             }
           case "Penalized Regression Spline (alglib)": {
@@ -371,9 +312,4 @@
 
       return new RegressionEnsembleModel(models);
-    }
-
-    public void BSplineSmoothing(double[] x, double[] y, string[] inputVars) {
-      Array.Sort(x, y);
-      CubicSpline2(x, y, inputVars);
     }
 
@@ -414,64 +350,4 @@
     }
 
-    public void CubicSpline2(double[] x, double[] y, string[] inputVars) {
-      // see Pollock: Smoothing Splines
-      int n = x.Length;
-      double lambda = 1.0;
-      double[] sigma = Enumerable.Repeat(1.0, n).ToArray();
-      SplineData[] S = new SplineData[x.Length + 1];
-      for (int i = 0; i < x.Length; i++) {
-        S[i].x = x[i];
-        S[i].y = y[i];
-      }
-      S[x.Length].x = S[x.Length - 1].x;
-      S[x.Length].y = S[x.Length - 1].y;
-
-      // var
-      double[] h = new double[n],
-        r = new double[n],
-        f = new double[n],
-        p = new double[n];
-      var q = new MyArray<double>(-1, n + 3);
-      var u = new MyArray<double>(-1, n + 3);
-      var v = new MyArray<double>(-1, n + 3);
-      var w = new MyArray<double>(-1, n + 3);
-      double mu;
-
-      mu = 2 * (1 - lambda) / (3 * lambda);
-      h[0] = S[1].x - S[0].x;
-      r[0] = 3 / h[0];
-      for (int i = 1; i < n - 1; i++) {
-        h[i] = S[i + 1].x - S[i].x;
-        r[i] = 3 / h[i];
-        f[i] = -(r[i - 1] + r[i]);
-        p[i] = 2 * (S[i + 1].x - S[i - 1].x);
-        q[i] = 3 * (S[i + 1].y - S[i].y) / h[i] - 3 * (S[i].y - S[i - 1].y) / h[i - 1];
-      }
-      for (int i = 1; i < n; i++) {
-        u[i] = Math.Pow(r[i - 1], 2) * sigma[i - 1] + Math.Pow(f[i], 2) * sigma[i] + Math.Pow(r[i], 2) * sigma[i + 1];    // TODO Sigma 1..n
-        u[i] = mu * u[i] + p[i];
-        v[i] = f[i] * r[i] * sigma[i] + r[i] * f[i + 1] * sigma[i + 1];
-        v[i] = mu * v[i] + h[i];
-        w[i] = mu * r[i] * r[i + 1] * sigma[i + 1];
-      }
-      Quincunx(n, u, v, w, q);
-      // { Spline P arameters}
-      S[0].d = S[0].y - mu * r[0] * q[1] * sigma[0];
-      S[1].d = S[1].y - mu * (f[1] * q[1] + r[1] * q[2]) * sigma[0];
-      S[0].a = q[1] / (3 * h[0]);
-      S[0].b = 0;
-      S[0].c = (S[1].d - S[0].d) / h[0] - q[1] * h[0] / 3;
-      r[0] = 0;
-      for (int j = 1; j < n; j++) {
-        S[j].a = (q[j + 1] - q[j]) / (3 * h[j]);
-        S[j].b = q[j];
-        S[j].c = (q[j] + q[j - 1]) * h[j - 1] + S[j - 1].c;
-        S[j].d = r[j - 1] * q[j - 1] + f[j] * q[j] + r[j] * q[j + 1];
-        S[j].d = y[j] - mu * S[j].d * sigma[j];
-      }
-
-      //   { SmoothingSpline}
-    }
-
     public void CalculateSmoothingSplineReinsch(double[] x, double[] y, string[] inputVars) {
       double s = x.Length + 1;
@@ -541,5 +417,5 @@
 
 
-    // automatically determines tolerance using loo cv, SLOW!
+    // automatically determines tolerance using cv, SLOW!
     public static IRegressionModel CalculateSmoothingSplineReinsch(double[] x, double[] y, string[] inputVars, string targetVar,
       out double optTolerance, out double cvRMSE) {
@@ -665,62 +541,4 @@
 
       return new ReinschSmoothingSplineModel(a, b, c, d, x, targetVar, inputVars);
-    }
-
-    private void CalculateSmoothingSpline(double[] x, double[] y, string[] inputVars) {
-      // see Smoothing and Non-Parametric Regression, Germán Rodríguez, 2001 2.3.1
-      double[] w = Enumerable.Repeat(1.0, x.Length).ToArray();  // weights necessary for sortAndBin but are ignored below (TODO)
-      SortAndBin(x, y, w, out x, out y, out w, 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;
-        }
-      }
-
-      // WInvD = W.Cholesky().Solve(delta);
-      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())));
     }
 
@@ -766,14 +584,4 @@
       w2 = wl.ToArray();
     }
-
-    private void AddAlglibSplineResult(alglib.spline1dinterpolant c, string[] inputVars) {
-      Results.Add(new Result("Solution", new RegressionSolution(new Spline1dModel(c, Problem.ProblemData.TargetVariable, inputVars),
-        (IRegressionProblemData)Problem.ProblemData.Clone())));
-
-    }
-    private void AddMathNetSplineResult(IInterpolation c, string[] inputVars) {
-      Results.Add(new Result("Solution", new RegressionSolution(new MathNetSplineModel(c, Problem.ProblemData.TargetVariable, inputVars),
-        (IRegressionProblemData)Problem.ProblemData.Clone())));
-    }
   }
 
@@ -905,47 +713,6 @@
 
   // 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
   public class Spline1dModel : NamedItem, IRegressionModel {
-    private alglib.spline1dinterpolant interpolant;
+    alglib.spline1dinterpolant interpolant;
 
     public string TargetVariable { get; set; }
@@ -958,8 +725,11 @@
       this.TargetVariable = orig.TargetVariable;
       this.VariablesUsedForPrediction = orig.VariablesUsedForPrediction.ToArray();
+
       this.interpolant = (alglib.spline1dinterpolant)orig.interpolant.make_copy();
     }
-    public Spline1dModel(alglib.spline1dinterpolant interpolant, string targetVar, string[] inputs) : base("SplineModel", "SplineModel") {
-      this.interpolant = interpolant;
+    internal Spline1dModel(
+      alglib.spline1dinterpolant interpolant,
+      string targetVar, string[] inputs, double offset = 0, double scale = 1) : base("SplineModel", "SplineModel") {
+      this.interpolant = (alglib.spline1dinterpolant)interpolant.make_copy();
       this.TargetVariable = targetVar;
       this.VariablesUsedForPrediction = inputs;
@@ -972,4 +742,8 @@
     public IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
       return new RegressionSolution(this, (IRegressionProblemData)problemData.Clone());
+    }
+
+    public double GetEstimatedValue(double xx) {
+      return alglib.spline1dcalc(interpolant, xx);
     }
 
@@ -979,44 +753,6 @@
         product = product.Zip(dataset.GetDoubleValues(factor, rows), (pi, fi) => pi * fi).ToArray();
       }
-
-      foreach (var x in product) {
-        yield return alglib.spline1dcalc(interpolant, x);
-      }
-    }
-  }
-
-
-  // UNFINISHED
-  public class MathNetSplineModel : NamedItem, IRegressionModel {
-    private IInterpolation interpolant;
-
-    public string TargetVariable { get; set; }
-
-    public IEnumerable<string> VariablesUsedForPrediction { get; private set; }
-
-    public event EventHandler TargetVariableChanged;
-
-    public MathNetSplineModel(MathNetSplineModel orig, Cloner cloner) : base(orig, cloner) {
-      this.TargetVariable = orig.TargetVariable;
-      this.VariablesUsedForPrediction = orig.VariablesUsedForPrediction.ToArray();
-      this.interpolant = orig.interpolant; // TODO COPY!
-    }
-    public MathNetSplineModel(IInterpolation interpolant, string targetVar, string[] inputs) : base("SplineModel", "SplineModel") {
-      this.interpolant = interpolant;
-      this.TargetVariable = targetVar;
-      this.VariablesUsedForPrediction = inputs;
-    }
-
-    public override IDeepCloneable Clone(Cloner cloner) {
-      return new MathNetSplineModel(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 interpolant.Interpolate(x);
+      foreach (var xx in product) {
+        yield return GetEstimatedValue(xx);
       }
     }
