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

Last change on this file since 15449 was 15449, checked in by gkronber, 4 years ago

#2789 created x64 build of CUBGCV algorithm, fixed some bugs

File size: 13.3 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using System;
23using System.Collections.Concurrent;
24using System.Collections.Generic;
25using System.Linq;
26using System.Threading;
27using System.Threading.Tasks;
28using HeuristicLab.Analysis;
29using HeuristicLab.Common;
30using HeuristicLab.Core;
31using HeuristicLab.Data;
32using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
33using HeuristicLab.Optimization;
34using HeuristicLab.Parameters;
35using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
36using HeuristicLab.Problems.DataAnalysis;
37using HeuristicLab.Problems.DataAnalysis.Symbolic;
38using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;
39
40namespace HeuristicLab.Algorithms.DataAnalysis.Experimental {
41  // UNFINISHED
42  [Item("Generalized Additive Modelling", "GAM")]
43  [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 102)]
44  [StorableClass]
45  public sealed class GAM : FixedDataAnalysisAlgorithm<IRegressionProblem> {
46
47    private const string LambdaParameterName = "Lambda";
48    private const string MaxIterationsParameterName = "Max iterations";
49    private const string MaxInteractionsParameterName = "Max interactions";
50
51    public IFixedValueParameter<DoubleValue> LambdaParameter {
52      get { return (IFixedValueParameter<DoubleValue>)Parameters[LambdaParameterName]; }
53    }
54    public IFixedValueParameter<IntValue> MaxIterationsParameter {
55      get { return (IFixedValueParameter<IntValue>)Parameters[MaxIterationsParameterName]; }
56    }
57    public IFixedValueParameter<IntValue> MaxInteractionsParameter {
58      get { return (IFixedValueParameter<IntValue>)Parameters[MaxInteractionsParameterName]; }
59    }
60
61    public double Lambda {
62      get { return LambdaParameter.Value.Value; }
63      set { LambdaParameter.Value.Value = value; }
64    }
65    public int MaxIterations {
66      get { return MaxIterationsParameter.Value.Value; }
67      set { MaxIterationsParameter.Value.Value = value; }
68    }
69    public int MaxInteractions {
70      get { return MaxInteractionsParameter.Value.Value; }
71      set { MaxInteractionsParameter.Value.Value = value; }
72    }
73
74    [StorableConstructor]
75    private GAM(bool deserializing) : base(deserializing) { }
76    [StorableHook(HookType.AfterDeserialization)]
77    private void AfterDeserialization() {
78    }
79
80    private GAM(GAM original, Cloner cloner)
81      : base(original, cloner) {
82    }
83    public override IDeepCloneable Clone(Cloner cloner) {
84      return new GAM(this, cloner);
85    }
86
87    public GAM()
88      : base() {
89      Problem = new RegressionProblem();
90      Parameters.Add(new FixedValueParameter<DoubleValue>(LambdaParameterName, "Regularization for smoothing splines", new DoubleValue(1.0)));
91      Parameters.Add(new FixedValueParameter<IntValue>(MaxIterationsParameterName, "", new IntValue(100)));
92      Parameters.Add(new FixedValueParameter<IntValue>(MaxInteractionsParameterName, "", new IntValue(1)));
93    }
94
95
96    protected override void Run(CancellationToken cancellationToken) {
97      double lambda = Lambda;
98      int maxIters = MaxIterations ;
99      int maxInteractions = MaxInteractions;
100      if (maxInteractions < 1 || maxInteractions > 5) throw new ArgumentException("Max interactions is outside the valid range [1 .. 5]");
101
102      // calculates a GAM model using a linear representation + independent non-linear functions of each variable
103      // using backfitting algorithm (see The Elements of Statistical Learning page 298)
104
105      var problemData = Problem.ProblemData;
106      var y = problemData.TargetVariableTrainingValues.ToArray();
107      var avgY = y.Average();
108      var inputVars = Problem.ProblemData.AllowedInputVariables.ToArray();
109      var nTerms = 0; // inputVars.Length; // LR
110      for (int i = 1; i <= maxInteractions; i++) {
111        nTerms += inputVars.Combinations(i).Count();
112      }
113      IRegressionModel[] f = new IRegressionModel[nTerms];
114      for (int i = 0; i < f.Length; i++) {
115        f[i] = new ConstantModel(0.0, problemData.TargetVariable);
116      }
117
118      var rmseTable = new DataTable("RMSE");
119      var rmseRow = new DataRow("RMSE (train)");
120      var rmseRowTest = new DataRow("RMSE (test)");
121      rmseTable.Rows.Add(rmseRow);
122      rmseTable.Rows.Add(rmseRowTest);
123
124      Results.Add(new Result("RMSE", rmseTable));
125      rmseRow.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TrainingIndices).StandardDeviation()); // -1 index to use all predictors
126      rmseRowTest.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TestIndices).StandardDeviation());
127
128      // for analytics
129      double[] rss = new double[f.Length];
130      string[] terms = new string[f.Length];
131      Results.Add(new Result("RSS Values", typeof(DoubleMatrix)));
132
133      // until convergence
134      int iters = 0;
135      var t = new double[y.Length];
136      while (iters++ < maxIters) {
137        int j = 0;
138        //foreach (var inputVar in inputVars) {
139        //  var res = CalculateResiduals(problemData, f, j, avgY, problemData.TrainingIndices);
140        //  rss[j] = res.Variance();
141        //  terms[j] = inputVar;
142        //  f[j] = RegressLR(problemData, inputVar, res);
143        //  j++;
144        //}
145
146        for (int interaction = 1; interaction <= maxInteractions; interaction++) {
147          var selectedVars = HeuristicLab.Common.EnumerableExtensions.Combinations(inputVars, interaction);
148
149          foreach (var element in selectedVars) {
150            var res = CalculateResiduals(problemData, f, j, avgY, problemData.TrainingIndices);
151            rss[j] = res.Variance();
152            terms[j] = string.Format("f({0})", string.Join(",", element));
153            f[j] = RegressSpline(problemData, element.ToArray(), res, lambda);
154            j++;
155          }
156        }
157
158        rmseRow.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TrainingIndices).StandardDeviation()); // -1 index to use all predictors
159        rmseRowTest.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TestIndices).StandardDeviation());
160
161        // calculate table with residual contributions of each term
162        var rssTable = new DoubleMatrix(rss.Length, 1, new string[] { "RSS" }, terms);
163        for (int i = 0; i < rss.Length; i++) rssTable[i, 0] = rss[i];
164        Results["RSS Values"].Value = rssTable;
165
166        if (cancellationToken.IsCancellationRequested) break;
167      }
168
169      var model = new RegressionEnsembleModel(f.Concat(new[] { new ConstantModel(avgY, problemData.TargetVariable) }));
170      model.AverageModelEstimates = false;
171      var solution = model.CreateRegressionSolution((IRegressionProblemData)problemData.Clone());
172      Results.Add(new Result("Ensemble solution", solution));
173    }
174
175    private double[] CalculateResiduals(IRegressionProblemData problemData, IRegressionModel[] f, int j, double avgY, IEnumerable<int> rows) {
176      var y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, rows);
177      double[] t = y.Select(yi => yi - avgY).ToArray();
178      // collect other predictions
179      for (int k = 0; k < f.Length; k++) {
180        if (k != j) {
181          var pred = f[k].GetEstimatedValues(problemData.Dataset, rows).ToArray();
182          // determine target for this smoother
183          for (int i = 0; i < t.Length; i++) {
184            t[i] -= pred[i];
185          }
186        }
187      }
188      return t;
189    }
190
191    private IRegressionModel RegressLR(IRegressionProblemData problemData, string inputVar, double[] target) {
192      // Umständlich!
193      var ds = ((Dataset)problemData.Dataset).ToModifiable();
194      ds.ReplaceVariable(problemData.TargetVariable, target.Concat(Enumerable.Repeat(0.0, ds.Rows - target.Length)).ToList<double>());
195      var pd = new RegressionProblemData(ds, new string[] { inputVar }, problemData.TargetVariable);
196      pd.TrainingPartition.Start = problemData.TrainingPartition.Start;
197      pd.TrainingPartition.End = problemData.TrainingPartition.End;
198      pd.TestPartition.Start = problemData.TestPartition.Start;
199      pd.TestPartition.End = problemData.TestPartition.End;
200      double rmsError, cvRmsError;
201      return LinearRegression.CreateLinearRegressionSolution(pd, out rmsError, out cvRmsError).Model;
202    }
203
204    // private IRegressionModel RegressSpline(IRegressionProblemData problemData, string inputVar, double[] target, double lambda) {
205    //   if (problemData.Dataset.VariableHasType<double>(inputVar)) {
206    //     // Umständlich!
207    //     return Splines.CalculatePenalizedRegressionSpline(
208    //       problemData.Dataset.GetDoubleValues(inputVar, problemData.TrainingIndices).ToArray(),
209    //       (double[])target.Clone(), lambda,
210    //       problemData.TargetVariable, new string[] { inputVar }
211    //       );
212    //   } else return new ConstantModel(target.Average(), problemData.TargetVariable);
213    // }
214    private IRegressionModel RegressSpline(IRegressionProblemData problemData, string[] inputVars, double[] target, double lambda) {
215      if (inputVars.All(problemData.Dataset.VariableHasType<double>)) {
216        var product = problemData.Dataset.GetDoubleValues(inputVars.First(), problemData.TrainingIndices).ToArray();
217        for (int i = 1; i < inputVars.Length; i++) {
218          product = product.Zip(problemData.Dataset.GetDoubleValues(inputVars[i], problemData.TrainingIndices), (pi, vi) => pi * vi).ToArray();
219        }
220        CubicSplineGCV.CubGcvReport report;
221        return CubicSplineGCV.CalculateCubicSpline(
222          product,
223          (double[])target.Clone(),
224          problemData.TargetVariable, inputVars, out report
225          );
226      } else return new ConstantModel(target.Average(), problemData.TargetVariable);
227    }
228
229    private IRegressionModel RegressRF(IRegressionProblemData problemData, string inputVar, double[] target, double lambda) {
230      if (problemData.Dataset.VariableHasType<double>(inputVar)) {
231        // Umständlich!
232        var ds = ((Dataset)problemData.Dataset).ToModifiable();
233        ds.ReplaceVariable(problemData.TargetVariable, target.Concat(Enumerable.Repeat(0.0, ds.Rows - target.Length)).ToList<double>());
234        var pd = new RegressionProblemData(ds, new string[] { inputVar }, problemData.TargetVariable);
235        pd.TrainingPartition.Start = problemData.TrainingPartition.Start;
236        pd.TrainingPartition.End = problemData.TrainingPartition.End;
237        pd.TestPartition.Start = problemData.TestPartition.Start;
238        pd.TestPartition.End = problemData.TestPartition.End;
239        double rmsError, oobRmsError;
240        double avgRelError, oobAvgRelError;
241        return RandomForestRegression.CreateRandomForestRegressionModel(pd, 100, 0.5, 0.5, 1234, out rmsError, out avgRelError, out oobRmsError, out oobAvgRelError);
242      } else return new ConstantModel(target.Average(), problemData.TargetVariable);
243    }
244  }
245
246
247  // UNFINISHED
248  public class RBFModel : NamedItem, IRegressionModel {
249    private alglib.rbfmodel model;
250
251    public string TargetVariable { get; set; }
252
253    public IEnumerable<string> VariablesUsedForPrediction { get; private set; }
254    private ITransformation<double>[] scaling;
255
256    public event EventHandler TargetVariableChanged;
257
258    public RBFModel(RBFModel orig, Cloner cloner) : base(orig, cloner) {
259      this.TargetVariable = orig.TargetVariable;
260      this.VariablesUsedForPrediction = orig.VariablesUsedForPrediction.ToArray();
261      this.model = (alglib.rbfmodel)orig.model.make_copy();
262      this.scaling = orig.scaling.Select(s => cloner.Clone(s)).ToArray();
263    }
264    public RBFModel(alglib.rbfmodel model, string targetVar, string[] inputs, IEnumerable<ITransformation<double>> scaling) : base("RBFModel", "RBFModel") {
265      this.model = model;
266      this.TargetVariable = targetVar;
267      this.VariablesUsedForPrediction = inputs;
268      this.scaling = scaling.ToArray();
269    }
270
271    public override IDeepCloneable Clone(Cloner cloner) {
272      return new RBFModel(this, cloner);
273    }
274
275    public IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
276      return new RegressionSolution(this, (IRegressionProblemData)problemData.Clone());
277    }
278
279    public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
280      double[] x = new double[VariablesUsedForPrediction.Count()];
281      double[] y;
282      foreach (var r in rows) {
283        int c = 0;
284        foreach (var v in VariablesUsedForPrediction) {
285          x[c] = scaling[c].Apply(dataset.GetDoubleValue(v, r).ToEnumerable()).First(); // OUCH!
286          c++;
287        }
288        alglib.rbfcalc(model, x, out y);
289        yield return y[0];
290      }
291    }
292  }
293}
Note: See TracBrowser for help on using the repository browser.