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

#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.Diagnostics;
using System.Linq;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Problems.DataAnalysis;

namespace HeuristicLab.Algorithms.DataAnalysis {
  /// <summary>
  /// Represents an RBF regression model.
  /// </summary>
  [StorableClass]
  [Item("RBFModel", "An RBF regression model")]
  public sealed class RadialBasisFunctionModel : RegressionModel, IConfidenceRegressionModel {
    public override IEnumerable<string> VariablesUsedForPrediction {
      get { return allowedInputVariables; }
    }

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

    [Storable]
    private readonly double[] alpha;

    [Storable]
    private readonly double[,] trainX; // it is better to store the original training dataset completely because this is more efficient in persistence

    [Storable]
    private readonly ITransformation<double>[] scaling;

    [Storable]
    private readonly ICovarianceFunction kernel;

    private double[,] gramInvert; // not storable as it can be large (recreate after deserialization as required)

    [Storable]
    private readonly double meanOffset; // implementation works for zero-mean target variables

    [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
      allowedInputVariables = original.allowedInputVariables;
      alpha = original.alpha;
      trainX = original.trainX;
      gramInvert = original.gramInvert;
      scaling = original.scaling;

      meanOffset = original.meanOffset;
      if (original.kernel != null)
        kernel = cloner.Clone(original.kernel);
    }
    public RadialBasisFunctionModel(IDataset dataset, string targetVariable, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows,
      bool scaleInputs, ICovarianceFunction kernel)
      : base(targetVariable) {
      if (kernel.GetNumberOfParameters(allowedInputVariables.Count()) > 0) throw new ArgumentException("All parameters in the kernel function must be specified.");
      name = ItemName;
      description = ItemDescription;
      this.allowedInputVariables = allowedInputVariables.ToArray();
      var trainingRows = rows.ToArray();
      this.kernel = (ICovarianceFunction)kernel.Clone();
      try {
        if (scaleInputs)
          scaling = CreateScaling(dataset, trainingRows);
        trainX = ExtractData(dataset, trainingRows, scaling);
        var y = dataset.GetDoubleValues(targetVariable, trainingRows).ToArray();
        meanOffset = y.Average();
        for (int i = 0; i < y.Length; i++) y[i] -= meanOffset;
        int info;
        // TODO: improve efficiency by decomposing matrix once instead of solving the system and then inverting the matrix
        alglib.densesolverlsreport denseSolveRep;
        gramInvert = BuildGramMatrix(trainX);
        int n = trainX.GetLength(0);
        alglib.rmatrixsolvels(gramInvert, n, n, y, 0.0, out info, out denseSolveRep, out alpha);
        if (info != 1) throw new ArgumentException("Could not create model.");

        alglib.matinvreport report;
        alglib.rmatrixinverse(ref gramInvert, out info, out report);
        if (info != 1) throw new ArgumentException("Could not invert matrix. Is it quadratic symmetric positive definite?");

      } 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 model", ae);
      }
    }

    private ITransformation<double>[] CreateScaling(IDataset dataset, int[] rows) {
      var trans = new ITransformation<double>[allowedInputVariables.Length];
      int i = 0;
      foreach (var variable in allowedInputVariables) {
        var lin = new LinearTransformation(allowedInputVariables);
        var max = dataset.GetDoubleValues(variable, rows).Max();
        var min = dataset.GetDoubleValues(variable, rows).Min();
        lin.Multiplier = 1.0 / (max - min);
        lin.Addend = -min / (max - min);
        trans[i] = lin;
        i++;
      }
      return trans;
    }

    private double[,] ExtractData(IDataset dataset, IEnumerable<int> rows, ITransformation<double>[] scaling = null) {
      double[][] variables;
      if (scaling != null) {
        variables =
          allowedInputVariables.Select((var, i) => scaling[i].Apply(dataset.GetDoubleValues(var, rows)).ToArray())
            .ToArray();
      } else {
        variables =
        allowedInputVariables.Select(var => dataset.GetDoubleValues(var, rows).ToArray()).ToArray();
      }
      int n = variables.First().Length;
      var res = new double[n, variables.Length];
      for (int r = 0; r < n; r++)
        for (int c = 0; c < variables.Length; c++) {
          res[r, c] = variables[c][r];
        }
      return res;
    }

    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 newX = ExtractData(dataset, rows, scaling);
      var dim = newX.GetLength(1);
      var cov = kernel.GetParameterizedCovarianceFunction(new double[0], Enumerable.Range(0, dim).ToArray());

      var pred = new double[newX.GetLength(0)];
      for (int i = 0; i < pred.Length; i++) {
        double sum = meanOffset;
        for (int j = 0; j < alpha.Length; j++) {
          sum += alpha[j] * cov.CrossCovariance(trainX, newX, j, i);
        }
        pred[i] = sum;
      }
      return pred;
    }
    public override IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
      return new ConfidenceRegressionSolution(this, new RegressionProblemData(problemData));
    }
    #endregion

    public IEnumerable<double> GetEstimatedVariances(IDataset dataset, IEnumerable<int> rows) {
      if (gramInvert == null) CreateAndInvertGramMatrix();
      int n = gramInvert.GetLength(0);
      var newData = ExtractData(dataset, rows, scaling);
      var dim = newData.GetLength(1);
      var cov = kernel.GetParameterizedCovarianceFunction(new double[0], Enumerable.Range(0, dim).ToArray());

      // TODO perf (matrix matrix multiplication)
      for (int i = 0; i < newData.GetLength(0); i++) {
        double[] p = new double[n];

        for (int j = 0; j < trainX.GetLength(0); j++) {
          p[j] = cov.CrossCovariance(trainX, newData, j, i);
        }

        var Ap = new double[n];
        alglib.ablas.rmatrixmv(n, n, gramInvert, 0, 0, 0, p, 0, ref Ap, 0);
        var res = 0.0;
        // dot product
        for (int j = 0; j < p.Length; j++) res += p[j] * Ap[j];
        yield return res > 0 ? res : 0;
      }
    }
    public double LeaveOneOutCrossValidationRootMeanSquaredError() {
      if (gramInvert == null) CreateAndInvertGramMatrix();
      var n = gramInvert.GetLength(0);
      var s = 1.0 / n;

      var sum = 0.0;
      for (int i = 0; i < alpha.Length; i++) {
        var x = alpha[i] / gramInvert[i, i];
        sum += x * x;
      }
      sum *= s;
      return Math.Sqrt(sum);
    }

    private void CreateAndInvertGramMatrix() {
      try {
        gramInvert = BuildGramMatrix(trainX);
        int info = 0;
        alglib.matinvreport report;
        alglib.rmatrixinverse(ref gramInvert, out info, out report);
        if (info != 1)
          throw new ArgumentException("Could not invert matrix. Is it quadratic symmetric positive definite?");
      } catch (alglib.alglibexception) {
        // wrap exception so that calling code doesn't have to know about alglib implementation
        throw new ArgumentException("Could not invert matrix. Is it quadratic symmetric positive definite?");
      }
    }
    #region helpers
    private double[,] BuildGramMatrix(double[,] data) {
      var n = data.GetLength(0);
      var dim = data.GetLength(1);
      var cov = kernel.GetParameterizedCovarianceFunction(new double[0], Enumerable.Range(0, dim).ToArray());
      var gram = new double[n, n];
      for (var i = 0; i < n; i++)
        for (var j = i; j < n; j++) {
          gram[j, i] = gram[i, j] = cov.Covariance(data, i, j); // symmetric matrix --> half of the work
        }
      return gram;
    }

    #endregion
  }
}