Changeset 8475 for trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess

Timestamp:

08/11/12 14:45:15 (12 years ago)

Author:

gkronber

Message:

#1902 fixed bug in calculation of variance in GPR model

File:

• trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs (modified) (6 diffs)

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs

@@ -130 +130 @@
       double[] m = meanFunction.GetMean(x);
       for (int i = 0; i < n; i++) {
-
         for (int j = i; j < n; j++) {
           l[j, i] = covarianceFunction.GetCovariance(i, j) / sqrSigmaNoise;
@@ -163 +162 @@
       int info;
       alglib.matinvreport matInvRep;
-
-      alglib.spdmatrixcholeskyinverse(ref l, n, false, out info, out 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);
       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++)
-          l[i, j] = l[i, j] / sqrSigmaNoise - alpha[i] * alpha[j];
-      }
-
-      double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => l[i, i]).Sum();
+          lCopy[i, j] = lCopy[i, j] / sqrSigmaNoise - alpha[i] * alpha[j];
+      }
+
+      double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => lCopy[i, i]).Sum();
 
       double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
@@ -184 +185 @@
           for (int k = 0; k < covGradients.Length; k++) {
             for (int j = 0; j < i; j++) {
-              covGradients[k] += l[i, j] * covarianceFunction.GetGradient(i, j, k);
+              covGradients[k] += lCopy[i, j] * covarianceFunction.GetGradient(i, j, k);
             }
-            covGradients[k] += 0.5 * l[i, i] * covarianceFunction.GetGradient(i, i, k);
+            covGradients[k] += 0.5 * lCopy[i, i] * covarianceFunction.GetGradient(i, i, k);
           }
         }
@@ -264 +265 @@
       double[,] sWKs = new double[n, newN];
 
-
       // for stddev 
       covarianceFunction.SetData(newX);
@@ -271 +271 @@
 
       covarianceFunction.SetData(x, newX);
-      for (int i = 0; i < n; i++) {
-        for (int j = 0; j < newN; j++) {
-          sWKs[i, j] = covarianceFunction.GetCovariance(i, j) / Math.Sqrt(sqrSigmaNoise);
+      for (int i = 0; i < newN; i++) {
+        for (int j = 0; j < n; j++) {
+          sWKs[j, i] = covarianceFunction.GetCovariance(j, i) / Math.Sqrt(sqrSigmaNoise);
         }
       }
@@ -281 +281 @@
       alglib.densesolverreport denseSolveRep;
       double[,] v;
-      double[,] lTrans = new double[l.GetLength(1), l.GetLength(0)];
-      for (int i = 0; i < lTrans.GetLength(0); i++)
-        for (int j = 0; j < lTrans.GetLength(1); j++)
-          lTrans[i, j] = l[j, i];
-      alglib.rmatrixsolvem(lTrans, n, sWKs, newN, true, out info, out denseSolveRep, out v); // not working!
-      // alglib.spdmatrixcholeskysolvem(lTrans, n, true, sWKs, newN, out info, out denseSolveRep, out v);
+
+      alglib.rmatrixsolvem(l, n, sWKs, newN, false, out info, out denseSolveRep, out v);
 
       for (int i = 0; i < newN; i++) {
-        var sumV2 = Util.ScalarProd(Util.GetCol(v, i), Util.GetCol(v, i));
-        yield return kss[i] - sumV2;
-      }
+        var sumV = Util.ScalarProd(Util.GetCol(v, i), Util.GetCol(v, i));
+        kss[i] -= sumV;
+        if (kss[i] < 0) kss[i] = 0;
+      }
+      return kss;
     }
   }