Index: /branches/HeuristicLab.Problems.GaussianProcessTuning/GaussianProcessDemo/Form1.cs
===================================================================
--- /branches/HeuristicLab.Problems.GaussianProcessTuning/GaussianProcessDemo/Form1.cs	(revision 10421)
+++ /branches/HeuristicLab.Problems.GaussianProcessTuning/GaussianProcessDemo/Form1.cs	(revision 10422)
@@ -4,4 +4,5 @@
 using System.Data;
 using System.Drawing;
+using System.Globalization;
 using System.Linq;
 using System.Text;
@@ -28,7 +29,7 @@
 
       var sum = new CovarianceSum();
-      var t = new CovarianceNeuralNetwork();
+      var t = new CovarianceSquaredExponentialIso();
       sum.Terms.Add(t);
-      sum.Terms.Add(new CovarianceNoise());
+      //sum.Terms.Add(new CovarianceNoise());
       this.covFunction = sum;
       UpdateSliders();
@@ -52,11 +53,11 @@
 
     private void InitData() {
-      int n = 200;
+      int n = 5000;
       data = new List<List<double>>();
-      data.Add(ValueGenerator.GenerateSteps(0, 1, 1.0 / n).ToList());
+      data.Add(Enumerable.Range(0, n).Select(e => e / 200.0).ToList());
 
       // sample from GP
       var normalRand = new NormalDistributedRandom(random, 0, 1);
-      alpha = (from i in Enumerable.Range(0, n + 1)
+      alpha = (from i in Enumerable.Range(0, n)
                select normalRand.NextDouble()).ToArray();
     }
@@ -72,31 +73,31 @@
       }
 
-      var trainingData = new List<List<double>>();
-      var trainingIndices = RandomEnumerable.SampleRandomWithoutRepetition(Enumerable.Range(0, y.Count), random, 10);
-      var trainingY = trainingIndices.Select(i => y[i]).ToList();
-      var trainingX = trainingIndices.Select(i => data[0][i]).ToList();
-      trainingData.Add(trainingY);
-      trainingData.Add(trainingX);
-      
-      //chart1.Series[2].Points.Clear();
-      //foreach (var p in trainingY.Zip(trainingX, (t, x) => new { t, x })) {
-      //  chart1.Series[2].Points.AddXY(p.x, p.t);
+      //var trainingData = new List<List<double>>();
+      //var trainingIndices = RandomEnumerable.SampleRandomWithoutRepetition(Enumerable.Range(0, y.Count), random, 10);
+      //var trainingY = trainingIndices.Select(i => y[i]).ToList();
+      //var trainingX = trainingIndices.Select(i => data[0][i]).ToList();
+      //trainingData.Add(trainingY);
+      //trainingData.Add(trainingX);
+      //
+      ////chart1.Series[2].Points.Clear();
+      ////foreach (var p in trainingY.Zip(trainingX, (t, x) => new { t, x })) {
+      ////  chart1.Series[2].Points.AddXY(p.x, p.t);
+      ////}
+      //
+      //var allData = new List<List<double>>();
+      //allData.Add(y);
+      //allData.Add(data[0]);
+      //var variableNames = new string[] { "y", "x" };
+      //var fullDataSet = new Dataset(variableNames, allData);
+      //var trainingDataSet = new Dataset(variableNames, trainingData);
+      //var trainingRows = Enumerable.Range(0, trainingIndices.Count());
+      //var fullRows = Enumerable.Range(0, data[0].Count);
+      //var correctModel = new GaussianProcessModel(fullDataSet, variableNames.First(), variableNames.Skip(1), fullRows, hyp, new MeanZero(),
+      //                                          (ICovarianceFunction)covFunction.Clone());
+      //var yPred = correctModel.GetEstimatedValues(fullDataSet, fullRows);
+      //chart1.Series[1].Points.Clear();
+      //foreach (var p in yPred.Zip(data[0], (t, x) => new { t, x })) {
+      //  chart1.Series[1].Points.AddXY(p.x, p.t);
       //}
-
-      var allData = new List<List<double>>();
-      allData.Add(y);
-      allData.Add(data[0]);
-      var variableNames = new string[] { "y", "x" };
-      var fullDataSet = new Dataset(variableNames, allData);
-      var trainingDataSet = new Dataset(variableNames, trainingData);
-      var trainingRows = Enumerable.Range(0, trainingIndices.Count());
-      var fullRows = Enumerable.Range(0, data[0].Count);
-      var correctModel = new GaussianProcessModel(fullDataSet, variableNames.First(), variableNames.Skip(1), fullRows, hyp, new MeanZero(),
-                                                (ICovarianceFunction)covFunction.Clone());
-      var yPred = correctModel.GetEstimatedValues(fullDataSet, fullRows);
-      chart1.Series[1].Points.Clear();
-      foreach (var p in yPred.Zip(data[0], (t, x) => new { t, x })) {
-        chart1.Series[1].Points.AddXY(p.x, p.t);
-      }
     }
 
