Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GBM/GradientBoostingRegressionAlgorithm.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GBM/GradientBoostingRegressionAlgorithm.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GBM/GradientBoostingRegressionAlgorithm.cs	(revision 13981)
@@ -343,4 +343,5 @@
       var features = models.Count;
       double[,] inputMatrix = new double[rows, features + 1];
+
       //add model estimates
       for (int m = 0; m < models.Count; m++) {
@@ -372,5 +373,5 @@
 
       var ensembleModel = new RegressionEnsembleModel(models, coefficients.Take(models.Count)) { AverageModelEstimates = false };
-      var ensembleSolution = ensembleModel.CreateRegressionSolution(problemData);
+      var ensembleSolution = ensembleModel.CreateRegressionSolution((IRegressionProblemData)problemData.Clone());
       return ensembleSolution;
     }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceConst.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceConst.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceConst.cs	(revision 13981)
@@ -21,6 +21,4 @@
 
 using System;
-using System.Collections.Generic;
-using System.Linq;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -83,5 +81,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double scale;
       GetParameterValues(p, out scale);
@@ -91,14 +89,10 @@
       cov.CrossCovariance = (x, xt, i, j) => scale;
       if (HasFixedScaleParameter) {
-        cov.CovarianceGradient = (x, i, j) => Enumerable.Empty<double>();
+        cov.CovarianceGradient = (x, i, j) => new double[0];
       } else {
-        cov.CovarianceGradient = (x, i, j) => GetGradient(x, i, j, scale, columnIndices);
+        cov.CovarianceGradient = (x, i, j) => new[] { 2.0 * scale };
       }
       return cov;
     }
-
-    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, double scale, IEnumerable<int> columnIndices) {
-      yield return 2.0 * scale;
-    }
   }
 }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceLinear.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceLinear.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceLinear.cs	(revision 13981)
@@ -21,6 +21,4 @@
 
 using System;
-using System.Collections.Generic;
-using System.Linq;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -52,11 +50,11 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       if (p.Length > 0) throw new ArgumentException("No parameters are allowed for the linear covariance function.");
       // create functions
       var cov = new ParameterizedCovarianceFunction();
-      cov.Covariance = (x, i, j) => Util.ScalarProd(x, i, j, 1, columnIndices);
-      cov.CrossCovariance = (x, xt, i, j) =>  Util.ScalarProd(x, i, xt, j, 1.0 , columnIndices);
-      cov.CovarianceGradient = (x, i, j) => Enumerable.Empty<double>();
+      cov.Covariance = (x, i, j) => Util.ScalarProd(x, i, j, columnIndices, 1.0);
+      cov.CrossCovariance = (x, xt, i, j) => Util.ScalarProd(x, i, xt, j, columnIndices, 1.0);
+      cov.CovarianceGradient = (x, i, j) => new double[0];
       return cov;
     }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceLinearArd.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceLinearArd.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceLinearArd.cs	(revision 13981)
@@ -81,5 +81,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double[] inverseLength;
       GetParameterValues(p, out inverseLength);
@@ -90,5 +90,5 @@
       cov.CrossCovariance = (x, xt, i, j) => Util.ScalarProd(x, i, xt, j, inverseLength, columnIndices);
       if (fixedInverseLength)
-        cov.CovarianceGradient = (x, i, j) => Enumerable.Empty<double>();
+        cov.CovarianceGradient = (x, i, j) => new double[0];
       else
         cov.CovarianceGradient = (x, i, j) => GetGradient(x, i, j, inverseLength, columnIndices);
@@ -96,10 +96,13 @@
     }
 
