source: branches/M5Regression/HeuristicLab.Algorithms.DataAnalysis/3.4/PCA/PrincipleComponentAnalysisStatic.cs @ 15430

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2017 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.Collections.Generic;
using System.Linq;
using HeuristicLab.Common;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Problems.DataAnalysis;

namespace HeuristicLab.Algorithms.DataAnalysis {
  [StorableClass]
  public class PrincipleComponentAnalysisStatic : IDeepCloneable {
    #region Properties
    [Storable]
    private double[,] Matrix { get; set; }
    [Storable]
    public double[] Variances { get; private set; }
    [Storable]
    public string[] Variables { get; private set; }
    [Storable]
    private double[] Deviations { get; set; }
    [Storable]
    private double[] Means { get; set; }
    public string[] Names {
      get { return Variables.Select((_, x) => "pc" + x).ToArray(); }
    }
    #endregion

    #region HLConstructors
    [StorableConstructor]
    protected PrincipleComponentAnalysisStatic(bool deserializing) { }
    protected PrincipleComponentAnalysisStatic(PrincipleComponentAnalysisStatic original, Cloner cloner) {
      if (original.Variances != null) Variances = original.Variances.ToArray();
      if (original.Variables != null) Variables = original.Variables.ToArray();
      if (original.Matrix == null) return;
      Matrix = new double[original.Matrix.GetLength(0), original.Matrix.GetLength(1)];
      for (var i = 0; i < original.Matrix.GetLength(0); i++)
      for (var j = 0; j < original.Matrix.GetLength(1); j++)
        Matrix[i, j] = original.Matrix[i, j];
    }
    private PrincipleComponentAnalysisStatic() { }
    public IDeepCloneable Clone(Cloner cloner) {
      return new PrincipleComponentAnalysisStatic(this, cloner);
    }
    public object Clone() {
      return new Cloner().Clone(this);
    }
    #endregion

    #region Static Interface
    public static PrincipleComponentAnalysisStatic Create(IDataset dataset, IEnumerable<int> rows, IEnumerable<string> variables, bool normalize = false) {
      var res = new PrincipleComponentAnalysisStatic();
      res.BuildPca(dataset, rows, variables, normalize);
      return res;
    }
    #endregion

    public IRegressionProblemData ProjectProblem(IRegressionProblemData pd) {
      var projected = ProjectData(pd.Dataset, pd.AllIndices);
      return CreateProblemData(pd, projected);
    }
    public IDataset ProjectDataset(IDataset data) {
      var projected = ProjectData(data, Enumerable.Range(0, data.Rows));
      return CreateDataset(data, projected);
    }
    public double[,] ProjectData(IDataset dataset, IEnumerable<int> rows) {
      var instances = rows.ToArray();
      var result = new double[instances.Length, Variables.Length];
      for (var r = 0; r < instances.Length; r++)
      for (var i = 0; i < Variables.Length; i++) {
        var val = (dataset.GetDoubleValue(Variables[i], instances[r]) - Means[i]) / Deviations[i];
        for (var j = 0; j < Variables.Length; j++)
          result[r, j] += val * Matrix[i, j];
      }
      return result;
    }
    public double[,] ReverseProjection(IDataset dataset, IEnumerable<int> rows, IEnumerable<string> pcaComponents) {
      var instances = rows.ToArray();
      var components = pcaComponents.ToArray();

      var result = new double[instances.Length, Variables.Length];
      for (var r = 0; r < instances.Length; r++)
      for (var i = 0; i < components.Length; i++) {
        var val = dataset.GetDoubleValue(components[i], instances[r]);
        for (var j = 0; j < Variables.Length; j++)
          result[r, j] += (val + Means[j]) * Deviations[j] * Matrix[j, i];
      }
      return result;
    }

    #region Helpers
    private IRegressionProblemData CreateProblemData(IRegressionProblemData pd, double[,] pcs) {
      var data = CreateDataset(pd.Dataset, pcs);
      var res = new RegressionProblemData(data, Names, pd.TargetVariable);
      res.TestPartition.Start = pd.TestPartition.Start;
      res.TestPartition.End = pd.TestPartition.End;
      res.TrainingPartition.Start = pd.TrainingPartition.Start;
      res.TrainingPartition.End = pd.TrainingPartition.End;
      return res;
    }
    private IDataset CreateDataset(IDataset data, double[,] pcs) {
      var n = Names;
      var n2 = data.DoubleVariables.Where(x => !Variables.Contains(x)).ToArray();
      return new Dataset(n.Concat(n2),
        n.Select((_, x) => Enumerable.Range(0, pcs.GetLength(0)).Select(r => pcs[r, x]).ToList())
          .Concat(n2.Select(x => data.GetDoubleValues(x).ToList())));
    }
    private void BuildPca(IDataset dataset, IEnumerable<int> rows, IEnumerable<string> variables, bool normalize) {
      var instances = rows.ToArray();
      var attributes = variables.ToArray();
      Means = normalize
        ? attributes.Select(v => dataset.GetDoubleValues(v, instances).Average()).ToArray()
        : attributes.Select(x => 0.0).ToArray();
      Deviations = normalize
        ? attributes.Select(v => dataset.GetDoubleValues(v, instances).StandardDeviationPop()).Select(x => x.IsAlmost(0.0) ? 1 : x).ToArray()
        : attributes.Select(x => 1.0).ToArray();

      var data = new double[instances.Length, attributes.Length];

      for (var j = 0; j < attributes.Length; j++) {
        var i = 0;
        foreach (var v in dataset.GetDoubleValues(attributes[j], instances)) {
          data[i, j] = (v - Means[j]) / Deviations[j];
          i++;
        }
      }

      int info;
      double[] variances;
      double[,] matrix;
      alglib.pcabuildbasis(data, instances.Length, attributes.Length, out info, out variances, out matrix);
      Matrix = matrix;
      Variances = variances;
      Variables = attributes;
    }
    #endregion
  }
}