Free cookie consent management tool by TermsFeed Policy Generator

source: branches/Breadcrumbs/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSpectralMixture.cs @ 11594

Last change on this file since 11594 was 11594, checked in by jkarder, 9 years ago

#2116: merged r10041-r11593 from trunk into branch

File size: 10.1 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2014 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 System.Collections.Generic;
24using System.Linq;
25using System.Linq.Expressions;
26using HeuristicLab.Common;
27using HeuristicLab.Core;
28using HeuristicLab.Data;
29using HeuristicLab.Parameters;
30using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
31
32namespace HeuristicLab.Algorithms.DataAnalysis {
33  [StorableClass]
34  [Item(Name = "CovarianceSpectralMixture",
35    Description = "The spectral mixture kernel described in Wilson A. G. and Adams R.P., Gaussian Process Kernels for Pattern Discovery and Exptrapolation, ICML 2013.")]
36  public sealed class CovarianceSpectralMixture : ParameterizedNamedItem, ICovarianceFunction {
37    public const string QParameterName = "Number of components (Q)";
38    public const string WeightParameterName = "Weight";
39    public const string FrequencyParameterName = "Component frequency (mu)";
40    public const string LengthScaleParameterName = "Length scale (nu)";
41    public IValueParameter<IntValue> QParameter {
42      get { return (IValueParameter<IntValue>)Parameters[QParameterName]; }
43    }
44
45    public IValueParameter<DoubleArray> WeightParameter {
46      get { return (IValueParameter<DoubleArray>)Parameters[WeightParameterName]; }
47    }
48    public IValueParameter<DoubleArray> FrequencyParameter {
49      get { return (IValueParameter<DoubleArray>)Parameters[FrequencyParameterName]; }
50    }
51
52    public IValueParameter<DoubleArray> LengthScaleParameter {
53      get { return (IValueParameter<DoubleArray>)Parameters[LengthScaleParameterName]; }
54    }
55
56    private bool HasFixedWeightParameter {
57      get { return WeightParameter.Value != null; }
58    }
59    private bool HasFixedFrequencyParameter {
60      get { return FrequencyParameter.Value != null; }
61    }
62    private bool HasFixedLengthScaleParameter {
63      get { return LengthScaleParameter.Value != null; }
64    }
65
66    [StorableConstructor]
67    private CovarianceSpectralMixture(bool deserializing)
68      : base(deserializing) {
69    }
70
71    private CovarianceSpectralMixture(CovarianceSpectralMixture original, Cloner cloner)
72      : base(original, cloner) {
73    }
74
75    public CovarianceSpectralMixture()
76      : base() {
77      Name = ItemName;
78      Description = ItemDescription;
79      Parameters.Add(new ValueParameter<IntValue>(QParameterName, "The number of Gaussians (Q) to use for the spectral mixture.", new IntValue(10)));
80      Parameters.Add(new OptionalValueParameter<DoubleArray>(WeightParameterName, "The weight of the component w (peak height of the Gaussian in spectrum)."));
81      Parameters.Add(new OptionalValueParameter<DoubleArray>(FrequencyParameterName, "The inverse component period parameter mu_q (location of the Gaussian in spectrum)."));
82      Parameters.Add(new OptionalValueParameter<DoubleArray>(LengthScaleParameterName, "The length scale parameter (nu_q) (variance of the Gaussian in the spectrum)."));
83    }
84
85    public override IDeepCloneable Clone(Cloner cloner) {
86      return new CovarianceSpectralMixture(this, cloner);
87    }
88
89    public int GetNumberOfParameters(int numberOfVariables) {
90      var q = QParameter.Value.Value;
91      return
92        (HasFixedWeightParameter ? 0 : q) +
93        (HasFixedFrequencyParameter ? 0 : q * numberOfVariables) +
94        (HasFixedLengthScaleParameter ? 0 : q * numberOfVariables);
95    }
96
97    public void SetParameter(double[] p) {
98      double[] weight, frequency, lengthScale;
99      GetParameterValues(p, out weight, out frequency, out lengthScale);
100      WeightParameter.Value = new DoubleArray(weight);
101      FrequencyParameter.Value = new DoubleArray(frequency);
102      LengthScaleParameter.Value = new DoubleArray(lengthScale);
103    }
104
105
106    private void GetParameterValues(double[] p, out double[] weight, out double[] frequency, out double[] lengthScale) {
107      // gather parameter values
108      int c = 0;
109      int q = QParameter.Value.Value;
110      // guess number of elements for frequency and length (=q * numberOfVariables)
111      int n = WeightParameter.Value == null ? ((p.Length - q) / 2) : (p.Length / 2);
112      if (HasFixedWeightParameter) {
113        weight = WeightParameter.Value.ToArray();
114      } else {
115        weight = p.Skip(c).Select(Math.Exp).Take(q).ToArray();
116        c += q;
117      }
118      if (HasFixedFrequencyParameter) {
119        frequency = FrequencyParameter.Value.ToArray();
120      } else {
121        frequency = p.Skip(c).Select(Math.Exp).Take(n).ToArray();
122        c += n;
123      }
124      if (HasFixedLengthScaleParameter) {
125        lengthScale = LengthScaleParameter.Value.ToArray();
126      } else {
127        lengthScale = p.Skip(c).Select(Math.Exp).Take(n).ToArray();
128        c += n;
129      }
130      if (p.Length != c) throw new ArgumentException("The length of the parameter vector does not match the number of free parameters for CovarianceSpectralMixture", "p");
131    }
132
133    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, IEnumerable<int> columnIndices) {
134      double[] weight, frequency, lengthScale;
135      GetParameterValues(p, out weight, out frequency, out lengthScale);
136      var fixedWeight = HasFixedWeightParameter;
137      var fixedFrequency = HasFixedFrequencyParameter;
138      var fixedLengthScale = HasFixedLengthScaleParameter;
139      // create functions
140      var cov = new ParameterizedCovarianceFunction();
141      cov.Covariance = (x, i, j) => {
142        return GetCovariance(x, x, i, j, QParameter.Value.Value, weight, frequency,
143                             lengthScale, columnIndices);
144      };
145      cov.CrossCovariance = (x, xt, i, j) => {
146        return GetCovariance(x, xt, i, j, QParameter.Value.Value, weight, frequency,
147                             lengthScale, columnIndices);
148      };
149      cov.CovarianceGradient = (x, i, j) => GetGradient(x, i, j, QParameter.Value.Value, weight, frequency,
150                             lengthScale, columnIndices, fixedWeight, fixedFrequency, fixedLengthScale);
151      return cov;
152    }
153
154    private static double GetCovariance(double[,] x, double[,] xt, int i, int j, int maxQ, double[] weight, double[] frequency, double[] lengthScale, IEnumerable<int> columnIndices) {
155      // tau = x - x' (only for selected variables)
156      double[] tau =
157        Util.GetRow(x, i, columnIndices).Zip(Util.GetRow(xt, j, columnIndices), (xi, xj) => xi - xj).ToArray();
158      int numberOfVariables = lengthScale.Length / maxQ;
159      double k = 0;
160      // for each component
161      for (int q = 0; q < maxQ; q++) {
162        double kc = weight[q]; // weighted kernel component
163
164        int idx = 0; // helper index for tau
165        // for each selected variable
166        foreach (var c in columnIndices) {
167          kc *= f1(tau[idx], lengthScale[q * numberOfVariables + c]) * f2(tau[idx], frequency[q * numberOfVariables + c]);
168          idx++;
169        }
170        k += kc;
171      }
172      return k;
173    }
174
175    public static double f1(double tau, double lengthScale) {
176      return Math.Exp(-2 * Math.PI * Math.PI * tau * tau * lengthScale);
177    }
178    public static double f2(double tau, double frequency) {
179      return Math.Cos(2 * Math.PI * tau * frequency);
180    }
181
182    // order of returned gradients must match the order in GetParameterValues!
183    private static IEnumerable<double> GetGradient(double[,] x, int i, int j, int maxQ, double[] weight, double[] frequency, double[] lengthScale, IEnumerable<int> columnIndices,
184      bool fixedWeight, bool fixedFrequency, bool fixedLengthScale) {
185      double[] tau = Util.GetRow(x, i, columnIndices).Zip(Util.GetRow(x, j, columnIndices), (xi, xj) => xi - xj).ToArray();
186      int numberOfVariables = lengthScale.Length / maxQ;
187
188      if (!fixedWeight) {
189        // weight
190        // for each component
191        for (int q = 0; q < maxQ; q++) {
192          double k = weight[q];
193          int idx = 0; // helper index for tau
194          // for each selected variable
195          foreach (var c in columnIndices) {
196            k *= f1(tau[idx], lengthScale[q * numberOfVariables + c]) * f2(tau[idx], frequency[q * numberOfVariables + c]);
197            idx++;
198          }
199          yield return k;
200        }
201      }
202
203      if (!fixedFrequency) {
204        // frequency
205        // for each component
206        for (int q = 0; q < maxQ; q++) {
207          int idx = 0; // helper index for tau
208          // for each selected variable
209          foreach (var c in columnIndices) {
210            double k = f1(tau[idx], lengthScale[q * numberOfVariables + c]) *
211                       -2 * Math.PI * tau[idx] * frequency[q * numberOfVariables + c] *
212                       Math.Sin(2 * Math.PI * tau[idx] * frequency[q * numberOfVariables + c]);
213            idx++;
214            yield return weight[q] * k;
215          }
216        }
217      }
218
219      if (!fixedLengthScale) {
220        // length scale
221        // for each component
222        for (int q = 0; q < maxQ; q++) {
223          int idx = 0; // helper index for tau
224          // for each selected variable
225          foreach (var c in columnIndices) {
226            double k = -2 * Math.PI * Math.PI * tau[idx] * tau[idx] * lengthScale[q * numberOfVariables + c] *
227                       f1(tau[idx], lengthScale[q * numberOfVariables + c]) *
228                       f2(tau[idx], frequency[q * numberOfVariables + c]);
229            idx++;
230            yield return weight[q] * k;
231          }
232        }
233      }
234    }
235  }
236}
Note: See TracBrowser for help on using the repository browser.