source: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/RadialBasisFunctions/RadialBasisFunctionModel.cs @ 14386

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2016 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.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Problems.DataAnalysis;

namespace HeuristicLab.Algorithms.DataAnalysis {
  /// <summary>
  /// Represents a Radial Basis Function regression model.
  /// </summary>
  [StorableClass]
  [Item("RBFModel", "Represents a Gaussian process posterior.")]
  public sealed class RadialBasisFunctionModel : RegressionModel, IConfidenceRegressionModel {
    public override IEnumerable<string> VariablesUsedForPrediction
    {
      get { return allowedInputVariables; }
    }

    [Storable]
    private string[] allowedInputVariables;
    public string[] AllowedInputVariables
    {
      get { return allowedInputVariables; }
    }

    [Storable]
    private double[] alpha;
    [Storable]
    private IDataset trainingDataset; // it is better to store the original training dataset completely because this is more efficient in persistence
    [Storable]
    private int[] trainingRows;
    [Storable]
    private IKernelFunction<double[]> kernel;
    [Storable]
    private DoubleMatrix gramInvert;


    [StorableConstructor]
    private RadialBasisFunctionModel(bool deserializing) : base(deserializing) { }
    private RadialBasisFunctionModel(RadialBasisFunctionModel original, Cloner cloner)
      : base(original, cloner) {
      // shallow copies of arrays because they cannot be modified
      trainingRows = original.trainingRows;
      allowedInputVariables = original.allowedInputVariables;
      alpha = original.alpha;
      trainingDataset = original.trainingDataset;
      kernel = original.kernel;
    }
    public RadialBasisFunctionModel(IDataset dataset, string targetVariable, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows, IKernelFunction<double[]> kernel)
      : base(targetVariable) {
      name = ItemName;
      description = ItemDescription;
      this.allowedInputVariables = allowedInputVariables.ToArray();
      trainingRows = rows.ToArray();
      trainingDataset = dataset;
      this.kernel = (IKernelFunction<double[]>)kernel.Clone();
      try {
        var data = ExtractData(dataset, trainingRows);
        var qualities = dataset.GetDoubleValues(targetVariable, trainingRows).ToArray();
        int info;
        alglib.densesolverlsreport denseSolveRep;
        var gr = BuildGramMatrix(data);
        alglib.rmatrixsolvels(gr, data.Length + 1, data.Length + 1, qualities.Concat(new[] { 0.0 }).ToArray(), 0.0, out info, out denseSolveRep, out alpha);
        if (info != 1) throw new ArgumentException("Could not create Model.");
        gramInvert = new DoubleMatrix(gr).Invert();
      }
      catch (alglib.alglibexception ae) {
        // wrap exception so that calling code doesn't have to know about alglib implementation
        throw new ArgumentException("There was a problem in the calculation of the RBF process model", ae);
      }
    }
    private double[][] ExtractData(IDataset dataset, IEnumerable<int> rows) {
      return rows.Select(r => allowedInputVariables.Select(v => dataset.GetDoubleValue(v, r)).ToArray()).ToArray();
    }

    public override IDeepCloneable Clone(Cloner cloner) {
      return new RadialBasisFunctionModel(this, cloner);
    }

    #region IRegressionModel Members
    public override IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
      var solutions = ExtractData(dataset, rows);
      var data = ExtractData(trainingDataset, trainingRows);
      return solutions.Select(solution => alpha.Zip(data, (a, d) => a * kernel.Get(solution, d)).Sum() + 1 * alpha[alpha.Length - 1]).ToArray();
    }
    public override IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
      return new RadialBasisFunctionRegressionSolution(this, new RegressionProblemData(problemData));
    }
    #endregion

    public IEnumerable<double> GetEstimatedVariances(IDataset dataset, IEnumerable<int> rows) {
      var data = ExtractData(trainingDataset, trainingRows);
      return ExtractData(dataset, rows).Select(x => GetVariance(x, data));
    }
    public double LeaveOneOutCrossValidationRootMeanSquaredError() {
      return Math.Sqrt(alpha.Select((t, i) => t / gramInvert[i, i]).Sum(d => d * d) / gramInvert.Rows);
    }

    #region helpers
    private double[,] BuildGramMatrix(double[][] data) {
      var size = data.Length + 1;
      var gram = new double[size, size];
      for (var i = 0; i < size; i++)
        for (var j = i; j < size; j++) {
          if (j == size - 1 && i == size - 1) gram[i, j] = 0;
          else if (j == size - 1 || i == size - 1) gram[j, i] = gram[i, j] = 1;
          else gram[j, i] = gram[i, j] = kernel.Get(data[i], data[j]); //symmteric Matrix --> half of the work
        }
      return gram;
    }
    private double GetVariance(double[] solution, IEnumerable<double[]> data) {
      var phiT = data.Select(x => kernel.Get(x, solution)).Concat(new[] { 1.0 }).ToColumnVector();
      var res = phiT.Transpose().Mul(gramInvert.Mul(phiT))[0, 0];
      return res > 0 ? res : 0;
    }
    #endregion
  }
}