-    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, double[] inverseLength, IEnumerable<int> columnIndices) {
+    private static IList<double> GetGradient(double[,] x, int i, int j, double[] inverseLength, int[] columnIndices) {
       int k = 0;
-      foreach (int columnIndex in columnIndices) {
-        yield return -2.0 * x[i, columnIndex] * x[j, columnIndex] * inverseLength[k] * inverseLength[k];
+      var g = new List<double>(columnIndices.Length);
+      for (int c = 0; c < columnIndices.Length; c++) {
+        var columnIndex = columnIndices[c];
+        g.Add(-2.0 * x[i, columnIndex] * x[j, columnIndex] * inverseLength[k] * inverseLength[k]);
         k++;
       }
+      return g;
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceMask.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceMask.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceMask.cs	(revision 13981)
@@ -20,8 +20,5 @@
 #endregion
 
-using System;
-using System.Collections.Generic;
 using System.Linq;
-using System.Linq.Expressions;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -73,9 +70,9 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       var cov = CovarianceFunctionParameter.Value;
       var selectedDimensions = SelectedDimensionsParameter.Value;
 
-      return cov.GetParameterizedCovarianceFunction(p, selectedDimensions.Intersect(columnIndices));
+      return cov.GetParameterizedCovarianceFunction(p, selectedDimensions.Intersect(columnIndices).ToArray());
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceMaternIso.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceMaternIso.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceMaternIso.cs	(revision 13981)
@@ -111,5 +111,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double inverseLength, scale;
       int d = DParameter.Value.Value;
@@ -122,9 +122,9 @@
         double dist = i == j
                        ? 0.0
-                       : Math.Sqrt(Util.SqrDist(x, i, j, Math.Sqrt(d) * inverseLength, columnIndices));
+                       : Math.Sqrt(Util.SqrDist(x, i, j, columnIndices, Math.Sqrt(d) * inverseLength));
         return scale * m(d, dist);
       };
       cov.CrossCovariance = (x, xt, i, j) => {
-        double dist = Math.Sqrt(Util.SqrDist(x, i, xt, j, Math.Sqrt(d) * inverseLength, columnIndices));
+        double dist = Math.Sqrt(Util.SqrDist(x, i, xt, j, columnIndices, Math.Sqrt(d) * inverseLength));
         return scale * m(d, dist);
       };
@@ -155,13 +155,14 @@
     }
 
-
-    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, int d, double scale, double inverseLength, IEnumerable<int> columnIndices,
+    private static IList<double> GetGradient(double[,] x, int i, int j, int d, double scale, double inverseLength, int[] columnIndices,
       bool fixedInverseLength, bool fixedScale) {
       double dist = i == j
                    ? 0.0
-                   : Math.Sqrt(Util.SqrDist(x, i, j, Math.Sqrt(d) * inverseLength, columnIndices));
+                   : Math.Sqrt(Util.SqrDist(x, i, j, columnIndices, Math.Sqrt(d) * inverseLength));
 
-      if (!fixedInverseLength) yield return scale * dm(d, dist);
-      if (!fixedScale) yield return 2 * scale * m(d, dist);
+      var g = new List<double>(2);
+      if (!fixedInverseLength) g.Add(scale * dm(d, dist));
+      if (!fixedScale) g.Add(2 * scale * m(d, dist));
+      return g;
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceNeuralNetwork.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceNeuralNetwork.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceNeuralNetwork.cs	(revision 13981)
@@ -102,5 +102,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double length, scale;
       GetParameterValues(p, out scale, out length);
@@ -113,5 +113,6 @@
         double s1 = 1.0;
         double s2 = 1.0;
-        foreach (var col in columnIndices) {
+        for (int c = 0; c < columnIndices.Length; c++) {
+          var col = columnIndices[c];
           sx += x[i, col] * x[j, col];
           s1 += x[i, col] * x[i, col];
@@ -125,5 +126,6 @@
         double s1 = 1.0;
         double s2 = 1.0;
-        foreach (var col in columnIndices) {
+        for (int c = 0; c < columnIndices.Length; c++) {
+          var col = columnIndices[c];
           sx += x[i, col] * xt[j, col];
           s1 += x[i, col] * x[i, col];
@@ -138,24 +140,22 @@
 
     // order of returned gradients must match the order in GetParameterValues!
-    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, double length, double scale, IEnumerable<int> columnIndices,
+    private static IList<double> GetGradient(double[,] x, int i, int j, double length, double scale, int[] columnIndices,
       bool fixedLength, bool fixedScale) {
-      {
-        double sx = 1.0;
-        double s1 = 1.0;
-        double s2 = 1.0;
-        foreach (var col in columnIndices) {
-          sx += x[i, col] * x[j, col];
-          s1 += x[i, col] * x[i, col];
-          s2 += x[j, col] * x[j, col];
-        }
-        var h = (length + s1) * (length + s2);
-        var f = sx / Math.Sqrt(h);
-        if (!fixedLength) {
-          yield return -scale / Math.Sqrt(1.0 - f * f) * ((length * sx * (2.0 * length + s1 + s2)) / Math.Pow(h, 3.0 / 2.0));
-        }
-        if (!fixedScale) {
-          yield return 2.0 * scale * Math.Asin(f);
-        }
+      double sx = 1.0;
+      double s1 = 1.0;
+      double s2 = 1.0;
+      for (int c = 0; c < columnIndices.Length; c++) {
+        var col = columnIndices[c];
+        sx += x[i, col] * x[j, col];
+        s1 += x[i, col] * x[i, col];
+        s2 += x[j, col] * x[j, col];
       }
+      var h = (length + s1) * (length + s2);
+      var f = sx / Math.Sqrt(h);
+
+      var g = new List<double>(2);
+      if (!fixedLength) g.Add(-scale / Math.Sqrt(1.0 - f * f) * ((length * sx * (2.0 * length + s1 + s2)) / Math.Pow(h, 3.0 / 2.0)));
+      if (!fixedScale) g.Add(2.0 * scale * Math.Asin(f));
+      return g;
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceNoise.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceNoise.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceNoise.cs	(revision 13981)
@@ -21,6 +21,4 @@
 
 using System;
-using System.Collections.Generic;
-using System.Linq;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -84,5 +82,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double scale;
       GetParameterValues(p, out scale);
@@ -91,9 +89,9 @@
       var cov = new ParameterizedCovarianceFunction();
       cov.Covariance = (x, i, j) => i == j ? scale : 0.0;
-      cov.CrossCovariance = (x, xt, i, j) => Util.SqrDist(x, i, xt, j, 1.0, columnIndices) < 1e-9 ? scale : 0.0;
+      cov.CrossCovariance = (x, xt, i, j) => Util.SqrDist(x, i, xt, j, columnIndices, 1.0) < 1e-9 ? scale : 0.0;
       if (fixedScale)
-        cov.CovarianceGradient = (x, i, j) => Enumerable.Empty<double>();
+        cov.CovarianceGradient = (x, i, j) => new double[0];
       else
-        cov.CovarianceGradient = (x, i, j) => Enumerable.Repeat(i == j ? 2.0 * scale : 0.0, 1);
+        cov.CovarianceGradient = (x, i, j) => new double[1] { i == j ? 2.0 * scale : 0.0 };
       return cov;
     }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovariancePeriodic.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovariancePeriodic.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovariancePeriodic.cs	(revision 13981)
@@ -22,5 +22,4 @@
 using System;
 using System.Collections.Generic;
-using System.Linq;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -117,5 +116,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double inverseLength, period, scale;
       GetParameterValues(p, out scale, out period, out inverseLength);
@@ -145,22 +144,23 @@
     }
 
-
-    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, IEnumerable<int> columnIndices, double scale, double period, double inverseLength,
+    private static IList<double> GetGradient(double[,] x, int i, int j, int[] columnIndices, double scale, double period, double inverseLength,
       bool fixedInverseLength, bool fixedPeriod, bool fixedScale) {
       double k = i == j ? 0.0 : Math.PI * GetDistance(x, x, i, j, columnIndices) / period;
       double gradient = Math.Sin(k) * inverseLength;
       gradient *= gradient;
-      if (!fixedInverseLength) yield return 4.0 * scale * Math.Exp(-2.0 * gradient) * gradient;
+      var g = new List<double>(3);
+      if (!fixedInverseLength) 
+        g.Add(4.0 * scale * Math.Exp(-2.0 * gradient) * gradient);
       if (!fixedPeriod) {
         double r = Math.Sin(k) * inverseLength;
-        yield return 2.0 * k * scale * Math.Exp(-2 * r * r) * Math.Sin(2 * k) * inverseLength * inverseLength;
+        g.Add(2.0 * k * scale * Math.Exp(-2 * r * r) * Math.Sin(2 * k) * inverseLength * inverseLength);
       }
       if (!fixedScale)
-        yield return 2.0 * scale * Math.Exp(-2 * gradient);
-
+        g.Add(2.0 * scale * Math.Exp(-2 * gradient));
+      return g;
     }
 
-    private static double GetDistance(double[,] x, double[,] xt, int i, int j, IEnumerable<int> columnIndices) {
-      return Math.Sqrt(Util.SqrDist(x, i, xt, j, 1, columnIndices));
+    private static double GetDistance(double[,] x, double[,] xt, int i, int j, int[] columnIndices) {
+      return Math.Sqrt(Util.SqrDist(x, i, xt, j, columnIndices, 1));
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovariancePiecewisePolynomial.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovariancePiecewisePolynomial.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovariancePiecewisePolynomial.cs	(revision 13981)
@@ -69,8 +69,8 @@
       Parameters.Add(new OptionalValueParameter<DoubleValue>("Scale", "The scale parameter of the piecewise polynomial covariance function."));
 
-      var validValues = new ItemSet<IntValue>(new IntValue[] { 
-        (IntValue)(new IntValue().AsReadOnly()), 
-        (IntValue)(new IntValue(1).AsReadOnly()), 
-        (IntValue)(new IntValue(2).AsReadOnly()), 
+      var validValues = new ItemSet<IntValue>(new IntValue[] {
+        (IntValue)(new IntValue().AsReadOnly()),
+        (IntValue)(new IntValue(1).AsReadOnly()),
+        (IntValue)(new IntValue(2).AsReadOnly()),
         (IntValue)(new IntValue(3).AsReadOnly()) });
       Parameters.Add(new ConstrainedValueParameter<IntValue>("V", "The v parameter of the piecewise polynomial function (allowed values 0, 1, 2, 3).", validValues, validValues.First()));
@@ -113,5 +113,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double length, scale;
       int v = VParameter.Value.Value;
@@ -148,9 +148,9 @@
       var cov = new ParameterizedCovarianceFunction();
       cov.Covariance = (x, i, j) => {
-        double k = Math.Sqrt(Util.SqrDist(x, i, x, j, 1.0 / length, columnIndices));
+        double k = Math.Sqrt(Util.SqrDist(x, i, x, j, columnIndices, 1.0 / length));
         return scale * Math.Pow(Math.Max(1 - k, 0), exp + v) * f(k);
       };
       cov.CrossCovariance = (x, xt, i, j) => {
-        double k = Math.Sqrt(Util.SqrDist(x, i, xt, j, 1.0 / length, columnIndices));
+        double k = Math.Sqrt(Util.SqrDist(x, i, xt, j, columnIndices, 1.0 / length));
         return scale * Math.Pow(Math.Max(1 - k, 0), exp + v) * f(k);
       };
@@ -159,9 +159,11 @@
     }
 
-    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, double length, double scale, int v, double exp, Func<double, double> f, Func<double, double> df, IEnumerable<int> columnIndices,
+    private static IList<double> GetGradient(double[,] x, int i, int j, double length, double scale, int v, double exp, Func<double, double> f, Func<double, double> df, int[] columnIndices,
       bool fixedLength, bool fixedScale) {
-      double k = Math.Sqrt(Util.SqrDist(x, i, x, j, 1.0 / length, columnIndices));
-      if (!fixedLength) yield return scale * Math.Pow(Math.Max(1.0 - k, 0), exp + v - 1) * k * ((exp + v) * f(k) - Math.Max(1 - k, 0) * df(k));
-      if (!fixedScale) yield return 2.0 * scale * Math.Pow(Math.Max(1 - k, 0), exp + v) * f(k);
+      double k = Math.Sqrt(Util.SqrDist(x, i, x, j, columnIndices, 1.0 / length));
+      var g = new List<double>(2);
+      if (!fixedLength) g.Add(scale * Math.Pow(Math.Max(1.0 - k, 0), exp + v - 1) * k * ((exp + v) * f(k) - Math.Max(1 - k, 0) * df(k)));
+      if (!fixedScale) g.Add(2.0 * scale * Math.Pow(Math.Max(1 - k, 0), exp + v) * f(k));
+      return g;
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovariancePolynomial.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovariancePolynomial.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovariancePolynomial.cs	(revision 13981)
@@ -22,5 +22,4 @@
 using System;
 using System.Collections.Generic;
-using System.Linq;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -107,5 +106,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double @const, scale;
       int degree = DegreeParameter.Value.Value;
@@ -116,15 +115,17 @@
       // create functions
       var cov = new ParameterizedCovarianceFunction();
-      cov.Covariance = (x, i, j) => scale * Math.Pow(@const + Util.ScalarProd(x, i, j, 1.0, columnIndices), degree);
-      cov.CrossCovariance = (x, xt, i, j) => scale * Math.Pow(@const + Util.ScalarProd(x, i, xt, j, 1.0, columnIndices), degree);
+      cov.Covariance = (x, i, j) => scale * Math.Pow(@const + Util.ScalarProd(x, i, j, columnIndices, 1.0), degree);
+      cov.CrossCovariance = (x, xt, i, j) => scale * Math.Pow(@const + Util.ScalarProd(x, i, xt, j, columnIndices, 1.0), degree);
       cov.CovarianceGradient = (x, i, j) => GetGradient(x, i, j, @const, scale, degree, columnIndices, fixedConst, fixedScale);
       return cov;
     }
 
-    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, double c, double scale, int degree, IEnumerable<int> columnIndices,
+    private static IList<double> GetGradient(double[,] x, int i, int j, double c, double scale, int degree, int[] columnIndices,
       bool fixedConst, bool fixedScale) {
-      double s = Util.ScalarProd(x, i, j, 1.0, columnIndices);
-      if (!fixedConst) yield return c * degree * scale * Math.Pow(c + s, degree - 1);
-      if (!fixedScale) yield return 2 * scale * Math.Pow(c + s, degree);
+      double s = Util.ScalarProd(x, i, j, columnIndices, 1.0);
+      var g = new List<double>(2);
+      if (!fixedConst) g.Add(c * degree * scale * Math.Pow(c + s, degree - 1));
+      if (!fixedScale) g.Add(2 * scale * Math.Pow(c + s, degree));
+      return g;
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceProduct.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceProduct.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceProduct.cs	(revision 13981)
@@ -23,5 +23,4 @@
 using System.Collections.Generic;
 using System.Linq;
-using System.Linq.Expressions;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -76,5 +75,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       if (factors.Count == 0) throw new ArgumentException("at least one factor is necessary for the product covariance function.");
       var functions = new List<ParameterizedCovarianceFunction>();
@@ -93,6 +92,7 @@
     }
 
-    public static IEnumerable<double> GetGradient(double[,] x, int i, int j, List<ParameterizedCovarianceFunction> factorFunctions) {
+    public static IList<double> GetGradient(double[,] x, int i, int j, List<ParameterizedCovarianceFunction> factorFunctions) {
       var covariances = factorFunctions.Select(f => f.Covariance(x, i, j)).ToArray();
+      var gr = new List<double>();
       for (int ii = 0; ii < factorFunctions.Count; ii++) {
         foreach (var g in factorFunctions[ii].CovarianceGradient(x, i, j)) {
@@ -100,7 +100,8 @@
           for (int jj = 0; jj < covariances.Length; jj++)
             if (ii != jj) res *= covariances[jj];
-          yield return res;
+          gr.Add(res);
         }
       }
+      return gr;
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceRationalQuadraticArd.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceRationalQuadraticArd.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceRationalQuadraticArd.cs	(revision 13981)
@@ -121,5 +121,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double scale, shape;
       double[] inverseLength;
@@ -144,5 +144,5 @@
     }
 
-    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, IEnumerable<int> columnIndices, double scale, double shape, double[] inverseLength,
+    private static IList<double> GetGradient(double[,] x, int i, int j, int[] columnIndices, double scale, double shape, double[] inverseLength,
       bool fixedInverseLength, bool fixedScale, bool fixedShape) {
       double d = i == j
@@ -151,14 +151,16 @@
       double b = 1 + 0.5 * d / shape;
       int k = 0;
+      var g = new List<double>(columnIndices.Length + 2);
       if (!fixedInverseLength) {
         foreach (var columnIndex in columnIndices) {
-          yield return
+          g.Add(
             scale * Math.Pow(b, -shape - 1) *
-            Util.SqrDist(x[i, columnIndex] * inverseLength[k], x[j, columnIndex] * inverseLength[k]);
+            Util.SqrDist(x[i, columnIndex] * inverseLength[k], x[j, columnIndex] * inverseLength[k]));
           k++;
         }
       }
-      if (!fixedScale) yield return 2 * scale * Math.Pow(b, -shape);
-      if (!fixedShape) yield return scale * Math.Pow(b, -shape) * (0.5 * d / b - shape * Math.Log(b));
+      if (!fixedScale) g.Add(2 * scale * Math.Pow(b, -shape));
+      if (!fixedShape) g.Add(scale * Math.Pow(b, -shape) * (0.5 * d / b - shape * Math.Log(b)));
+      return g;
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceRationalQuadraticIso.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceRationalQuadraticIso.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceRationalQuadraticIso.cs	(revision 13981)
@@ -117,5 +117,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double scale, shape, inverseLength;
       GetParameterValues(p, out scale, out shape, out inverseLength);
@@ -128,9 +128,9 @@
         double d = i == j
                     ? 0.0
-                    : Util.SqrDist(x, i, j, inverseLength, columnIndices);
+                    : Util.SqrDist(x, i, j, columnIndices, inverseLength);
         return scale * Math.Pow(1 + 0.5 * d / shape, -shape);
       };
       cov.CrossCovariance = (x, xt, i, j) => {
-        double d = Util.SqrDist(x, i, xt, j, inverseLength, columnIndices);
+        double d = Util.SqrDist(x, i, xt, j, columnIndices, inverseLength);
         return scale * Math.Pow(1 + 0.5 * d / shape, -shape);
       };
@@ -139,14 +139,16 @@
     }
 
-    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, IEnumerable<int> columnIndices, double scale, double shape, double inverseLength,
+    private static IList<double> GetGradient(double[,] x, int i, int j, int[] columnIndices, double scale, double shape, double inverseLength,
       bool fixedInverseLength, bool fixedScale, bool fixedShape) {
       double d = i == j
                    ? 0.0
-                   : Util.SqrDist(x, i, j, inverseLength, columnIndices);
+                   : Util.SqrDist(x, i, j, columnIndices, inverseLength);
 
       double b = 1 + 0.5 * d / shape;
-      if (!fixedInverseLength) yield return scale * Math.Pow(b, -shape - 1) * d;
-      if (!fixedScale) yield return 2 * scale * Math.Pow(b, -shape);
-      if (!fixedShape) yield return scale * Math.Pow(b, -shape) * (0.5 * d / b - shape * Math.Log(b));
+      var g = new List<double>(3);
+      if (!fixedInverseLength) g.Add(scale * Math.Pow(b, -shape - 1) * d);
+      if (!fixedScale) g.Add(2 * scale * Math.Pow(b, -shape));
+      if (!fixedShape) g.Add(scale * Math.Pow(b, -shape) * (0.5 * d / b - shape * Math.Log(b)));
+      return g;
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceScale.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceScale.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceScale.cs	(revision 13981)
@@ -87,5 +87,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double scale;
       GetParameterValues(p, out scale);
@@ -100,11 +100,13 @@
     }
 
-    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, IEnumerable<int> columnIndices, double scale, ParameterizedCovarianceFunction cov,
+    private static IList<double> GetGradient(double[,] x, int i, int j, int[] columnIndices, double scale, ParameterizedCovarianceFunction cov,
       bool fixedScale) {
+      var gr = new List<double>((!fixedScale ? 1 : 0) + cov.CovarianceGradient(x, i, j).Count);
       if (!fixedScale) {
-        yield return 2 * scale * cov.Covariance(x, i, j);
+        gr.Add(2 * scale * cov.Covariance(x, i, j));
       }
       foreach (var g in cov.CovarianceGradient(x, i, j))
-        yield return scale * g;
+        gr.Add(scale * g);
+      return gr;
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSpectralMixture.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSpectralMixture.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSpectralMixture.cs	(revision 13981)
@@ -23,5 +23,4 @@
 using System.Collections.Generic;
 using System.Linq;
-using System.Linq.Expressions;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -131,5 +130,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double[] weight, frequency, lengthScale;
       GetParameterValues(p, out weight, out frequency, out lengthScale);
@@ -152,5 +151,5 @@
     }
 
-    private static double GetCovariance(double[,] x, double[,] xt, int i, int j, int maxQ, double[] weight, double[] frequency, double[] lengthScale, IEnumerable<int> columnIndices) {
+    private static double GetCovariance(double[,] x, double[,] xt, int i, int j, int maxQ, double[] weight, double[] frequency, double[] lengthScale, int[] columnIndices) {
       // tau = x - x' (only for selected variables)
       double[] tau =
@@ -164,6 +163,7 @@
         int idx = 0; // helper index for tau
         // for each selected variable
-        foreach (var c in columnIndices) {
-          kc *= f1(tau[idx], lengthScale[q * numberOfVariables + c]) * f2(tau[idx], frequency[q * numberOfVariables + c]);
+        for (int c = 0; c < columnIndices.Length; c++) {
+          var col = columnIndices[c];
+          kc *= f1(tau[idx], lengthScale[q * numberOfVariables + col]) * f2(tau[idx], frequency[q * numberOfVariables + col]);
           idx++;
         }
@@ -181,9 +181,10 @@
 
     // order of returned gradients must match the order in GetParameterValues!
-    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, int maxQ, double[] weight, double[] frequency, double[] lengthScale, IEnumerable<int> columnIndices,
+    private static IList<double> GetGradient(double[,] x, int i, int j, int maxQ, double[] weight, double[] frequency, double[] lengthScale, int[] columnIndices,
       bool fixedWeight, bool fixedFrequency, bool fixedLengthScale) {
       double[] tau = Util.GetRow(x, i, columnIndices).Zip(Util.GetRow(x, j, columnIndices), (xi, xj) => xi - xj).ToArray();
       int numberOfVariables = lengthScale.Length / maxQ;
 
+      var g = new List<double>((!fixedWeight ? maxQ : 0) + (!fixedFrequency ? maxQ * columnIndices.Length : 0) + (!fixedLengthScale ? maxQ * columnIndices.Length : 0));
       if (!fixedWeight) {
         // weight
@@ -193,9 +194,10 @@
           int idx = 0; // helper index for tau
           // for each selected variable
-          foreach (var c in columnIndices) {
-            k *= f1(tau[idx], lengthScale[q * numberOfVariables + c]) * f2(tau[idx], frequency[q * numberOfVariables + c]);
+          for (int c = 0; c < columnIndices.Length; c++) {
+            var col = columnIndices[c];
+            k *= f1(tau[idx], lengthScale[q * numberOfVariables + col]) * f2(tau[idx], frequency[q * numberOfVariables + col]);
             idx++;
           }
-          yield return k;
+          g.Add(k);
         }
       }
@@ -212,5 +214,5 @@
                        Math.Sin(2 * Math.PI * tau[idx] * frequency[q * numberOfVariables + c]);
             idx++;
-            yield return weight[q] * k;
+            g.Add(weight[q] * k);
           }
         }
@@ -228,8 +230,10 @@
                        f2(tau[idx], frequency[q * numberOfVariables + c]);
             idx++;
-            yield return weight[q] * k;
+            g.Add(weight[q] * k);
           }
         }
       }
+
+      return g;
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSquaredExponentialArd.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSquaredExponentialArd.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSquaredExponentialArd.cs	(revision 13981)
@@ -99,5 +99,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double scale;
       double[] inverseLength;
@@ -122,5 +122,5 @@
 
     // order of returned gradients must match the order in GetParameterValues!
-    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, IEnumerable<int> columnIndices, double scale, double[] inverseLength,
+    private static IList<double> GetGradient(double[,] x, int i, int j, int[] columnIndices, double scale, double[] inverseLength,
       bool fixedInverseLength, bool fixedScale) {
       double d = i == j
@@ -129,12 +129,15 @@
 
       int k = 0;
+      var g = new List<double>((!fixedInverseLength ? columnIndices.Length : 0) + (!fixedScale ? 1 : 0));
       if (!fixedInverseLength) {
-        foreach (var columnIndex in columnIndices) {
+        for (int c = 0; c < columnIndices.Length; c++) {
+          var columnIndex = columnIndices[c];
           double sqrDist = Util.SqrDist(x[i, columnIndex] * inverseLength[k], x[j, columnIndex] * inverseLength[k]);
-          yield return scale * Math.Exp(-d / 2.0) * sqrDist;
+          g.Add(scale * Math.Exp(-d / 2.0) * sqrDist);
           k++;
         }
       }
-      if (!fixedScale) yield return 2.0 * scale * Math.Exp(-d / 2.0);
+      if (!fixedScale) g.Add(2.0 * scale * Math.Exp(-d / 2.0));
+      return g;
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSquaredExponentialIso.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSquaredExponentialIso.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSquaredExponentialIso.cs	(revision 13981)
@@ -22,5 +22,4 @@
 using System;
 using System.Collections.Generic;
-using System.Linq.Expressions;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -104,5 +103,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       double inverseLength, scale;
       GetParameterValues(p, out scale, out inverseLength);
@@ -114,9 +113,9 @@
         double d = i == j
                 ? 0.0
-                : Util.SqrDist(x, i, j, inverseLength, columnIndices);
+                : Util.SqrDist(x, i, j, columnIndices, inverseLength);
         return scale * Math.Exp(-d / 2.0);
       };
       cov.CrossCovariance = (x, xt, i, j) => {
-        double d = Util.SqrDist(x, i, xt, j, inverseLength, columnIndices);
+        double d = Util.SqrDist(x, i, xt, j, columnIndices, inverseLength);
         return scale * Math.Exp(-d / 2.0);
       };
@@ -127,12 +126,14 @@
 
     // order of returned gradients must match the order in GetParameterValues!
-    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, double sf2, double inverseLength, IEnumerable<int> columnIndices, 
+    private static IList<double> GetGradient(double[,] x, int i, int j, double sf2, double inverseLength, int[] columnIndices,
       bool fixedInverseLength, bool fixedScale) {
       double d = i == j
                    ? 0.0
-                   : Util.SqrDist(x, i, j, inverseLength, columnIndices);
+                   : Util.SqrDist(x, i, j, columnIndices, inverseLength);
       double g = Math.Exp(-d / 2.0);
-      if (!fixedInverseLength) yield return sf2 * g * d;
-      if (!fixedScale) yield return 2.0 * sf2 * g;
+      var gr = new List<double>(2);
+      if (!fixedInverseLength) gr.Add(sf2 * g * d);
+      if (!fixedScale) gr.Add(2.0 * sf2 * g);
+      return gr;
     }
   }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSum.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSum.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSum.cs	(revision 13981)
@@ -23,5 +23,4 @@
 using System.Collections.Generic;
 using System.Linq;
-using System.Linq.Expressions;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -76,5 +75,5 @@
     }
 
-    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
       if (terms.Count == 0) throw new ArgumentException("at least one term is necessary for the product covariance function.");
       var functions = new List<ParameterizedCovarianceFunction>();
@@ -88,5 +87,10 @@
       sum.Covariance = (x, i, j) => functions.Select(e => e.Covariance(x, i, j)).Sum();
       sum.CrossCovariance = (x, xt, i, j) => functions.Select(e => e.CrossCovariance(x, xt, i, j)).Sum();
-      sum.CovarianceGradient = (x, i, j) => functions.Select(e => e.CovarianceGradient(x, i, j)).Aggregate(Enumerable.Concat);
+      sum.CovarianceGradient = (x, i, j) => {
+        var g = new List<double>();
+        foreach (var e in functions)
+          g.AddRange(e.CovarianceGradient(x, i, j));
+        return g;
+      };
       return sum;
     }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs	(revision 13981)
@@ -167,5 +167,6 @@
       try {
         CalculateModel(ds, rows, scaleInputs);
-      } catch (alglib.alglibexception ae) {
+      }
+      catch (alglib.alglibexception ae) {
         // wrap exception so that calling code doesn't have to know about alglib implementation
         throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
@@ -185,10 +186,11 @@
       int n = x.GetLength(0);
 
+      var columns = Enumerable.Range(0, x.GetLength(1)).ToArray();
       // calculate cholesky decomposed (lower triangular) covariance matrix
-      var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
+      var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
       this.l = CalculateL(x, cov, sqrSigmaNoise);
 
       // calculate mean
-      var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, x.GetLength(1)));
+      var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, columns);
       double[] m = Enumerable.Range(0, x.GetLength(0))
         .Select(r => mean.Mean(x, r))
@@ -227,6 +229,7 @@
       double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
       for (int k = 0; k < meanGradients.Length; k++) {
-        var meanGrad = Enumerable.Range(0, alpha.Length)
-        .Select(r => mean.Gradient(x, r, k));
+        var meanGrad = new double[alpha.Length];
+        for (int g = 0; g < meanGrad.Length; g++)
+          meanGrad[g] = mean.Gradient(x, g, k);
         meanGradients[k] = -Util.ScalarProd(meanGrad, alpha);
       }
@@ -236,5 +239,5 @@
         for (int i = 0; i < n; i++) {
           for (int j = 0; j < i; j++) {
-            var g = cov.CovarianceGradient(x, i, j).ToArray();
+            var g = cov.CovarianceGradient(x, i, j);
             for (int k = 0; k < covGradients.Length; k++) {
               covGradients[k] += lCopy[i, j] * g[k];
@@ -242,5 +245,5 @@
           }
 
-          var gDiag = cov.CovarianceGradient(x, i, i).ToArray();
+          var gDiag = cov.CovarianceGradient(x, i, i);
           for (int k = 0; k < covGradients.Length; k++) {
             // diag
@@ -320,19 +323,22 @@
         int newN = newX.GetLength(0);
 
-        var Ks = new double[newN, n];
-        var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, newX.GetLength(1)));
+        var Ks = new double[newN][];
+        var columns = Enumerable.Range(0, newX.GetLength(1)).ToArray();
+        var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, columns);
         var ms = Enumerable.Range(0, newX.GetLength(0))
         .Select(r => mean.Mean(newX, r))
         .ToArray();
-        var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, newX.GetLength(1)));
+        var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
         for (int i = 0; i < newN; i++) {
+          Ks[i] = new double[n];
           for (int j = 0; j < n; j++) {
-            Ks[i, j] = cov.CrossCovariance(x, newX, j, i);
+            Ks[i][j] = cov.CrossCovariance(x, newX, j, i);
           }
         }
 
         return Enumerable.Range(0, newN)
-          .Select(i => ms[i] + Util.ScalarProd(Util.GetRow(Ks, i), alpha));
-      } catch (alglib.alglibexception ae) {
+          .Select(i => ms[i] + Util.ScalarProd(Ks[i], alpha));
+      }
+      catch (alglib.alglibexception ae) {
         // wrap exception so that calling code doesn't have to know about alglib implementation
         throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
@@ -352,5 +358,6 @@
         var kss = new double[newN];
         double[,] sWKs = new double[n, newN];
-        var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
+        var columns = Enumerable.Range(0, newX.GetLength(1)).ToArray();
+        var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
 
         if (l == null) {
@@ -372,5 +379,6 @@
 
         for (int i = 0; i < newN; i++) {
-          var sumV = Util.ScalarProd(Util.GetCol(sWKs, i), Util.GetCol(sWKs, i));
+          var col = Util.GetCol(sWKs, i).ToArray();
+          var sumV = Util.ScalarProd(col, col);
           kss[i] += sqrSigmaNoise; // kss is V(f), add noise variance of predictive distibution to get V(y)
           kss[i] -= sumV;
@@ -378,5 +386,6 @@
         }
         return kss;
-      } catch (alglib.alglibexception ae) {
+      }
+      catch (alglib.alglibexception ae) {
         // wrap exception so that calling code doesn't have to know about alglib implementation
         throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/ICovarianceFunction.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/ICovarianceFunction.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/ICovarianceFunction.cs	(revision 13981)
@@ -20,7 +20,5 @@
 #endregion
 
-using System;
 using System.Collections.Generic;
-using System.Linq.Expressions;
 using HeuristicLab.Core;
 
@@ -29,5 +27,5 @@
   public delegate double CovarianceFunctionDelegate(double[,] x, int i, int j);
   public delegate double CrossCovarianceFunctionDelegate(double[,] x, double[,] xt, int i, int j);
-  public delegate IEnumerable<double> CovarianceGradientFunctionDelegate(double[,] x, int i, int j);
+  public delegate IList<double> CovarianceGradientFunctionDelegate(double[,] x, int i, int j);
 
   public class ParameterizedCovarianceFunction {
@@ -40,5 +38,5 @@
     int GetNumberOfParameters(int numberOfVariables);
     void SetParameter(double[] p);
-    ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices);
+    ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices);
   }
 }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/IMeanFunction.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/IMeanFunction.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/IMeanFunction.cs	(revision 13981)
