source: branches/DataAnalysis Refactoring/HeuristicLab.Algorithms.DataAnalysis/3.4/Linear/LinearDiscriminantAnalysis.cs @ 5678

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2011 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.Linq;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Optimization;
using HeuristicLab.Parameters;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Problems.DataAnalysis;
using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
using System.Collections.Generic;
using HeuristicLab.Problems.DataAnalysis.Symbolic;
using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;
using HeuristicLab.Problems.DataAnalysis.Symbolic.Classification;

namespace HeuristicLab.Algorithms.DataAnalysis {
  /// <summary>
  /// Linear discriminant analysis classification algorithm.
  /// </summary>
  [Item("Linear Discriminant Analysis", "Linear discriminant analysis classification algorithm.")]
  [Creatable("Data Analysis")]
  [StorableClass]
  public sealed class LinearDiscriminantAnalysis : FixedDataAnalysisAlgorithm<IClassificationProblem> {
    private const string LinearDiscriminantAnalysisSolutionResultName = "Linear discriminant analysis solution";

    [StorableConstructor]
    private LinearDiscriminantAnalysis(bool deserializing) : base(deserializing) { }
    private LinearDiscriminantAnalysis(LinearDiscriminantAnalysis original, Cloner cloner)
      : base(original, cloner) {
    }
    public LinearDiscriminantAnalysis()
      : base() {
      Problem = new ClassificationProblem();
    }
    [StorableHook(HookType.AfterDeserialization)]
    private void AfterDeserialization() { }

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

    #region Fisher LDA
    protected override void Run() {
      var solution = CreateLinearDiscriminantAnalysisSolution(Problem.ProblemData);
      Results.Add(new Result(LinearDiscriminantAnalysisSolutionResultName, "The linear discriminant analysis.", solution));
    }

    public static IClassificationSolution CreateLinearDiscriminantAnalysisSolution(IClassificationProblemData problemData) {
      Dataset dataset = problemData.Dataset;
      string targetVariable = problemData.TargetVariable;
      IEnumerable<string> allowedInputVariables = problemData.AllowedInputVariables;
      int samplesStart = problemData.TrainingPartitionStart.Value;
      int samplesEnd = problemData.TrainingPartitionEnd.Value;
      IEnumerable<int> rows = Enumerable.Range(samplesStart, samplesEnd - samplesStart);
      int nClasses = problemData.ClassNames.Count();
      double[,] inputMatrix = AlglibUtil.PrepareInputMatrix(dataset, allowedInputVariables.Concat(new string[] { targetVariable }), rows);

      // change class values into class index
      int targetVariableColumn = inputMatrix.GetLength(1) - 1;
      List<double> classValues = problemData.ClassValues.OrderBy(x => x).ToList();
      for (int row = 0; row < inputMatrix.GetLength(0); row++) {
        inputMatrix[row, targetVariableColumn] = classValues.IndexOf(inputMatrix[row, targetVariableColumn]);
      }
      int info;
      double[] w;
      alglib.fisherlda(inputMatrix, inputMatrix.GetLength(0), allowedInputVariables.Count(), nClasses, out info, out w);
      if (info < 1) throw new ArgumentException("Error in calculation of linear discriminant analysis solution");

      ISymbolicExpressionTree tree = new SymbolicExpressionTree(new ProgramRootSymbol().CreateTreeNode());
      ISymbolicExpressionTreeNode startNode = new StartSymbol().CreateTreeNode();
      tree.Root.AddSubTree(startNode);
      ISymbolicExpressionTreeNode addition = new Addition().CreateTreeNode();
      startNode.AddSubTree(addition);

      int col = 0;
      foreach (string column in allowedInputVariables) {
        VariableTreeNode vNode = (VariableTreeNode)new HeuristicLab.Problems.DataAnalysis.Symbolic.Variable().CreateTreeNode();
        vNode.VariableName = column;
        vNode.Weight = w[col];
        addition.AddSubTree(vNode);
        col++;
      }

      ConstantTreeNode cNode = (ConstantTreeNode)new Constant().CreateTreeNode();
      cNode.Value = w[w.Length - 1];
      addition.AddSubTree(cNode);


      var model = LinearDiscriminantAnalysis.CreateDiscriminantFunctionModel(tree, new SymbolicDataAnalysisExpressionTreeInterpreter(), problemData, rows);
      SymbolicDiscriminantFunctionClassificationSolution solution = new SymbolicDiscriminantFunctionClassificationSolution(model, problemData);

      return solution;
    }
    #endregion

