source: branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Util.cs @ 12130

#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

using System;
using System.Collections.Generic;
using System.Linq;
using HeuristicLab.Algorithms.DataAnalysis;
using HeuristicLab.Core;
using HeuristicLab.Random;

namespace HeuristicLab.Problems.Instances.DataAnalysis {
  public class Util {


    public static List<double> SampleGaussianProcess(IRandom random, ParameterizedCovarianceFunction covFunction, List<List<double>> data) {
      int n = data[0].Count;

      var normalRand = new NormalDistributedRandom(random, 0, 1);
      var alpha = (from i in Enumerable.Range(0, n)
                   select normalRand.NextDouble()).ToArray();
      return SampleGaussianProcess(random, covFunction, data, alpha);
    }

    public static List<double> SampleGaussianProcess(IRandom random, ParameterizedCovarianceFunction covFunction, List<List<double>> data, double[] alpha) {
      if (alpha.Length != data[0].Count) throw new ArgumentException();

      double[,] x = new double[data[0].Count, data.Count];
      for (int i = 0; i < x.GetLength(0); i++)
        for (int j = 0; j < x.GetLength(1); j++)
          x[i, j] = data[j][i];
      double[,] K = new double[x.GetLength(0), x.GetLength(0)];
      for (int i = 0; i < K.GetLength(0); i++)
        for (int j = i; j < K.GetLength(1); j++)
          K[i, j] = covFunction.Covariance(x, i, j);

      // if (!alglib.spdmatrixcholesky(ref K, K.GetLength(0), true)) throw new ArgumentException();
      K = toeplitz_cholesky_lower(K.GetLength(0), K);
      List<double> target = new List<double>(K.GetLength(0));
      for (int i = 0; i < K.GetLength(0); i++) {
        double s = 0.0;
        for (int j = K.GetLength(0) - 1; j >= 0; j--) {

          s += K[j, i] * alpha[j];
        }

        target.Add(s);
      }

      return target;
    }

    //****************************************************************************80

    public static double[,] toeplitz_cholesky_lower(int n, double[,] a)

    //****************************************************************************80
      //
      //  Purpose:
      //
      //    TOEPLITZ_CHOLESKY_LOWER: lower Cholesky factor of a Toeplitz matrix.
      //
      //  Discussion:
      //
      //    The Toeplitz matrix must be positive semi-definite.
      //
      //    After factorization, A = L * L'.
      //
      //  Licensing:
      //
      //    This code is distributed under the GNU LGPL license.
      //
      //  Modified:
      //
      //    13 November 2012
      //    29 January  2014: adapted to C# by Gabriel Kronberger
      //  Author:
      //
      //    John Burkardt
      //
      //  Reference:
      //
      //    Michael Stewart,
      //    Cholesky factorization of semi-definite Toeplitz matrices.
      //
      //  Parameters:
      //
      //    Input, int N, the order of the matrix.
      //
      //    Input, double A[N,N], the Toeplitz matrix.
      //
      //    Output, double TOEPLITZ_CHOLESKY_LOWER[N,N], the lower Cholesky factor.
      //
    {
      double div;
      double[] g;
      double g1j;
      double g2j;
      int i;
      int j;
      double[,] l;
      double rho;

      l = new double[n, n];

      for (j = 0; j < n; j++) {
        for (i = 0; i < n; i++) {
          l[i, j] = 0.0;
        }
      }
      g = new double[2 * n];

      for (j = 0; j < n; j++) {
        g[0 + j * 2] = a[0, j];
      }
      g[1 + 0 * 2] = 0.0;
      for (j = 1; j < n; j++) {
        g[1 + j * 2] = a[j, 0];
      }

      for (i = 0; i < n; i++) {
        l[i, 0] = g[0 + i * 2];
      }
      for (j = n - 1; 1 <= j; j--) {
        g[0 + j * 2] = g[0 + (j - 1) * 2];
      }
      g[0 + 0 * 2] = 0.0;

      for (i = 1; i < n; i++) {
        rho = -g[1 + i * 2] / g[0 + i * 2];
        div = Math.Sqrt((1.0 - rho) * (1.0 + rho));

        for (j = i; j < n; j++) {
          g1j = g[0 + j * 2];
          g2j = g[1 + j * 2];
          g[0 + j * 2] = (g1j + rho * g2j) / div;
          g[1 + j * 2] = (rho * g1j + g2j) / div;
        }

        for (j = i; j < n; j++) {
          l[j, i] = g[0 + j * 2];
        }
        for (j = n - 1; i < j; j--) {
          g[0 + j * 2] = g[0 + (j - 1) * 2];
        }
        g[0 + i * 2] = 0.0;
      }


      return l;
    }
    //****************************************************************************80

    public static double[,] toeplitz_cholesky_upper(int n, double[,] a)

    //****************************************************************************80
      //
      //  Purpose:
      //
      //    TOEPLITZ_CHOLESKY_UPPER: upper Cholesky factor of a Toeplitz matrix.
      //
      //  Discussion:
      //
      //    The Toeplitz matrix must be positive semi-definite.
      //
      //    After factorization, A = R' * R.
      //
      //  Licensing:
      //
      //    This code is distributed under the GNU LGPL license.
      //
      //  Modified:
      //
      //    14 November 2012
      //    29 January  2014: adapted to C# by Gabriel Kronberger
      //
      //  Author:
      //
      //    John Burkardt
      //
      //  Reference:
      //
      //    Michael Stewart,
      //    Cholesky factorization of semi-definite Toeplitz matrices.
      //
      //  Parameters:
      //
      //    Input, int N, the order of the matrix.
      //
      //    Input, double A[N,N], the Toeplitz matrix.
      //
      //    Output, double TOEPLITZ_CHOLESKY_UPPER[N,N], the upper Cholesky factor.
      //
    {
      double div;
      double[] g;
      double g1j;
      double g2j;
      int i;
      int j;
      double[,] r;
      double rho;

      r = new double[n, n];

      for (j = 0; j < n; j++) {
        for (i = 0; i < n; i++) {
          r[i, j] = 0.0;
        }
      }
      g = new double[2 * n];

      for (j = 0; j < n; j++) {
        g[0 + j * 2] = a[0, j];
      }

      g[1 + 0 * 2] = 0.0;
      for (j = 1; j < n; j++) {
        g[1 + j * 2] = a[j, 0];
      }
      for (j = 0; j < n; j++) {
        r[0, j] = g[0 + j * 2];
      }
      for (j = n - 1; 1 <= j; j--) {
        g[0 + j * 2] = g[0 + (j - 1) * 2];
      }
      g[0 + 0 * 2] = 0.0;

      for (i = 1; i < n; i++) {
        rho = -g[1 + i * 2] / g[0 + i * 2];
        div = Math.Sqrt((1.0 - rho) * (1.0 + rho));
        for (j = i; j < n; j++) {
          g1j = g[0 + j * 2];
          g2j = g[1 + j * 2];
          g[0 + j * 2] = (g1j + rho * g2j) / div;
          g[1 + j * 2] = (rho * g1j + g2j) / div;
        }
        for (j = i; j < n; j++) {
          r[i, j] = g[0 + j * 2];
        }
        for (j = n - 1; i < j; j--) {
          g[0 + j * 2] = g[0 + (j - 1) * 2];
        }
        g[0 + i * 2] = 0.0;
      }

      return r;
    }
  }
}