@@ -20,6 +20,4 @@
 #endregion
 
-using System;
-using System.Collections.Generic;
 using HeuristicLab.Core;
 
@@ -36,5 +34,5 @@
     int GetNumberOfParameters(int numberOfVariables);
     void SetParameter(double[] p);
-    ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, IEnumerable<int> columnIndices);
+    ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, int[] columnIndices);
   }
 }
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanConst.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanConst.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanConst.cs	(revision 13981)
@@ -76,5 +76,5 @@
     }
 
-    public ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, int[] columnIndices) {
       double c;
       GetParameters(p, out c);
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanLinear.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanLinear.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanLinear.cs	(revision 13981)
@@ -21,5 +21,4 @@
 
 using System;
-using System.Collections.Generic;
 using System.Linq;
 using HeuristicLab.Common;
@@ -70,7 +69,7 @@
     }
 
-    public ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, int[] columnIndices) {
       double[] weights;
-      int[] columns = columnIndices.ToArray();
+      int[] columns = columnIndices;
       GetParameter(p, out weights);
       var mf = new ParameterizedMeanFunction();
@@ -78,5 +77,5 @@
         // sanity check
         if (weights.Length != columns.Length) throw new ArgumentException("The number of rparameters must match the number of variables for the linear mean function.");
-        return Util.ScalarProd(weights, Util.GetRow(x, i, columns));
+        return Util.ScalarProd(weights, Util.GetRow(x, i, columns).ToArray());
       };
       mf.Gradient = (x, i, k) => {
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanModel.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanModel.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanModel.cs	(revision 13981)
@@ -73,9 +73,9 @@
     }
 