    private static SymbolicDiscriminantFunctionClassificationModel CreateDiscriminantFunctionModel(ISymbolicExpressionTree tree,
      ISymbolicDataAnalysisExpressionTreeInterpreter interpreter,
      IClassificationProblemData problemData,
      IEnumerable<int> rows) {
      string targetVariable = problemData.TargetVariable;
      List<double> originalClasses = problemData.ClassValues.ToList();
      int nClasses = problemData.Classes;
      List<double> estimatedValues = interpreter.GetSymbolicExpressionTreeValues(tree, problemData.Dataset, rows).ToList();
      double maxEstimatedValue = estimatedValues.Max();
      double minEstimatedValue = estimatedValues.Min();
      var estimatedTargetValues =
         (from row in problemData.TrainingIndizes
          select new { EstimatedValue = estimatedValues[row], TargetValue = problemData.Dataset[targetVariable, row] })
         .ToList();

      Dictionary<double, double> classMean = new Dictionary<double, double>();
      Dictionary<double, double> classStdDev = new Dictionary<double, double>();
      // calculate moments per class
      foreach (var classValue in originalClasses) {
        var estimatedValuesForClass = from x in estimatedTargetValues
                                      where x.TargetValue == classValue
                                      select x.EstimatedValue;
        double mean, variance;
        OnlineMeanAndVarianceCalculator.Calculate(estimatedValuesForClass, out mean, out variance);
        classMean[classValue] = mean;
        classStdDev[classValue] = Math.Sqrt(variance);
      }
      List<double> thresholds = new List<double>();
      for (int i = 0; i < nClasses - 1; i++) {
        for (int j = i + 1; j < nClasses; j++) {
          double x1, x2;
          double class0 = originalClasses[i];
          double class1 = originalClasses[j];
          // calculate all thresholds
          CalculateCutPoints(classMean[class0], classStdDev[class0], classMean[class1], classStdDev[class1], out x1, out x2);
          if (!thresholds.Any(x => x.IsAlmost(x1))) thresholds.Add(x1);
          if (!thresholds.Any(x => x.IsAlmost(x2))) thresholds.Add(x2);
        }
      }
      thresholds.Sort();
      thresholds.Insert(0, double.NegativeInfinity);
      thresholds.Add(double.PositiveInfinity);
      List<double> classValues = new List<double>();
      for (int i = 0; i < thresholds.Count - 1; i++) {
        double m;
        if (double.IsNegativeInfinity(thresholds[i])) {
          m = thresholds[i + 1] - 1.0;
        } else if (double.IsPositiveInfinity(thresholds[i + 1])) {
          m = thresholds[i] + 1.0;
        } else {
          m = thresholds[i] + (thresholds[i + 1] - thresholds[i]) / 2.0;
        }

        double maxDensity = 0;
        double maxDensityClassValue = -1;
        foreach (var classValue in originalClasses) {
          double density = NormalDensity(m, classMean[classValue], classStdDev[classValue]);
          if (density > maxDensity) {
            maxDensity = density;
            maxDensityClassValue = classValue;
          }
        }
        classValues.Add(maxDensityClassValue);
      }
      List<double> filteredThresholds = new List<double>();
      List<double> filteredClassValues = new List<double>();
      filteredThresholds.Add(thresholds[0]);
      filteredClassValues.Add(classValues[0]);
      for (int i = 0; i < classValues.Count - 1; i++) {
        if (classValues[i] != classValues[i + 1]) {
          filteredThresholds.Add(thresholds[i + 1]);
          filteredClassValues.Add(classValues[i + 1]);
        }
      }
      filteredThresholds.Add(double.PositiveInfinity);

      return new SymbolicDiscriminantFunctionClassificationModel(tree, interpreter, filteredClassValues, filteredThresholds);
    }

    private static double NormalDensity(double x, double mu, double sigma) {
      return (1.0 / Math.Sqrt(2.0 * Math.PI * sigma * sigma)) * Math.Exp(-((x - mu) * (x - mu)) / (2.0 * sigma * sigma));
    }

    private static void CalculateCutPoints(double m1, double s1, double m2, double s2, out double x1, out double x2) {
      double a = (s1 * s1 - s2 * s2);
      double b = (m1 * s2 * s2 - m2 * s1 * s1);
      double c = 2 * s1 * s1 * s2 * s2 * Math.Log(s2) - 2 * s1 * s1 * s2 * s2 * Math.Log(s1) - s1 * s1 * m2 * m2 + s2 * s2 * m1 * m1;
      x1 = -(-m2 * s1 * s1 + m1 * s2 * s2 + Math.Sqrt(s1 * s1 * s2 * s2 * ((m1 - m2) * (m1 - m2) + 2.0 * (-s1 * s1 + s2 * s2) * Math.Log(s2 / s1)))) / a;
      x2 = (m2 * s1 * s1 - m1 * s2 * s2 + Math.Sqrt(s1 * s1 * s2 * s2 * ((m1 - m2) * (m1 - m2) + 2.0 * (-s1 * s1 + s2 * s2) * Math.Log(s2 / s1)))) / a;
    }
  }
}