source: branches/2789_MathNetNumerics-Exploration/HeuristicLab.Algorithms.DataAnalysis.Experimental/LinearRegressionCombinations.cs @ 16310

#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.Concurrent;
using System.Collections.Generic;
using System.Linq;
using System.Threading;
using System.Threading.Tasks;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
using HeuristicLab.Optimization;
using HeuristicLab.Parameters;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Problems.DataAnalysis;
using HeuristicLab.Problems.DataAnalysis.Symbolic;
using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;

namespace HeuristicLab.Algorithms.DataAnalysis.Experimental {
  /// <summary>
  /// Linear regression data analysis algorithm.
  /// </summary>
  [Item("Linear Regression Combinations (LR)", "Calculates all possible LR solutions.")]
  [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 102)]
  [StorableClass]
  public sealed class LinearRegressionCombinations : FixedDataAnalysisAlgorithm<IRegressionProblem> {
    public IFixedValueParameter<IntValue> MaximumInputsParameter {
      get { return (IFixedValueParameter<IntValue>)Parameters["Maximum Inputs"]; }
    }
    public int MaximumInputs {
      get { return MaximumInputsParameter.Value.Value; }
      set { MaximumInputsParameter.Value.Value = value; }
    }

    public IFixedValueParameter<BoolValue> CreateSolutionParameter {
      get { return (IFixedValueParameter<BoolValue>)Parameters["Create Solution"]; }
    }
    public bool CreateSolution {
      get { return CreateSolutionParameter.Value.Value; }
      set { CreateSolutionParameter.Value.Value = value; }
    }

    public IFixedValueParameter<IntValue> MaximumSolutionsParameter {
      get { return (IFixedValueParameter<IntValue>)Parameters["Maximum Solutions stored"]; }
    }
    public int MaximumSolutions {
      get { return MaximumSolutionsParameter.Value.Value; }
      set { MaximumSolutionsParameter.Value.Value = value; }
    }

    private IntValue CalculatedModelsResults {
      get {
        if (!Results.ContainsKey("Calculated Models")) Results.Add(new Result("Calculated Models", "The number of calculated linear models ", new IntValue(0)));
        return (IntValue)Results["Calculated Models"].Value;
      }
    }
    public int CalculatedModels {
      get { return CalculatedModelsResults.Value; }
      set { CalculatedModelsResults.Value = value; }
    }

    private IntValue TotalModelsResult {
      get {
        if (!Results.ContainsKey("Total Models")) Results.Add(new Result("Total Models", "The total number of linear models to calculate", new IntValue(0)));
        return (IntValue)Results["Total Models"].Value;
      }
    }
    public int TotalModels {
      get { return TotalModelsResult.Value; }
      set { TotalModelsResult.Value = value; }
    }

    private IntValue CalculatedInputResults {
      get {
        if (!Results.ContainsKey("Calculated Inputs")) Results.Add(new Result("Calculated Inputs", "The maximum of already calculated input combinations.", new IntValue(0)));
        return (IntValue)Results["Calculated Inputs"].Value;
      }
    }
    public int CalculatedInputs {
      get { return CalculatedInputResults.Value; }
      set { CalculatedInputResults.Value = value; }
    }

    [StorableConstructor]
    private LinearRegressionCombinations(bool deserializing) : base(deserializing) { }
    [StorableHook(HookType.AfterDeserialization)]
    private void AfterDeserialization() {
      RegisterEventHandlers();
    }

    private LinearRegressionCombinations(LinearRegressionCombinations original, Cloner cloner)
      : base(original, cloner) {
      RegisterEventHandlers();
    }
    public override IDeepCloneable Clone(Cloner cloner) {
      return new LinearRegressionCombinations(this, cloner);
    }

    public LinearRegressionCombinations()
      : base() {
      Parameters.Add(new FixedValueParameter<IntValue>("Maximum Inputs", "The maximum number of input variables used in the linear models.", new IntValue(1)));
      Parameters.Add(new FixedValueParameter<IntValue>("Maximum Solutions stored", "The maximum number of solutions that are stored per number of inputs.", new IntValue(1000)));
      Parameters.Add(new FixedValueParameter<BoolValue>("Create Solution", "Flag that indicates if a solution should be produced at the end of the run", new BoolValue(false)));

      Problem = new RegressionProblem();
      RegisterEventHandlers();
    }

    private void RegisterEventHandlers() {
      Problem.ProblemDataChanged += (o, e) => { MaximumInputs = Problem.ProblemData.InputVariables.CheckedItems.Count(); };
    }
    protected override void OnProblemChanged() {
      base.OnProblemChanged();
      MaximumInputs = Problem.ProblemData.InputVariables.CheckedItems.Count();
    }


    private static long CalculateCombinations(int maximumInputs, int totalVariables) {
      long combinations = 0;

      for (int i = 1; i <= maximumInputs; i++)
        combinations += Common.EnumerableExtensions.BinomialCoefficient(totalVariables, i);
      return combinations;
    }

    protected override void Run(CancellationToken cancellationToken) {
      double[,] inputMatrix = Problem.ProblemData.Dataset.ToArray(Problem.ProblemData.AllowedInputVariables.Concat(new string[] { Problem.ProblemData.TargetVariable }), Problem.ProblemData.TrainingIndices);
      if (inputMatrix.Cast<double>().Any(x => double.IsNaN(x) || double.IsInfinity(x)))
        throw new NotSupportedException("Linear regression does not support NaN or infinity values in the input dataset.");

      var templateProblemData = (IRegressionProblemData)Problem.ProblemData.Clone();
      foreach (var variable in templateProblemData.InputVariables)
        templateProblemData.InputVariables.SetItemCheckedState(variable, false);

      var inputVariableNames = Problem.ProblemData.InputVariables.CheckedItems.Select(i => i.Value.Value).ToList();
      var createSolution = CreateSolution;
      var maximumInputs = MaximumInputs;
      var maximumSolutions = MaximumSolutions;

      var storedRuns = new List<IRun>[maximumInputs];
      var runs = new ConcurrentBag<IRun>();

      TotalModels = (int)CalculateCombinations(MaximumInputs, inputVariableNames.Count);
      CalculatedModels = 0;
      CalculatedInputs = 0;

      for (int inputs = 1; inputs <= MaximumInputs; inputs++) {
        Parallel.ForEach(inputVariableNames.Combinations(inputs).ToList(), inputCombination => {
          var problemData = new RegressionProblemData(templateProblemData.Dataset, inputCombination, templateProblemData.TargetVariable);
          problemData.TrainingPartition.Start = templateProblemData.TrainingPartition.Start;
          problemData.TrainingPartition.End = templateProblemData.TrainingPartition.End;
          problemData.TestPartition.Start = templateProblemData.TestPartition.Start;
          problemData.TestPartition.End = templateProblemData.TestPartition.End;

          double trainRmsError, testRmsError;
          var solution = CreateLinearRegressionSolution(problemData, createSolution, out trainRmsError, out testRmsError);

          var run = new Run();
          run.Name = string.Format("Run - Inputs {0}/{1}", inputCombination.Count(), MaximumInputs);
          if (solution != null) run.Results.Add("Solution", solution);
          run.Results.Add("RMSE train", new DoubleValue(trainRmsError));
          run.Results.Add("RMSE test", new DoubleValue(testRmsError));
          run.Results.Add("Inputs", new IntValue(inputCombination.Count()));
          run.Results.Add("Input names", new StringValue(string.Join(" ", inputCombination)));
          runs.Add(run);
        });

        CalculatedModels += runs.Count;
        CalculatedInputs = inputs;
        storedRuns[inputs - 1] = runs.OrderBy(r => ((DoubleValue)r.Results["RMSE test"]).Value).Take(maximumSolutions).ToList();
        runs = new ConcurrentBag<IRun>();

        if (cancellationToken.IsCancellationRequested) {
          Results.Add(new Result("Runs", new RunCollection(storedRuns.SelectMany(r => r))));
          cancellationToken.ThrowIfCancellationRequested();
        }
      }


      Results.Add(new Result("Runs", new RunCollection(storedRuns.SelectMany(r => r))));
    }


    public static ISymbolicRegressionSolution CreateLinearRegressionSolution(IRegressionProblemData problemData, bool buildSolution, out double trainRmsError, out double testRmsError) {
      var dataset = problemData.Dataset;
      string targetVariable = problemData.TargetVariable;
      IEnumerable<string> allowedInputVariables = problemData.AllowedInputVariables;

      double[,] inputMatrix = dataset.ToArray(allowedInputVariables.Concat(new string[] { targetVariable }), problemData.TrainingIndices);
      double[,] testMatrix = dataset.ToArray(allowedInputVariables.Concat(new string[] { targetVariable }), problemData.TestIndices);

      alglib.linearmodel lm = new alglib.linearmodel();
      alglib.lrreport ar = new alglib.lrreport();
      int nRows = inputMatrix.GetLength(0);
      int nFeatures = inputMatrix.GetLength(1) - 1;
      double[] coefficients = new double[nFeatures + 1]; // last coefficient is for the constant

      int retVal = 1;
      alglib.lrbuild(inputMatrix, nRows, nFeatures, out retVal, out lm, out ar);
      if (retVal != 1) throw new ArgumentException("Error in calculation of linear regression solution");
      trainRmsError = ar.rmserror;

      alglib.lrunpack(lm, out coefficients, out nFeatures);
      testRmsError = alglib.lrrmserror(lm, testMatrix, testMatrix.GetLength(0));

      if (!buildSolution) return null;

      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 = coefficients[col];
        addition.AddSubtree(vNode);
        col++;
      }

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

      SymbolicRegressionSolution solution = new SymbolicRegressionSolution(new SymbolicRegressionModel(problemData.TargetVariable, tree, new SymbolicDataAnalysisExpressionTreeInterpreter()), problemData);
      solution.Model.Name = "Linear Regression Model";
      solution.Name = "Linear Regression Solution";
      return solution;
    }
  }
}