-    public ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, int[] columnIndices) {
       if (p.Length > 0) throw new ArgumentException("No parameters allowed for model-based mean function.", "p");
       var solution = RegressionSolution;
       var variableNames = solution.ProblemData.AllowedInputVariables.ToArray();
-      if (variableNames.Length != columnIndices.Count())
+      if (variableNames.Length != columnIndices.Length)
         throw new ArgumentException("The number of input variables does not match in MeanModel");
       var variableValues = variableNames.Select(_ => new List<double>() { 0.0 }).ToArray(); // or of zeros
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanProduct.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanProduct.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanProduct.cs	(revision 13981)
@@ -73,5 +73,5 @@
 
 
-    public ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, int[] columnIndices) {
       var factorMf = new List<ParameterizedMeanFunction>();
       int totalNumberOfParameters = GetNumberOfParameters(numberOfVariables);
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanSum.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanSum.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanSum.cs	(revision 13981)
@@ -68,5 +68,5 @@
     }
 
-    public ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, int[] columnIndices) {
       var termMf = new List<ParameterizedMeanFunction>();
       int totalNumberOfParameters = GetNumberOfParameters(numberOfVariables);
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanZero.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanZero.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanFunctions/MeanZero.cs	(revision 13981)
@@ -20,6 +20,4 @@
 #endregion
 using System;
-using System.Collections.Generic;
-using System.Linq;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -50,5 +48,5 @@
     }
 
