source: branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs @ 15433

#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.Analysis;
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 {
  // UNFINISHED
  [Item("Generalized Additive Modelling", "GAM")]
  [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 102)]
  [StorableClass]
  public sealed class GAM : FixedDataAnalysisAlgorithm<IRegressionProblem> {
    [StorableConstructor]
    private GAM(bool deserializing) : base(deserializing) { }
    [StorableHook(HookType.AfterDeserialization)]
    private void AfterDeserialization() {      
    }

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

    public GAM()
      : base() {
      Problem = new RegressionProblem();
      Parameters.Add(new ValueParameter<DoubleValue>("Lambda", "Regularization for smoothing splines", new DoubleValue(1.0)));
      Parameters.Add(new ValueParameter<IntValue>("Max iterations", "", new IntValue(100)));
      Parameters.Add(new ValueParameter<IntValue>("Max interactions", "", new IntValue(1)));
    }    

    protected override void Run(CancellationToken cancellationToken) {

      double lambda = ((IValueParameter<DoubleValue>)Parameters["Lambda"]).Value.Value;
      int maxIters = ((IValueParameter<IntValue>)Parameters["Max iterations"]).Value.Value;
      int maxInteractions = ((IValueParameter<IntValue>)Parameters["Max interactions"]).Value.Value;
      if (maxInteractions < 1 || maxInteractions > 5) throw new ArgumentException("Max interactions is outside the valid range [1 .. 5]");

      // calculates a GAM model using a linear representation + independent non-linear functions of each variable
      // using backfitting algorithm (see The Elements of Statistical Learning page 298)

      var problemData = Problem.ProblemData;
      var y = problemData.TargetVariableTrainingValues.ToArray();
      var avgY = y.Average();
      var inputVars = Problem.ProblemData.AllowedInputVariables.ToArray();
      var nTerms = inputVars.Length; // LR
      for(int i=1;i<=maxInteractions;i++) {
        nTerms += inputVars.Combinations(i).Count();
      }
      IRegressionModel[] f = new IRegressionModel[nTerms]; 
      for(int i=0;i<f.Length;i++) {
        f[i] = new ConstantModel(0.0, problemData.TargetVariable);
      }

      var rmseTable = new DataTable("RMSE");
      var rmseRow = new DataRow("RMSE (train)");
      var rmseRowTest = new DataRow("RMSE (test)");
      rmseTable.Rows.Add(rmseRow);
      rmseTable.Rows.Add(rmseRowTest);

      Results.Add(new Result("RMSE", rmseTable));
      rmseRow.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TrainingIndices).StandardDeviation()); // -1 index to use all predictors
      rmseRowTest.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TestIndices).StandardDeviation());

      // for analytics
      double[] rss = new double[f.Length];
      string[] terms = new string[f.Length];
      Results.Add(new Result("RSS Values", typeof(DoubleMatrix)));

      // until convergence
      int iters = 0;
      var t = new double[y.Length];
      while (iters++ < maxIters) {
        int j = 0;
        foreach (var inputVar in inputVars) {
          var res = CalculateResiduals(problemData, f, j, avgY, problemData.TrainingIndices);
          rss[j] = res.Variance();
          terms[j] = inputVar;
          f[j] = RegressLR(problemData, inputVar, res);
          j++;
        }

        for(int interaction = 1; interaction <= maxInteractions;interaction++) {
          var selectedVars = HeuristicLab.Common.EnumerableExtensions.Combinations(inputVars, interaction);

          foreach (var element in selectedVars) {
            var res = CalculateResiduals(problemData, f, j, avgY, problemData.TrainingIndices);
            rss[j] = res.Variance();
            terms[j] = string.Format("f({0})", string.Join(",", element));
            f[j] = RegressSpline(problemData, element.ToArray(), res, lambda);
            j++;
          }
        }

        rmseRow.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TrainingIndices).StandardDeviation()); // -1 index to use all predictors
        rmseRowTest.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TestIndices).StandardDeviation());

        // calculate table with residual contributions of each term
        var rssTable = new DoubleMatrix(rss.Length, 1, new string[] { "RSS" }, terms);
        for (int i = 0; i < rss.Length; i++) rssTable[i, 0] = rss[i];
        Results["RSS Values"].Value = rssTable;

        if (cancellationToken.IsCancellationRequested) break;
      }

      var model = new RegressionEnsembleModel(f.Concat(new[] { new ConstantModel(avgY, problemData.TargetVariable) }));
      model.AverageModelEstimates = false;
      Results.Add(new Result("Ensemble solution", model.CreateRegressionSolution((IRegressionProblemData)problemData.Clone())));

    }

    private double[] CalculateResiduals(IRegressionProblemData problemData, IRegressionModel[] f, int j, double avgY, IEnumerable<int> rows) {
      var y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, rows);
      double[] t = y.Select(yi => yi - avgY).ToArray();
      // collect other predictions
      for (int k = 0; k < f.Length; k++) {
        if (k != j) {
          var pred = f[k].GetEstimatedValues(problemData.Dataset, rows).ToArray();
          // determine target for this smoother
          for (int i = 0; i < t.Length; i++) {
            t[i] -= pred[i];
          }
        }
      }
      return t;
    }

    private IRegressionModel RegressLR(IRegressionProblemData problemData, string inputVar, double[] target) {
      // Umständlich!
      var ds = ((Dataset)problemData.Dataset).ToModifiable();
      ds.ReplaceVariable(problemData.TargetVariable, target.Concat(Enumerable.Repeat(0.0, ds.Rows - target.Length)).ToList<double>());
      var pd = new RegressionProblemData(ds, new string[] { inputVar }, problemData.TargetVariable);
      pd.TrainingPartition.Start = problemData.TrainingPartition.Start;
      pd.TrainingPartition.End = problemData.TrainingPartition.End;
      pd.TestPartition.Start = problemData.TestPartition.Start;
      pd.TestPartition.End = problemData.TestPartition.End;
      double rmsError, cvRmsError;
      return LinearRegression.CreateLinearRegressionSolution(pd, out rmsError, out cvRmsError).Model;
    }

    private IRegressionModel RegressSpline(IRegressionProblemData problemData, string inputVar, double[] target, double lambda) {
      if (problemData.Dataset.VariableHasType<double>(inputVar)) {
        // Umständlich!
        return Splines.CalculatePenalizedRegressionSpline(
          problemData.Dataset.GetDoubleValues(inputVar, problemData.TrainingIndices).ToArray(),
          (double[])target.Clone(), lambda,
          problemData.TargetVariable, new string[] { inputVar }
          );
      } else return new ConstantModel(target.Average(), problemData.TargetVariable);
    }
    private IRegressionModel RegressSpline(IRegressionProblemData problemData, string[] inputVars, double[] target, double lambda) {
      if (inputVars.All(problemData.Dataset.VariableHasType<double>)) {
        var product = problemData.Dataset.GetDoubleValues(inputVars.First(), problemData.TrainingIndices).ToArray();
        for(int i = 1;i<inputVars.Length;i++) {
          product = product.Zip(problemData.Dataset.GetDoubleValues(inputVars[i], problemData.TrainingIndices), (pi, vi) => pi * vi).ToArray();
        }
        // Umständlich!
        return Splines.CalculatePenalizedRegressionSpline(
          product,
          (double[])target.Clone(), lambda,
          problemData.TargetVariable, inputVars
          );
      } else return new ConstantModel(target.Average(), problemData.TargetVariable);
    }

    private IRegressionModel RegressRF(IRegressionProblemData problemData, string inputVar, double[] target, double lambda) {
      if (problemData.Dataset.VariableHasType<double>(inputVar)) {
        // Umständlich!
        var ds = ((Dataset)problemData.Dataset).ToModifiable();
        ds.ReplaceVariable(problemData.TargetVariable, target.Concat(Enumerable.Repeat(0.0, ds.Rows - target.Length)).ToList<double>());
        var pd = new RegressionProblemData(ds, new string[] { inputVar }, problemData.TargetVariable);
        pd.TrainingPartition.Start = problemData.TrainingPartition.Start;
        pd.TrainingPartition.End = problemData.TrainingPartition.End;
        pd.TestPartition.Start = problemData.TestPartition.Start;
        pd.TestPartition.End = problemData.TestPartition.End;
        double rmsError, oobRmsError;
        double avgRelError, oobAvgRelError; 
        return RandomForestRegression.CreateRandomForestRegressionModel(pd, 100, 0.5, 0.5, 1234, out rmsError, out avgRelError, out oobRmsError, out oobAvgRelError);
      } else return new ConstantModel(target.Average(), problemData.TargetVariable);
    }
  }


  // UNFINISHED
  public class RBFModel : NamedItem, IRegressionModel {
    private alglib.rbfmodel model;

    public string TargetVariable { get; set; }

    public IEnumerable<string> VariablesUsedForPrediction { get; private set; }
    private ITransformation<double>[] scaling;

    public event EventHandler TargetVariableChanged;

    public RBFModel(RBFModel orig, Cloner cloner) : base(orig, cloner) {
      this.TargetVariable = orig.TargetVariable;
      this.VariablesUsedForPrediction = orig.VariablesUsedForPrediction.ToArray();
      this.model = (alglib.rbfmodel)orig.model.make_copy();
      this.scaling = orig.scaling.Select(s => cloner.Clone(s)).ToArray();
    }
    public RBFModel(alglib.rbfmodel model, string targetVar, string[] inputs, IEnumerable<ITransformation<double>> scaling) : base("RBFModel", "RBFModel") {
      this.model = model;
      this.TargetVariable = targetVar;
      this.VariablesUsedForPrediction = inputs;
      this.scaling = scaling.ToArray();
    }

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

    public IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
      return new RegressionSolution(this, (IRegressionProblemData)problemData.Clone());
    }

    public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
      double[] x = new double[VariablesUsedForPrediction.Count()];
      double[] y;
      foreach (var r in rows) {
        int c = 0;
        foreach (var v in VariablesUsedForPrediction) {
          x[c] = scaling[c].Apply(dataset.GetDoubleValue(v, r).ToEnumerable()).First(); // OUCH!
          c++;
        }
        alglib.rbfcalc(model, x, out y);
        yield return y[0];
      }
    }
  }
}