Changeset 8417 for trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovariancePeriodic.cs

Timestamp:

08/06/12 16:16:28 (12 years ago)

Author:

gkronber

Message:

#1902 added periodic covariance function

File:

• trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovariancePeriodic.cs (modified) (4 diffs)

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

@@ -1 @@
-using System;
-using System.Collections.Generic;
-using System.Linq;
+#region License Information
+/* HeuristicLab
+ * Copyright (C) 2002-2012 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
 
-namespace HeuristicLab.Algorithms.DataAnalysis.GaussianProcess {
-  public class CovariancePeriodic : ICovarianceFunction {
+using System;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  [StorableClass]
+  [Item(Name = "CovariancePeriodic", Description = "Periodic covariance function for Gaussian processes.")]
+  public class CovariancePeriodic : Item, ICovarianceFunction {
+    [Storable]
     private double[,] x;
+    [Storable]
     private double[,] xt;
+    [Storable]
     private double sf2;
+    [Storable]
     private double l;
-    private double[,] sd;
+    [Storable]
     private double p;
 
-    public int NumberOfParameters {
-      get { return 2; }
+    private bool symmetric;
+
+    private double[,] sd;
+    public int GetNumberOfParameters(int numberOfVariables) {
+      return 3;
+    }
+    [StorableConstructor]
+    protected CovariancePeriodic(bool deserializing) : base(deserializing) { }
+    protected CovariancePeriodic(CovariancePeriodic original, Cloner cloner)
+      : base(original, cloner) {
+      if (original.x != null) {
+        x = new double[original.x.GetLength(0), original.x.GetLength(1)];
+        Array.Copy(original.x, x, x.Length);
+        xt = new double[original.xt.GetLength(0), original.xt.GetLength(1)];
+        Array.Copy(original.xt, xt, xt.Length);
+      }
+      sf2 = original.sf2;
+      l = original.l;
+      p = original.p;
+      symmetric = original.symmetric;
+    }
+    public CovariancePeriodic()
+      : base() {
     }
 
-    public CovariancePeriodic(double p) {
-      this.p = p;
+    public override IDeepCloneable Clone(Cloner cloner) {
+      return new CovariancePeriodic(this, cloner);
     }
 
-    public void SetMatrix(double[,] x) {
-      SetMatrix(x, x);
+    public void SetParameter(double[] hyp) {
+      if (hyp.Length != 3) throw new ArgumentException();
+      this.l = Math.Exp(hyp[0]);
+      this.p = Math.Exp(hyp[1]);
+      this.sf2 = Math.Exp(2 * hyp[2]);
+
+      sf2 = Math.Min(10E6, sf2); // upper limit for the scale
+
+      sd = null;
+    }
+    public void SetData(double[,] x) {
+      SetData(x, x);
+      this.symmetric = true;
     }
 
-    public void SetMatrix(double[,] x, double[,] xt) {
+    public void SetData(double[,] x, double[,] xt) {
       this.x = x;
       this.xt = xt;
-      sd = null;
-    }
+      this.symmetric = false;
 
-    public void SetHyperparamter(double[] hyp) {
-      if (hyp.Length != 2) throw new ArgumentException();
-      this.l = Math.Exp(hyp[0]);
-      this.sf2 = Math.Exp(2 * hyp[1]);
       sd = null;
     }
@@ -53 +108 @@
       var cov = new double[rows];
       for (int i = 0; i < rows; i++) {
-        double k = Math.Sqrt(SqrDist(GetRow(x, i), GetRow(xt, i)));
+        double k = Math.Sqrt(Util.SqrDist(Util.GetRow(x, i), Util.GetRow(xt, i)));
         k = Math.PI * k / p;
         k = Math.Sin(k) / l;
@@ -62 +117 @@
     }
 
-    public double[] GetDerivatives(int i, int j) {
+    public double[] GetGradient(int i, int j) {
 
-      var res = new double[2];
+      var res = new double[3];
       double k = sd[i, j];
       k = Math.PI * k / p;
-      k = Math.Sin(k) / l;
-      k = k * k;
-      res[0] = 4 * sf2 * Math.Exp(-2 * k) * k;
-      res[1] = 2 * sf2 * Math.Exp(-2 * k);
+      {
+        double newK = Math.Sin(k) / l;
+        newK = newK * newK;
+        res[0] = 4 * sf2 * Math.Exp(-2 * newK) * newK;
+      }
+      {
+        double r = Math.Sin(k) / l;
+        res[1] = 4 * sf2 / l * Math.Exp(-2 * r * r) * r * Math.Cos(k) * k;
+      }
+      {
+        double newK = Math.Sin(k) / l;
+        newK = newK * newK;
+        res[2] = 2 * sf2 * Math.Exp(-2 * newK);
+      }
+
       return res;
     }
@@ -79 +145 @@
       int cols = xt.GetLength(0);
       sd = new double[rows, cols];
-      bool symmetric = x == xt;
-      for (int i = 0; i < rows; i++) {
-        for (int j = i; j < rows; j++) {
-          sd[i, j] = Math.Sqrt(SqrDist(GetRow(x, i), GetRow(xt, j)));
-          if (symmetric) {
+
+      if (symmetric) {
+        for (int i = 0; i < rows; i++) {
+          for (int j = i; j < cols; j++) {
+            sd[i, j] = Math.Sqrt(Util.SqrDist(Util.GetRow(x, i), Util.GetRow(x, j)));
             sd[j, i] = sd[i, j];
-          } else {
-            sd[j, i] = Math.Sqrt(SqrDist(GetRow(x, j), GetRow(xt, i)));
+          }
+        }
+      } else {
+        for (int i = 0; i < rows; i++) {
+          for (int j = 0; j < cols; j++) {
+            sd[i, j] = Math.Sqrt(Util.SqrDist(Util.GetRow(x, i), Util.GetRow(xt, j)));
           }
         }
       }
     }
-
-
-    private double SqrDist(IEnumerable<double> x, IEnumerable<double> y) {
-      var d0 = x.Zip(y, (a, b) => (a - b) * (a - b));
-      return Math.Max(0, d0.Sum());
-    }
-    private static IEnumerable<double> GetRow(double[,] x, int r) {
-      int cols = x.GetLength(1);
-      return Enumerable.Range(0, cols).Select(c => x[r, c]);
-    }
   }
 }