Index: /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs	(revision 15468)
+++ /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs	(revision 15469)
@@ -223,12 +223,12 @@
           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
-          );
-
-        double optTolerance; double cvRMSE;
+        // CubicSplineGCV.CubGcvReport report;
+        // return CubicSplineGCV.CalculateCubicSpline(
+        //   product,
+        //   (double[])target.Clone(), 
+        //   problemData.TargetVariable, inputVars, out report
+        //   );
+        // 
+        // double optTolerance; double cvRMSE;
         // find tolerance
         // var ensemble = Splines.CalculateSmoothingSplineReinsch(product, (double[])target.Clone(), inputVars, problemData.TargetVariable, out optTolerance, out cvRMSE);
@@ -239,7 +239,7 @@
         // find tolerance
         //var bestLambda = double.NaN;
-        double bestCVRMSE = target.StandardDeviation();
-        double avgTrainRMSE = double.PositiveInfinity;
-        double[] bestPredictions = new double[target.Length]; // zero
+        // double bestCVRMSE = target.StandardDeviation();
+        // double avgTrainRMSE = double.PositiveInfinity;
+        // double[] bestPredictions = new double[target.Length]; // zero
 
 
@@ -269,5 +269,6 @@
         // 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, out rep);
+        var w = product.Select(_ => 1.0).ToArray();
+        var model = SBART.CalculateSBART(product, (double[])target.Clone(), w, 10, 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 15468)
+++ /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/SBART.cs	(revision 15469)
@@ -184,101 +184,117 @@
       int info;
       alglib.kmeansgenerate(xy, x.Length, 1, nKnots, 10, out info, out c, out xyc);
+      var g = x.Zip(xyc, (double xi, int ci) => Tuple.Create(xi,ci)).GroupBy(t => t.Item2).Select(gr => HeuristicLab.Common.EnumerableStatisticExtensions.Median(gr.Select(gi=>gi.Item1))).ToArray();
+      // find medians
       double[] knots = new double[nKnots];
-      for (int i = 0; i < knots.Length; i++) knots[i] = c[0, i];
+      for (int i = 0; i < g.Length; i++) knots[i] = g[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;
+      int ier = 99;
+      int tries = 0;
       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);
+
+      do {
+        tries++;
+        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;
+        int j = 1;
+        foreach (var knot in knots.OrderBy(ki => ki)) {
+          regKnots[i++] = xs[j * nx / (knots.Length + 1)];  // ((float)knot - min) / range;
+          j++;
+        }
+        // right boundary
+        regKnots[i++] = 1.0f;
+        regKnots[i++] = 1.0f;
+        regKnots[i++] = 1.0f;
+        regKnots[i++] = 1.0f;
+        nk = i - 4;
+
+        float criterion = -99.0f; // GCV
+        int icrit = 1; // calculate GCV
+        float smoothingParameter = -99.0f;
+        int smoothingParameterIndicator = 0;
+        float lowerSmoothingParameter = 0.0f;
+        float upperSmoothingParameter = 1.0f;
+        int isetup = 0; // not setup?
+
+        // results
+        float[] coeff = new float[nk];
+        float[] leverage = new float[nx];
+        float[] y_smoothed = new float[nx];
+
+        // 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) {
+          Console.WriteLine("ERROR {0} smooth {1}  criterion {2}", ier, smoothingParameter, criterion);
+          tol *= 2;
+          tol = Math.Min(tol, 1.0f);
+        } else {
+          if (tries > 1) {
+            Console.WriteLine("Success {0} smooth {1}  criterion {2}", ier, smoothingParameter, criterion);
+          }
+          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);
+        }
+      } while (ier > 0);
+      throw new ArgumentException();
     }
 
