Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovariancePeriodic.cs @ 8473

Last change on this file since 8473 was 8473, checked in by gkronber, 12 years ago

#1902 worked on GPR: added line chart, made parameters of mean and covariance functions readable, removed target variable scaling, moved noise hyperparameter for likelihood function to the end of the parameter list, added methods to calculate the predicted variance, removed limits for scale of covariance functions and introduced exception handling to catch non-spd or singular cov matrixes, implemented rational quadratic covariance function, added unit test case from GBML book (however it does not work as the book seemingly uses a noise-less likelihood function)

File size: 4.7 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2012 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using System;
23using HeuristicLab.Common;
24using HeuristicLab.Core;
25using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
26
27namespace HeuristicLab.Algorithms.DataAnalysis {
28  [StorableClass]
29  [Item(Name = "CovariancePeriodic", Description = "Periodic covariance function for Gaussian processes.")]
30  public class CovariancePeriodic : Item, ICovarianceFunction {
31    [Storable]
32    private double[,] x;
33    [Storable]
34    private double[,] xt;
35    [Storable]
36    private double sf2;
37    public double Scale { get { return sf2; } }
38    [Storable]
39    private double l;
40    public double Length { get { return l; } }
41    [Storable]
42    private double p;
43    public double Period { get { return p; } }
44
45    private bool symmetric;
46
47    private double[,] sd;
48    public int GetNumberOfParameters(int numberOfVariables) {
49      return 3;
50    }
51    [StorableConstructor]
52    protected CovariancePeriodic(bool deserializing) : base(deserializing) { }
53    protected CovariancePeriodic(CovariancePeriodic original, Cloner cloner)
54      : base(original, cloner) {
55      if (original.x != null) {
56        x = new double[original.x.GetLength(0), original.x.GetLength(1)];
57        Array.Copy(original.x, x, x.Length);
58        xt = new double[original.xt.GetLength(0), original.xt.GetLength(1)];
59        Array.Copy(original.xt, xt, xt.Length);
60      }
61      sf2 = original.sf2;
62      l = original.l;
63      p = original.p;
64      symmetric = original.symmetric;
65    }
66    public CovariancePeriodic()
67      : base() {
68    }
69
70    public override IDeepCloneable Clone(Cloner cloner) {
71      return new CovariancePeriodic(this, cloner);
72    }
73
74    public void SetParameter(double[] hyp) {
75      if (hyp.Length != 3) throw new ArgumentException();
76      this.l = Math.Exp(hyp[0]);
77      this.p = Math.Exp(hyp[1]);
78      this.sf2 = Math.Exp(2 * hyp[2]);
79      // sf2 = Math.Min(10E6, sf2); // upper limit for the scale
80
81      sd = null;
82    }
83    public void SetData(double[,] x) {
84      SetData(x, x);
85      this.symmetric = true;
86    }
87
88    public void SetData(double[,] x, double[,] xt) {
89      this.x = x;
90      this.xt = xt;
91      this.symmetric = false;
92
93      sd = null;
94    }
95
96    public double GetCovariance(int i, int j) {
97      if (sd == null) CalculateSquaredDistances();
98      double k = sd[i, j];
99      k = Math.PI * k / p;
100      k = Math.Sin(k) / l;
101      k = k * k;
102
103      return sf2 * Math.Exp(-2.0 * k);
104    }
105
106    public double GetGradient(int i, int j, int k) {
107      double v = Math.PI * sd[i, j] / p;
108      switch (k) {
109        case 0: {
110            double newK = Math.Sin(v) / l;
111            newK = newK * newK;
112            return 4 * sf2 * Math.Exp(-2 * newK) * newK;
113          }
114        case 1: {
115            double r = Math.Sin(v) / l;
116            return 4 * sf2 / l * Math.Exp(-2 * r * r) * r * Math.Cos(v) * v;
117          }
118        case 2: {
119            double newK = Math.Sin(v) / l;
120            newK = newK * newK;
121            return 2 * sf2 * Math.Exp(-2 * newK);
122
123          }
124        default: {
125            throw new ArgumentException("CovariancePeriodic only has three hyperparameters.", "k");
126          }
127      }
128    }
129
130    private void CalculateSquaredDistances() {
131      if (x.GetLength(1) != xt.GetLength(1)) throw new InvalidOperationException();
132      int rows = x.GetLength(0);
133      int cols = xt.GetLength(0);
134      sd = new double[rows, cols];
135
136      if (symmetric) {
137        for (int i = 0; i < rows; i++) {
138          for (int j = i; j < cols; j++) {
139            sd[i, j] = Math.Sqrt(Util.SqrDist(Util.GetRow(x, i), Util.GetRow(x, j)));
140            sd[j, i] = sd[i, j];
141          }
142        }
143      } else {
144        for (int i = 0; i < rows; i++) {
145          for (int j = 0; j < cols; j++) {
146            sd[i, j] = Math.Sqrt(Util.SqrDist(Util.GetRow(x, i), Util.GetRow(xt, j)));
147          }
148        }
149      }
150    }
151  }
152}
Note: See TracBrowser for help on using the repository browser.