@@ -124,6 +125,6 @@
     private string DataToText() {
       var str = new StringBuilder();
-      foreach (var p in chart1.Series[1].Points) {
-        str.AppendLine(p.XValue + "\t" + p.YValues.First());
+      foreach (var p in chart1.Series[0].Points) {
+        str.AppendLine(p.XValue.ToString(CultureInfo.InvariantCulture) + "\t" + p.YValues.First().ToString(CultureInfo.InvariantCulture));
       }
       return str.ToString();
Index: /branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Algorithms.DataAnalysis.Experimental/HeuristicLab.Algorithms.DataAnalysis.Experimental.csproj
===================================================================
--- /branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Algorithms.DataAnalysis.Experimental/HeuristicLab.Algorithms.DataAnalysis.Experimental.csproj	(revision 10421)
+++ /branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Algorithms.DataAnalysis.Experimental/HeuristicLab.Algorithms.DataAnalysis.Experimental.csproj	(revision 10422)
@@ -72,7 +72,7 @@
   </PropertyGroup>
   <ItemGroup>
-    <Reference Include="ALGLIB-3.6.0, Version=3.6.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=x86">
+    <Reference Include="ALGLIB-3.7.0, Version=3.7.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
       <SpecificVersion>False</SpecificVersion>
-      <HintPath>..\..\..\trunk\sources\bin\ALGLIB-3.6.0.dll</HintPath>
+      <HintPath>..\..\..\trunk\sources\bin\ALGLIB-3.7.0.dll</HintPath>
     </Reference>
     <Reference Include="HeuristicLab.Algorithms.DataAnalysis-3.4">
Index: /branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Instances.DataAnalysis.GaussianProcessRegression-3.3.csproj
===================================================================
--- /branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Instances.DataAnalysis.GaussianProcessRegression-3.3.csproj	(revision 10421)
+++ /branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Instances.DataAnalysis.GaussianProcessRegression-3.3.csproj	(revision 10422)
@@ -98,8 +98,7 @@
   </PropertyGroup>
   <ItemGroup>
-    <Reference Include="ALGLIB-3.6.0, Version=3.6.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
-      <SpecificVersion>False</SpecificVersion>
-      <HintPath>..\..\..\trunk\sources\bin\ALGLIB-3.6.0.dll</HintPath>
-      <Private>False</Private>
+    <Reference Include="ALGLIB-3.7.0, Version=3.7.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
+      <SpecificVersion>False</SpecificVersion>
+      <HintPath>..\..\..\trunk\sources\bin\ALGLIB-3.7.0.dll</HintPath>
     </Reference>
     <Reference Include="HeuristicLab.Algorithms.DataAnalysis-3.4, Version=3.4.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
Index: /branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Plugin.cs.frame
===================================================================
--- /branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Plugin.cs.frame	(revision 10421)
+++ /branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Plugin.cs.frame	(revision 10422)
@@ -25,5 +25,5 @@
   [Plugin("HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression", "3.3.7.$WCREV$")]
   [PluginFile("HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression-3.3.dll", PluginFileType.Assembly)]
-  [PluginDependency("HeuristicLab.ALGLIB", "3.6.0")]
+  [PluginDependency("HeuristicLab.ALGLIB", "3.7.0")]
   [PluginDependency("HeuristicLab.Algorithms.DataAnalysis", "3.4")]
   [PluginDependency("HeuristicLab.Collections", "3.3")]
Index: /branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Util.cs
===================================================================
--- /branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Util.cs	(revision 10421)
+++ /branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Util.cs	(revision 10422)
@@ -52,6 +52,6 @@
           K[i, j] = covFunction.Covariance(x, i, j);
 
-      if (!alglib.spdmatrixcholesky(ref K, K.GetLength(0), true)) throw new ArgumentException();
-
+      // if (!alglib.spdmatrixcholesky(ref K, K.GetLength(0), true)) throw new ArgumentException();
+      K = toeplitz_cholesky_lower(K.GetLength(0), K);
       List<double> target = new List<double>(K.GetLength(0));
       for (int i = 0; i < K.GetLength(0); i++) {
@@ -67,4 +67,201 @@
       return target;
     }
+
+    //****************************************************************************80
+
+    public static double[,] toeplitz_cholesky_lower(int n, double[,] a)
+
+    //****************************************************************************80
+      //
+      //  Purpose:
+      //
+      //    TOEPLITZ_CHOLESKY_LOWER: lower Cholesky factor of a Toeplitz matrix.
+      //
+      //  Discussion:
+      //
+      //    The Toeplitz matrix must be positive semi-definite.
+      //
+      //    After factorization, A = L * L'.
+      //
+      //  Licensing:
+      //
+      //    This code is distributed under the GNU LGPL license.
+      //
+      //  Modified:
+      //
+      //    13 November 2012
+      //    29 January  2014: adapted to C# by Gabriel Kronberger
+      //  Author:
+      //
+      //    John Burkardt
+      //
+      //  Reference:
+      //
+      //    Michael Stewart,
+      //    Cholesky factorization of semi-definite Toeplitz matrices.
+      //
+      //  Parameters:
+      //
+      //    Input, int N, the order of the matrix.
+      //
+      //    Input, double A[N,N], the Toeplitz matrix.
+      //
+      //    Output, double TOEPLITZ_CHOLESKY_LOWER[N,N], the lower Cholesky factor.
+      //
+    {
+      double div;
+      double[] g;
+      double g1j;
+      double g2j;
+      int i;
+      int j;
+      double[,] l;
+      double rho;
+
+      l = new double[n, n];
+
+      for (j = 0; j < n; j++) {
+        for (i = 0; i < n; i++) {
+          l[i, j] = 0.0;
+        }
+      }
+      g = new double[2 * n];
+
+      for (j = 0; j < n; j++) {
+        g[0 + j * 2] = a[0, j];
+      }
+      g[1 + 0 * 2] = 0.0;
+      for (j = 1; j < n; j++) {
+        g[1 + j * 2] = a[j, 0];
+      }
+
+      for (i = 0; i < n; i++) {
+        l[i, 0] = g[0 + i * 2];
+      }
+      for (j = n - 1; 1 <= j; j--) {
+        g[0 + j * 2] = g[0 + (j - 1) * 2];
+      }
+      g[0 + 0 * 2] = 0.0;
+
+      for (i = 1; i < n; i++) {
+        rho = -g[1 + i * 2] / g[0 + i * 2];
+        div = Math.Sqrt((1.0 - rho) * (1.0 + rho));
+
+        for (j = i; j < n; j++) {
+          g1j = g[0 + j * 2];
+          g2j = g[1 + j * 2];
+          g[0 + j * 2] = (g1j + rho * g2j) / div;
+          g[1 + j * 2] = (rho * g1j + g2j) / div;
+        }
+
+        for (j = i; j < n; j++) {
+          l[j, i] = g[0 + j * 2];
+        }
+        for (j = n - 1; i < j; j--) {
+          g[0 + j * 2] = g[0 + (j - 1) * 2];
+        }
+        g[0 + i * 2] = 0.0;
+      }
+
+
+      return l;
+    }
+    //****************************************************************************80
+
+    public static double[,] toeplitz_cholesky_upper(int n, double[,] a)
+
+    //****************************************************************************80
+      //
+      //  Purpose:
+      //
+      //    TOEPLITZ_CHOLESKY_UPPER: upper Cholesky factor of a Toeplitz matrix.
+      //
+      //  Discussion:
+      //
+      //    The Toeplitz matrix must be positive semi-definite.
+      //
+      //    After factorization, A = R' * R.
+      //
+      //  Licensing:
+      //
+      //    This code is distributed under the GNU LGPL license.
+      //
+      //  Modified:
+      //
+      //    14 November 2012
+      //    29 January  2014: adapted to C# by Gabriel Kronberger
+      //
+      //  Author:
+      //
+      //    John Burkardt
+      //
+      //  Reference:
+      //
+      //    Michael Stewart,
+      //    Cholesky factorization of semi-definite Toeplitz matrices.
+      //
+      //  Parameters:
+      //
+      //    Input, int N, the order of the matrix.
+      //
+      //    Input, double A[N,N], the Toeplitz matrix.
+      //
+      //    Output, double TOEPLITZ_CHOLESKY_UPPER[N,N], the upper Cholesky factor.
+      //
+    {
+      double div;
+      double[] g;
+      double g1j;
+      double g2j;
+      int i;
+      int j;
+      double[,] r;
+      double rho;
+
+      r = new double[n, n];
+
+      for (j = 0; j < n; j++) {
+        for (i = 0; i < n; i++) {
+          r[i, j] = 0.0;
+        }
+      }
+      g = new double[2 * n];
+
+      for (j = 0; j < n; j++) {
+        g[0 + j * 2] = a[0, j];
+      }
+
+      g[1 + 0 * 2] = 0.0;
+      for (j = 1; j < n; j++) {
+        g[1 + j * 2] = a[j, 0];
+      }
+      for (j = 0; j < n; j++) {
+        r[0, j] = g[0 + j * 2];
+      }
+      for (j = n - 1; 1 <= j; j--) {
+        g[0 + j * 2] = g[0 + (j - 1) * 2];
+      }
+      g[0 + 0 * 2] = 0.0;
+
+      for (i = 1; i < n; i++) {
+        rho = -g[1 + i * 2] / g[0 + i * 2];
+        div = Math.Sqrt((1.0 - rho) * (1.0 + rho));
+        for (j = i; j < n; j++) {
+          g1j = g[0 + j * 2];
+          g2j = g[1 + j * 2];
+          g[0 + j * 2] = (g1j + rho * g2j) / div;
+          g[1 + j * 2] = (rho * g1j + g2j) / div;
+        }
+        for (j = i; j < n; j++) {
+          r[i, j] = g[0 + j * 2];
+        }
+        for (j = n - 1; i < j; j--) {
+          g[0 + j * 2] = g[0 + (j - 1) * 2];
+        }
+        g[0 + i * 2] = 0.0;
+      }
+
+      return r;
+    }
   }
 }