Index: /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs	(revision 15468)
+++ /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs	(revision 15469)
@@ -223,5 +223,5 @@
           case "Penalized Regression Spline (alglib)": {
               var lambda = ((IValueParameter<DoubleValue>)Parameters["Lambda"]).Value.Value; // controls extent of smoothing
-              var model = CalculatePenalizedRegressionSpline(x, y, lambda, Problem.ProblemData.TargetVariable, inputVars);
+              var model = CalculatePenalizedRegressionSpline(x, y, lambda, x.Length, Problem.ProblemData.TargetVariable, inputVars);
               var targetVar = Problem.ProblemData.TargetVariable;
               var problemData = (IRegressionProblemData)Problem.ProblemData.Clone();
@@ -237,12 +237,10 @@
 
 
-    public static IRegressionModel CalculatePenalizedRegressionSpline(double[] x, double[] y, double lambda, string targetVar, string[] inputVars) {
+    public static IRegressionModel CalculatePenalizedRegressionSpline(double[] x, double[] y, double lambda, int numKnots, string targetVar, string[] inputVars) {
       int info;
       alglib.spline1dinterpolant interpolant;
       alglib.spline1dfitreport rep;
-      int n = x.Length;
-      n = (int)Math.Max(50, 3 * Math.Sqrt(n));
-
-      alglib.spline1dfitpenalized(x, y, n, lambda, out info, out interpolant, out rep);
+
+      alglib.spline1dfitpenalized(x, y, numKnots, lambda, out info, out interpolant, out rep);
       return new Spline1dModel(interpolant, targetVar, inputVars);
     }
@@ -728,5 +726,5 @@
       this.interpolant = (alglib.spline1dinterpolant)orig.interpolant.make_copy();
     }
-    internal Spline1dModel(
+    public Spline1dModel(
       alglib.spline1dinterpolant interpolant,
       string targetVar, string[] inputs, double offset = 0, double scale = 1) : base("SplineModel", "SplineModel") {
Index: /branches/MathNetNumerics-Exploration-2789/Main/Program.cs
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/Main/Program.cs	(revision 15468)
+++ /branches/MathNetNumerics-Exploration-2789/Main/Program.cs	(revision 15469)
@@ -13,13 +13,14 @@
     static void Main(string[] args) {
 
-      var xs = HeuristicLab.Common.SequenceGenerator.GenerateSteps(-3.5, 3.5, 0.02, includeEnd: true).ToArray();
-      var ys = xs.Select(xi => (1.0 / (Math.Sqrt(2 * Math.PI)) * Math.Exp(-0.5 * xi * xi))).ToArray();
-
-      alglib.hqrndstate state;
-      alglib.hqrndseed(1234, 5678, out state);
-      var ys_noise = ys.Select(yi => yi + alglib.hqrndnormal(state) * 0.01).ToArray();
-
-      SBART.SBART_Report rep;
-      SBART.CalculateSBART(xs, ys_noise, "y", new string[] { "x" }, 1.0f, out rep);
+       var xs = HeuristicLab.Common.SequenceGenerator.GenerateSteps(-3.5, 3.5, 0.02, includeEnd: true).ToArray();
+       var ys = xs.Select(xi => (1.0 / (Math.Sqrt(2 * Math.PI)) * Math.Exp(-0.5 * xi * xi))).ToArray();
+       var ws = xs.Select(_ => 1.0).ToArray();
+ 
+       alglib.hqrndstate state;
+       alglib.hqrndseed(1234, 5678, out state);
+       var ys_noise = ys.Select(yi => yi + alglib.hqrndnormal(state) * 0.01).ToArray();
+ 
+       SBART.SBART_Report rep;
+       SBART.CalculateSBART(xs, ys_noise, ws, 2, "y", new string[] { "x" }, out rep);
 
 
Index: /branches/MathNetNumerics-Exploration-2789/Test/UnitTest1.cs
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/Test/UnitTest1.cs	(revision 15468)
+++ /branches/MathNetNumerics-Exploration-2789/Test/UnitTest1.cs	(revision 15469)
@@ -46,5 +46,6 @@
         double loo_rmse;
         double avgTrainRmse;
-        Splines.CalculatePenalizedRegressionSpline(xs, ys_noise, rho, "y", new string[] { "x" }, out avgTrainRmse, out loo_rmse);
+        double[] predictions;
+        Splines.CalculatePenalizedRegressionSpline(xs, ys_noise, rho, "y", new string[] { "x" }, out avgTrainRmse, out loo_rmse, out predictions);
         iters++;
         Console.WriteLine("{0} {1} {2}", rho, avgTrainRmse, loo_rmse);