-    public ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, IEnumerable<int> columnIndices) {
+    public ParameterizedMeanFunction GetParameterizedMeanFunction(double[] p, int[] columnIndices) {
       if (p.Length > 0) throw new ArgumentException("No parameters allowed for zero mean function.", "p");
       var mf = new ParameterizedMeanFunction();
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/StudentTProcessModel.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/StudentTProcessModel.cs	(revision 13981)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/StudentTProcessModel.cs	(revision 13981)
@@ -0,0 +1,411 @@
+#region License Information
+/* HeuristicLab
+ * Copyright (C) 2002-2015 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.Generic;
+using System.Linq;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+using HeuristicLab.Problems.DataAnalysis;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  /// <summary>
+  /// Represents a Gaussian process model.
+  /// </summary>
+  [StorableClass]
+  [Item("StudentTProcessModel", "Represents a Student-t process posterior.")]
+  public sealed class StudentTProcessModel : NamedItem, IGaussianProcessModel {
+    [Storable]
+    private double negativeLogLikelihood;
+    public double NegativeLogLikelihood {
+      get { return negativeLogLikelihood; }
+    }
+
+    [Storable]
+    private double[] hyperparameterGradients;
+    public double[] HyperparameterGradients {
+      get {
+        var copy = new double[hyperparameterGradients.Length];
+        Array.Copy(hyperparameterGradients, copy, copy.Length);
+        return copy;
+      }
+    }
+
+    [Storable]
+    private ICovarianceFunction covarianceFunction;
+    public ICovarianceFunction CovarianceFunction {
+      get { return covarianceFunction; }
+    }
+    [Storable]
+    private IMeanFunction meanFunction;
+    public IMeanFunction MeanFunction {
+      get { return meanFunction; }
+    }
+    [Storable]
+    private string targetVariable;
+    public string TargetVariable {
+      get { return targetVariable; }
+    }
+    [Storable]
+    private string[] allowedInputVariables;
+    public string[] AllowedInputVariables {
+      get { return allowedInputVariables; }
+    }
+
+    [Storable]
+    private double[] alpha;
+    [Storable]
+    private double beta;
+
+    public double SigmaNoise {
+      get { return 0; }
+    }
+
+    [Storable]
+    private double[] meanParameter;
+    [Storable]
+    private double[] covarianceParameter;
+    [Storable]
+    private double nu;
+
+    private double[,] l; // used to be storable in previous versions (is calculated lazily now)
+    private double[,] x; // scaled training dataset, used to be storable in previous versions (is calculated lazily now)
+
+    // BackwardsCompatibility3.4
+    #region Backwards compatible code, remove with 3.5
+    [Storable(Name = "l")] // restore if available but don't store anymore
+    private double[,] l_storable {
+      set { this.l = value; }
+      get {
+        if (trainingDataset == null) return l; // this model has been created with an old version 
+        else return null; // if the training dataset is available l should not be serialized
+      }
+    }
+    [Storable(Name = "x")] // restore if available but don't store anymore
+    private double[,] x_storable {
+      set { this.x = value; }
+      get {
+        if (trainingDataset == null) return x; // this model has been created with an old version 
+        else return null; // if the training dataset is available x should not be serialized
+      }
+    }
+    #endregion
+
+
+    [Storable]
+    private IDataset trainingDataset; // it is better to store the original training dataset completely because this is more efficient in persistence
+    [Storable]
+    private int[] trainingRows;
+
+    [Storable]
+    private Scaling inputScaling;
+
+
+    [StorableConstructor]
+    private StudentTProcessModel(bool deserializing) : base(deserializing) { }
+    private StudentTProcessModel(StudentTProcessModel original, Cloner cloner)
+      : base(original, cloner) {
+      this.meanFunction = cloner.Clone(original.meanFunction);
+      this.covarianceFunction = cloner.Clone(original.covarianceFunction);
+      if (original.inputScaling != null)
+        this.inputScaling = cloner.Clone(original.inputScaling);
+      this.trainingDataset = cloner.Clone(original.trainingDataset);
+      this.negativeLogLikelihood = original.negativeLogLikelihood;
+      this.targetVariable = original.targetVariable;
+      if (original.meanParameter != null) {
+        this.meanParameter = (double[])original.meanParameter.Clone();
+      }
+      if (original.covarianceParameter != null) {
+        this.covarianceParameter = (double[])original.covarianceParameter.Clone();
+      }
+      nu = original.nu;
+
+      // shallow copies of arrays because they cannot be modified
+      this.trainingRows = original.trainingRows;
+      this.allowedInputVariables = original.allowedInputVariables;
+      this.alpha = original.alpha;
+      this.beta = original.beta;
+      this.l = original.l;
+      this.x = original.x;
+    }
+    public StudentTProcessModel(IDataset ds, string targetVariable, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows,
+      IEnumerable<double> hyp, IMeanFunction meanFunction, ICovarianceFunction covarianceFunction,
+      bool scaleInputs = true)
+      : base() {
+      this.name = ItemName;
+      this.description = ItemDescription;
+      this.meanFunction = (IMeanFunction)meanFunction.Clone();
+      this.covarianceFunction = (ICovarianceFunction)covarianceFunction.Clone();
+      this.targetVariable = targetVariable;
+      this.allowedInputVariables = allowedInputVariables.ToArray();
+
+
+      int nVariables = this.allowedInputVariables.Length;
+      meanParameter = hyp
+        .Take(this.meanFunction.GetNumberOfParameters(nVariables))
+        .ToArray();
+
+      covarianceParameter = hyp.Skip(meanParameter.Length)
+                                             .Take(this.covarianceFunction.GetNumberOfParameters(nVariables))
+                                             .ToArray();
+      nu = Math.Exp(hyp.Skip(meanParameter.Length + covarianceParameter.Length).First()) + 2; //TODO check gradient
+      try {
+        CalculateModel(ds, rows, scaleInputs);
+      }
+      catch (alglib.alglibexception ae) {
+        // wrap exception so that calling code doesn't have to know about alglib implementation
+        throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
+      }
+    }
+
+    private void CalculateModel(IDataset ds, IEnumerable<int> rows, bool scaleInputs = true) {
+      this.trainingDataset = (IDataset)ds.Clone();
+      this.trainingRows = rows.ToArray();
+      this.inputScaling = scaleInputs ? new Scaling(ds, allowedInputVariables, rows) : null;
+
+      x = GetData(ds, this.allowedInputVariables, this.trainingRows, this.inputScaling);
+
+      IEnumerable<double> y;
+      y = ds.GetDoubleValues(targetVariable, rows);
+
+      int n = x.GetLength(0);
+      var columns = Enumerable.Range(0, x.GetLength(1)).ToArray();
+
+      // calculate cholesky decomposed (lower triangular) covariance matrix
+      var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
+      this.l = CalculateL(x, cov);
+
+      // calculate mean
+      var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, columns);
+      double[] m = Enumerable.Range(0, x.GetLength(0))
+        .Select(r => mean.Mean(x, r))
+        .ToArray();
+
+      // calculate sum of diagonal elements for likelihood
+      double diagSum = Enumerable.Range(0, n).Select(i => Math.Log(l[i, i])).Sum();
+
+      // solve for alpha
+      double[] ym = y.Zip(m, (a, b) => a - b).ToArray();
+
+      int info;
+      alglib.densesolverreport denseSolveRep;
+
+      alglib.spdmatrixcholeskysolve(l, n, false, ym, out info, out denseSolveRep, out alpha);
+
+      beta = Util.ScalarProd(ym, alpha);
+      double sign0, sign1;
+      double lngamma0 = alglib.lngamma(0.5 * (nu + n), out sign0);
+      lngamma0 *= sign0;
+      double lngamma1 = alglib.lngamma(0.5 * nu, out sign1);
+      lngamma1 *= sign1;
+      negativeLogLikelihood =
+        0.5 * n * Math.Log((nu - 2) * Math.PI) +
+        diagSum +
+        -lngamma0 + lngamma1 +
+        //-Math.Log(alglib.gammafunction((n + nu) / 2) / alglib.gammafunction(nu / 2)) +
+        0.5 * (nu + n) * Math.Log(1 + beta / (nu - 2));
+
+      // derivatives
+      int nAllowedVariables = x.GetLength(1);
+
+      alglib.matinvreport matInvRep;
+      double[,] lCopy = new double[l.GetLength(0), l.GetLength(1)];
+      Array.Copy(l, lCopy, lCopy.Length);
+
+      alglib.spdmatrixcholeskyinverse(ref lCopy, n, false, out info, out matInvRep);
+      double c = (nu + n) / (nu + beta - 2);
+      if (info != 1) throw new ArgumentException("Can't invert matrix to calculate gradients.");
+      for (int i = 0; i < n; i++) {
+        for (int j = 0; j <= i; j++)
+          lCopy[i, j] = lCopy[i, j] - c * alpha[i] * alpha[j];
+      }
+
+      double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
+      for (int k = 0; k < meanGradients.Length; k++) {
+        var meanGrad = new double[alpha.Length];
+        for (int g = 0; g < meanGrad.Length; g++)
+          meanGrad[g] = mean.Gradient(x, g, k);
+        meanGradients[k] = -Util.ScalarProd(meanGrad, alpha);//TODO not working yet, try to fix with gradient check
+      }
+
+      double[] covGradients = new double[covarianceFunction.GetNumberOfParameters(nAllowedVariables)];
+      if (covGradients.Length > 0) {
+        for (int i = 0; i < n; i++) {
+          for (int j = 0; j < i; j++) {
+            var g = cov.CovarianceGradient(x, i, j);
+            for (int k = 0; k < covGradients.Length; k++) {
+              covGradients[k] += lCopy[i, j] * g[k];
+            }
+          }
+
+          var gDiag = cov.CovarianceGradient(x, i, i);
+          for (int k = 0; k < covGradients.Length; k++) {
+            // diag
+            covGradients[k] += 0.5 * lCopy[i, i] * gDiag[k];
+          }
+        }
+      }
+
+      double nuGradient = 0.5 * n
+        - 0.5 * (nu - 2) * alglib.psi((n + nu) / 2) + 0.5 * (nu - 2) * alglib.psi(nu / 2)
+        + 0.5 * (nu - 2) * Math.Log(1 + beta / (nu - 2)) - beta * (n + nu) / (2 * (beta + (nu - 2)));
+
+      //nuGradient = (nu-2) * nuGradient;
+      hyperparameterGradients =
+        meanGradients
+        .Concat(covGradients)
+        .Concat(new double[] { nuGradient }).ToArray();
+
+    }
+
+    private static double[,] GetData(IDataset ds, IEnumerable<string> allowedInputs, IEnumerable<int> rows, Scaling scaling) {
+      if (scaling != null) {
+        return AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputs, rows, scaling);
+      } else {
+        return AlglibUtil.PrepareInputMatrix(ds, allowedInputs, rows);
+      }
+    }
+
+    private static double[,] CalculateL(double[,] x, ParameterizedCovarianceFunction cov) {
+      int n = x.GetLength(0);
+      var l = new double[n, n];
+
+      // calculate covariances
+      for (int i = 0; i < n; i++) {
+        for (int j = i; j < n; j++) {
+          l[j, i] = cov.Covariance(x, i, j);
+        }
+      }
+
+      // cholesky decomposition
+      var res = alglib.trfac.spdmatrixcholesky(ref l, n, false);
+      if (!res) throw new ArgumentException("Matrix is not positive semidefinite");
+      return l;
+    }
+
+
+    public override IDeepCloneable Clone(Cloner cloner) {
+      return new StudentTProcessModel(this, cloner);
+    }
+
+    // is called by the solution creator to set all parameter values of the covariance and mean function
+    // to the optimized values (necessary to make the values visible in the GUI)
+    public void FixParameters() {
+      covarianceFunction.SetParameter(covarianceParameter);
+      meanFunction.SetParameter(meanParameter);
+      covarianceParameter = new double[0];
+      meanParameter = new double[0];
+    }
+
+    #region IRegressionModel Members
+    public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
+      return GetEstimatedValuesHelper(dataset, rows);
+    }
+    public GaussianProcessRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
+      return new GaussianProcessRegressionSolution(this, new RegressionProblemData(problemData));
+    }
+    IRegressionSolution IRegressionModel.CreateRegressionSolution(IRegressionProblemData problemData) {
+      return CreateRegressionSolution(problemData);
+    }
+    #endregion
+
+
+    private IEnumerable<double> GetEstimatedValuesHelper(IDataset dataset, IEnumerable<int> rows) {
+      try {
+        if (x == null) {
+          x = GetData(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
+        }
+        int n = x.GetLength(0);
+
+        double[,] newX = GetData(dataset, allowedInputVariables, rows, inputScaling);
+        int newN = newX.GetLength(0);
+        var columns = Enumerable.Range(0, newX.GetLength(1)).ToArray();
+
+        var Ks = new double[newN][];
+        var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, columns);
+        var ms = Enumerable.Range(0, newX.GetLength(0))
+        .Select(r => mean.Mean(newX, r))
+        .ToArray();
+        var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
+        for (int i = 0; i < newN; i++) {
+          Ks[i] = new double[n];
+          for (int j = 0; j < n; j++) {
+            Ks[i][j] = cov.CrossCovariance(x, newX, j, i);
+          }
+        }
+
+        return Enumerable.Range(0, newN)
+          .Select(i => ms[i] + Util.ScalarProd(Ks[i], alpha));
+      }
+      catch (alglib.alglibexception ae) {
+        // wrap exception so that calling code doesn't have to know about alglib implementation
+        throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
+      }
+    }
+
+    public IEnumerable<double> GetEstimatedVariance(IDataset dataset, IEnumerable<int> rows) {
+      try {
+        if (x == null) {
+          x = GetData(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
+        }
+        int n = x.GetLength(0);
+
+        var newX = GetData(dataset, allowedInputVariables, rows, inputScaling);
+        int newN = newX.GetLength(0);
+
+        var kss = new double[newN];
+        double[,] sWKs = new double[n, newN];
+        var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)).ToArray());
+
+        if (l == null) {
+          l = CalculateL(x, cov);
+        }
+
+        // for stddev 
+        for (int i = 0; i < newN; i++)
+          kss[i] = cov.Covariance(newX, i, i);
+
+        for (int i = 0; i < newN; i++) {
+          for (int j = 0; j < n; j++) {
+            sWKs[j, i] = cov.CrossCovariance(x, newX, j, i);
+          }
+        }
+
+        // for stddev 
+        alglib.ablas.rmatrixlefttrsm(n, newN, l, 0, 0, false, false, 0, ref sWKs, 0, 0);
+
+        for (int i = 0; i < newN; i++) {
+          var col = Util.GetCol(sWKs, i).ToArray();
+          var sumV = Util.ScalarProd(col, col);
+          kss[i] -= sumV;
+          kss[i] *= (nu + beta - 2) / (nu + n - 2);
+          if (kss[i] < 0) kss[i] = 0;
+        }
+        return kss;
+      }
+      catch (alglib.alglibexception ae) {
+        // wrap exception so that calling code doesn't have to know about alglib implementation
+        throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
+      }
+    }
+  }
+}
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/StudentTProcessRegressionModelCreator.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/StudentTProcessRegressionModelCreator.cs	(revision 13981)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/StudentTProcessRegressionModelCreator.cs	(revision 13981)
@@ -0,0 +1,81 @@
+﻿#region License Information
+/* HeuristicLab
+ * Copyright (C) 2002-2015 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.Linq;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Data;
+using HeuristicLab.Encodings.RealVectorEncoding;
+using HeuristicLab.Parameters;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+using HeuristicLab.Problems.DataAnalysis;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  [StorableClass]
+  [Item(Name = "StudentTProcessRegressionModelCreator",
+    Description = "Creates a Student-t process model for regression given the data, the hyperparameters, a mean function, and a covariance function.")]
+  public sealed class StudentTProcessRegressionModelCreator : GaussianProcessModelCreator, IGaussianProcessRegressionModelCreator {
+    private const string ProblemDataParameterName = "ProblemData";
+
+    #region Parameter Properties
+    public ILookupParameter<IRegressionProblemData> ProblemDataParameter {
+      get { return (ILookupParameter<IRegressionProblemData>)Parameters[ProblemDataParameterName]; }
+    }
+    #endregion
+
+    #region Properties
+    private IRegressionProblemData ProblemData {
+      get { return ProblemDataParameter.ActualValue; }
+    }
+    #endregion
+    [StorableConstructor]
+    private StudentTProcessRegressionModelCreator(bool deserializing) : base(deserializing) { }
+    private StudentTProcessRegressionModelCreator(StudentTProcessRegressionModelCreator original, Cloner cloner) : base(original, cloner) { }
+    public StudentTProcessRegressionModelCreator()
+      : base() {
+      Parameters.Add(new LookupParameter<IRegressionProblemData>(ProblemDataParameterName, "The regression problem data for the Gaussian process model."));
+    }
+
+    public override IDeepCloneable Clone(Cloner cloner) {
+      return new StudentTProcessRegressionModelCreator(this, cloner);
+    }
+
+    public override IOperation Apply() {
+      try {
+        var model = Create(ProblemData, Hyperparameter.ToArray(), MeanFunction, CovarianceFunction, ScaleInputValues);
+        ModelParameter.ActualValue = model;
+        NegativeLogLikelihoodParameter.ActualValue = new DoubleValue(model.NegativeLogLikelihood);
+        HyperparameterGradientsParameter.ActualValue = new RealVector(model.HyperparameterGradients);
+        return base.Apply();
+      }
+      catch (ArgumentException) { }
+      catch (alglib.alglibexception) { }
+      NegativeLogLikelihoodParameter.ActualValue = new DoubleValue(1E300);
+      HyperparameterGradientsParameter.ActualValue = new RealVector(Hyperparameter.Count());
+      return base.Apply();
+    }
+
+    public static IGaussianProcessModel Create(IRegressionProblemData problemData, double[] hyperparameter, IMeanFunction meanFunction, ICovarianceFunction covarianceFunction, bool scaleInputs = true) {
+      return new StudentTProcessModel(problemData.Dataset, problemData.TargetVariable, problemData.AllowedInputVariables, problemData.TrainingIndices, hyperparameter, meanFunction, covarianceFunction, scaleInputs);
+    }
+  }
+}
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/Util.cs
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/Util.cs	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/Util.cs	(revision 13981)
@@ -23,11 +23,13 @@
 using System.Collections.Generic;
 using System.Linq;
-using HeuristicLab.Core;
-using HeuristicLab.Data;
 
 namespace HeuristicLab.Algorithms.DataAnalysis {
   internal static class Util {
-    public static double ScalarProd(IEnumerable<double> v, IEnumerable<double> u) {
-      return v.Zip(u, (vi, ui) => vi * ui).Sum();
+    public static double ScalarProd(double[] v, double[] u) {
+      if (v.Length != u.Length) throw new InvalidOperationException();
+      double prod = 0.0;
+      for (int i = 0; i < v.Length; i++)
+        prod += v[i] * u[i];
+      return prod;
     }
 
@@ -41,12 +43,13 @@
     }
 
-    public static double SqrDist(double[,] x, int i, int j, double scale = 1.0, IEnumerable<int> columnIndices = null) {
-      return SqrDist(x, i, x, j, scale, columnIndices);
+    public static double SqrDist(double[,] x, int i, int j, int[] columnIndices, double scale = 1.0) {
+      return SqrDist(x, i, x, j, columnIndices, scale);
     }
 
-    public static double SqrDist(double[,] x, int i, double[,] xt, int j, double scale = 1.0, IEnumerable<int> columnIndices = null) {
+    public static double SqrDist(double[,] x, int i, double[,] xt, int j, int[] columnIndices, double scale = 1.0) {
       double ss = 0.0;
-      if (columnIndices == null) columnIndices = Enumerable.Range(0, x.GetLength(1));
-      foreach (int columnIndex in columnIndices) {
+      if (columnIndices == null) columnIndices = Enumerable.Range(0, x.GetLength(1)).ToArray();
+      for (int c = 0; c < columnIndices.Length; c++) {
+        var columnIndex = columnIndices[c];
         double d = x[i, columnIndex] - xt[j, columnIndex];
         ss += d * d;
@@ -55,13 +58,13 @@
     }
 
-    public static double SqrDist(double[,] x, int i, int j, double[] scale, IEnumerable<int> columnIndices = null) {
+    public static double SqrDist(double[,] x, int i, int j, double[] scale, int[] columnIndices) {
       return SqrDist(x, i, x, j, scale, columnIndices);
     }
 
-    public static double SqrDist(double[,] x, int i, double[,] xt, int j, double[] scale, IEnumerable<int> columnIndices = null) {
+    public static double SqrDist(double[,] x, int i, double[,] xt, int j, double[] scale, int[] columnIndices) {
       double ss = 0.0;
-      if (columnIndices == null) columnIndices = Enumerable.Range(0, x.GetLength(1));
       int scaleIndex = 0;
-      foreach (int columnIndex in columnIndices) {
+      for (int c = 0; c < columnIndices.Length; c++) {
+        var columnIndex = columnIndices[c];
         double d = x[i, columnIndex] - xt[j, columnIndex];
         ss += d * d * scale[scaleIndex] * scale[scaleIndex];
@@ -73,27 +76,26 @@
       return ss;
     }
-    public static double ScalarProd(double[,] x, int i, int j, double scale = 1.0, IEnumerable<int> columnIndices = null) {
-      return ScalarProd(x, i, x, j, scale, columnIndices);
+    public static double ScalarProd(double[,] x, int i, int j, int[] columnIndices, double scale = 1.0) {
+      return ScalarProd(x, i, x, j, columnIndices, scale);
     }
 
-    public static double ScalarProd(double[,] x, int i, double[,] xt, int j, double scale = 1.0, IEnumerable<int> columnIndices = null) {
+    public static double ScalarProd(double[,] x, int i, double[,] xt, int j, int[] columnIndices, double scale = 1.0) {
       double sum = 0.0;
-      if (columnIndices == null) columnIndices = Enumerable.Range(0, x.GetLength(1));
-      foreach (int columnIndex in columnIndices) {
+      for (int c = 0; c < columnIndices.Length; c++) {
+        var columnIndex = columnIndices[c];
         sum += x[i, columnIndex] * xt[j, columnIndex];
       }
       return scale * scale * sum;
     }
-    public static double ScalarProd(double[,] x, int i, int j, double[] scale, IEnumerable<int> columnIndices = null) {
+    public static double ScalarProd(double[,] x, int i, int j, double[] scale, int[] columnIndices) {
       return ScalarProd(x, i, x, j, scale, columnIndices);
     }
 
-    public static double ScalarProd(double[,] x, int i, double[,] xt, int j, double[] scale, IEnumerable<int> columnIndices = null) {
+    public static double ScalarProd(double[,] x, int i, double[,] xt, int j, double[] scale, int[] columnIndices) {
       double sum = 0.0;
-      if (columnIndices == null) columnIndices = Enumerable.Range(0, x.GetLength(1));
       int scaleIndex = 0;
-      foreach (int columnIndex in columnIndices) {
+      for (int c = 0; c < columnIndices.Length; c++, scaleIndex++) {
+        var columnIndex = columnIndices[c];
         sum += x[i, columnIndex] * scale[scaleIndex] * xt[j, columnIndex] * scale[scaleIndex];
-        scaleIndex++;
       }
       // must be at the end of scale after iterating over columnIndices
Index: /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/HeuristicLab.Algorithms.DataAnalysis-3.4.csproj
===================================================================
--- /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/HeuristicLab.Algorithms.DataAnalysis-3.4.csproj	(revision 13980)
+++ /stable/HeuristicLab.Algorithms.DataAnalysis/3.4/HeuristicLab.Algorithms.DataAnalysis-3.4.csproj	(revision 13981)
@@ -130,4 +130,6 @@
     </Compile>
     <Compile Include="FixedDataAnalysisAlgorithm.cs" />
+    <Compile Include="GaussianProcess\StudentTProcessRegressionModelCreator.cs" />
+    <Compile Include="GaussianProcess\StudentTProcessModel.cs" />
     <Compile Include="GaussianProcess\CovarianceFunctions\CovarianceSpectralMixture.cs" />
     <Compile Include="GaussianProcess\CovarianceFunctions\CovariancePiecewisePolynomial.cs" />
Index: /stable/HeuristicLab.Tests/HeuristicLab.Algorithms.DataAnalysis-3.4/GaussianProcessFunctionsTest.cs
===================================================================
--- /stable/HeuristicLab.Tests/HeuristicLab.Algorithms.DataAnalysis-3.4/GaussianProcessFunctionsTest.cs	(revision 13980)
+++ /stable/HeuristicLab.Tests/HeuristicLab.Algorithms.DataAnalysis-3.4/GaussianProcessFunctionsTest.cs	(revision 13981)
@@ -1747,5 +1747,5 @@
       int rows1 = xt.GetLength(0);
       var actualCov = new double[rows0, rows1];
-      var covFunction = cf.GetParameterizedCovarianceFunction(hyp, Enumerable.Range(0, x.GetLength(1)));
+      var covFunction = cf.GetParameterizedCovarianceFunction(hyp, Enumerable.Range(0, x.GetLength(1)).ToArray());
       for (int i = 0; i < rows0; i++)
         for (int j = 0; j < rows1; j++)
@@ -1769,5 +1769,5 @@
       int nHyp = mf.GetNumberOfParameters(x.GetLength(1));
       var hyp = Enumerable.Repeat(hypValue, nHyp).ToArray();
-      var meanFunction = mf.GetParameterizedMeanFunction(hyp, Enumerable.Range(0, x.GetLength(1)));
+      var meanFunction = mf.GetParameterizedMeanFunction(hyp, Enumerable.Range(0, x.GetLength(1)).ToArray());
 
       var m = Enumerable.Range(0, xt.GetLength(0)).Select(i => meanFunction.Mean(xt, i)).ToArray();
Index: /stable/HeuristicLab.Tests/HeuristicLab.Algorithms.DataAnalysis-3.4/GaussianProcessRegressionTest.cs
===================================================================
--- /stable/HeuristicLab.Tests/HeuristicLab.Algorithms.DataAnalysis-3.4/GaussianProcessRegressionTest.cs	(revision 13980)
+++ /stable/HeuristicLab.Tests/HeuristicLab.Algorithms.DataAnalysis-3.4/GaussianProcessRegressionTest.cs	(revision 13981)
@@ -50,4 +50,5 @@
       var alg = new GaussianProcessRegression();
       alg.Engine = new HeuristicLab.SequentialEngine.SequentialEngine();
+      alg.SetSeedRandomly = false;
 
       alg.Problem = new RegressionProblem();
