Free cookie consent management tool by TermsFeed Policy Generator

source: branches/Persistence Test/HeuristicLab.LinearRegression/3.2/LinearRegressionOperator.cs @ 4593

Last change on this file since 4593 was 2445, checked in by gkronber, 15 years ago

Fixed #787 (LinearRegressionOperator uses leastsquares function of ALGLIB instead of linearregression function)

File size: 9.4 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2008 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
21using System;
22using System.Collections.Generic;
23using System.Linq;
24using System.Text;
25using HeuristicLab.Core;
26using HeuristicLab.Common;
27using HeuristicLab.Data;
28using HeuristicLab.DataAnalysis;
29using HeuristicLab.Modeling;
30using HeuristicLab.GP;
31using HeuristicLab.GP.StructureIdentification;
32using HeuristicLab.GP.Interfaces;
33
34namespace HeuristicLab.LinearRegression {
35  public class LinearRegressionOperator : OperatorBase {
36    private static double constant = 1.0;
37
38    public LinearRegressionOperator() {
39      AddVariableInfo(new VariableInfo("TargetVariable", "Name of the target variable", typeof(StringData), VariableKind.In));
40      AddVariableInfo(new VariableInfo("Dataset", "Dataset with all samples on which to apply the function", typeof(Dataset), VariableKind.In));
41      AddVariableInfo(new VariableInfo("SamplesStart", "Start index of samples in dataset to evaluate", typeof(IntData), VariableKind.In));
42      AddVariableInfo(new VariableInfo("SamplesEnd", "End index of samples in dataset to evaluate", typeof(IntData), VariableKind.In));
43      AddVariableInfo(new VariableInfo("MaxTimeOffset", "(optional) Maximal time offset for time-series prognosis", typeof(IntData), VariableKind.In));
44      AddVariableInfo(new VariableInfo("MinTimeOffset", "(optional) Minimal time offset for time-series prognosis", typeof(IntData), VariableKind.In));
45      AddVariableInfo(new VariableInfo("LinearRegressionModel", "Formula that was calculated by linear regression", typeof(IGeneticProgrammingModel), VariableKind.Out | VariableKind.New));
46    }
47
48    public override IOperation Apply(IScope scope) {
49      Dataset dataset = GetVariableValue<Dataset>("Dataset", scope, true);
50      string targetVariable = GetVariableValue<StringData>("TargetVariable", scope, true).Data;
51      int targetVariableIndex = dataset.GetVariableIndex(targetVariable);
52      int start = GetVariableValue<IntData>("SamplesStart", scope, true).Data;
53      int end = GetVariableValue<IntData>("SamplesEnd", scope, true).Data;
54      IntData maxTimeOffsetData = GetVariableValue<IntData>("MaxTimeOffset", scope, true, false);
55      int maxTimeOffset = maxTimeOffsetData == null ? 0 : maxTimeOffsetData.Data;
56      IntData minTimeOffsetData = GetVariableValue<IntData>("MinTimeOffset", scope, true, false);
57      int minTimeOffset = minTimeOffsetData == null ? 0 : minTimeOffsetData.Data;
58
59      List<int> allowedColumns = CalculateAllowedColumns(dataset, targetVariableIndex, start, end);
60      List<int> allowedRows = CalculateAllowedRows(dataset, targetVariableIndex, allowedColumns, start, end, minTimeOffset, maxTimeOffset);
61
62      double[,] inputMatrix = PrepareInputMatrix(dataset, allowedColumns, allowedRows, minTimeOffset, maxTimeOffset);
63      double[] targetVector = PrepareTargetVector(dataset, targetVariableIndex, allowedRows);
64      double[] coefficients = CalculateCoefficients(inputMatrix, targetVector);
65      IFunctionTree tree = CreateModel(coefficients, allowedColumns.Select(i => dataset.GetVariableName(i)).ToList(), minTimeOffset, maxTimeOffset);
66
67      scope.AddVariable(new HeuristicLab.Core.Variable(scope.TranslateName("LinearRegressionModel"), new GeneticProgrammingModel(tree)));
68      return null;
69    }
70
71    private IFunctionTree CreateModel(double[] coefficients, List<string> allowedVariables, int minTimeOffset, int maxTimeOffset) {
72      IFunctionTree root = new Addition().GetTreeNode();
73      IFunctionTree actNode = root;
74      int timeOffsetRange = (maxTimeOffset - minTimeOffset + 1);
75
76      Queue<IFunctionTree> nodes = new Queue<IFunctionTree>();
77      for (int i = 0; i < allowedVariables.Count; i++) {
78        for (int timeOffset = minTimeOffset; timeOffset <= maxTimeOffset; timeOffset++) {
79          var vNode = (VariableFunctionTree)new GP.StructureIdentification.Variable().GetTreeNode();
80          vNode.VariableName = allowedVariables[i];
81          vNode.Weight = coefficients[(i * timeOffsetRange) + (timeOffset - minTimeOffset)];
82          vNode.SampleOffset = timeOffset;
83          nodes.Enqueue(vNode);
84        }
85      }
86      var cNode = (ConstantFunctionTree)new Constant().GetTreeNode();
87
88      cNode.Value = coefficients[coefficients.Length - 1];
89      nodes.Enqueue(cNode);
90
91      IFunctionTree newTree;
92      while (nodes.Count != 1) {
93        newTree = new Addition().GetTreeNode();
94        newTree.AddSubTree(nodes.Dequeue());
95        newTree.AddSubTree(nodes.Dequeue());
96        nodes.Enqueue(newTree);
97      }
98
99      return nodes.Dequeue();
100    }
101
102    private double[] CalculateCoefficients(double[,] inputMatrix, double[] targetVector) {
103      int retVal = 0;
104      alglib.linreg.linearmodel lm = new alglib.linreg.linearmodel();
105      alglib.linreg.lrreport ar = new alglib.linreg.lrreport();
106      int n = targetVector.Length;
107      int p = inputMatrix.GetLength(1);
108      double[,] dataset = new double[n, p];
109      for (int row = 0; row < n; row++) {
110        for (int column = 0; column < p-1; column++) {
111          dataset[row, column] = inputMatrix[row, column];
112        }
113        dataset[row, p-1] = targetVector[row];
114      }
115      alglib.linreg.lrbuild(ref dataset, n, p-1, ref retVal, ref lm, ref ar);
116      if (retVal != 1) throw new ArgumentException("Error in calculation of linear regression model");
117      Console.Out.WriteLine("ALGLIB Linear Regression: Estimated generalization RMS = {0}", ar.cvrmserror);
118
119      double[] coefficients = new double[p];
120      for (int i = 0; i < p; i++) {
121        coefficients[i] = lm.w[i+4];
122      }
123      return coefficients;
124    }
125
126    //returns list of valid row indexes (rows without NaN values)
127    private List<int> CalculateAllowedRows(Dataset dataset, int targetVariable, IList<int> allowedColumns, int start, int end, int minTimeOffset, int maxTimeOffset) {
128      List<int> allowedRows = new List<int>();
129      bool add;
130      for (int row = start; row < end; row++) {
131        add = true;
132        for (int colIndex = 0; colIndex < allowedColumns.Count && add == true; colIndex++) {
133          for (int timeOffset = minTimeOffset; timeOffset <= maxTimeOffset; timeOffset++) {
134            if (
135              row + timeOffset < 0 ||
136              row + timeOffset > dataset.Rows ||
137              double.IsNaN(dataset.GetValue(row + timeOffset, allowedColumns[colIndex])) ||
138              double.IsInfinity(dataset.GetValue(row + timeOffset, allowedColumns[colIndex])) ||
139              double.IsNaN(dataset.GetValue(row + timeOffset, targetVariable))) {
140              add = false;
141            }
142          }
143        }
144        if (add)
145          allowedRows.Add(row);
146        add = true;
147      }
148      return allowedRows;
149    }
150
151    //returns list of valid column indexes (columns which contain max. 10% NaN (or infinity) and contain at least two different values)
152    private List<int> CalculateAllowedColumns(Dataset dataset, int targetVariable, int start, int end) {
153      List<int> allowedColumns = new List<int>();
154      double n = end - start;
155      for (int i = 0; i < dataset.Columns; i++) {
156        double nanRatio = dataset.CountMissingValues(i, start, end) / n;
157        if (i != targetVariable && nanRatio < 0.1 && dataset.GetRange(i, start, end) > 0.0) {
158          allowedColumns.Add(i);
159        }
160      }
161      return allowedColumns;
162    }
163
164    private double[,] PrepareInputMatrix(Dataset dataset, List<int> allowedColumns, List<int> allowedRows, int minTimeOffset, int maxTimeOffset) {
165      int rowCount = allowedRows.Count;
166      int timeOffsetRange = (maxTimeOffset - minTimeOffset + 1);
167      double[,] matrix = new double[rowCount, (allowedColumns.Count * timeOffsetRange) + 1];
168      for (int row = 0; row < allowedRows.Count; row++)
169        for (int col = 0; col < allowedColumns.Count; col++) {
170          for (int timeOffset = minTimeOffset; timeOffset <= maxTimeOffset; timeOffset++)
171            matrix[row, (col * timeOffsetRange) + (timeOffset - minTimeOffset)] = dataset.GetValue(allowedRows[row] + timeOffset, allowedColumns[col]);
172        }
173      //add constant 1.0 in last column
174      for (int i = 0; i < rowCount; i++)
175        matrix[i, allowedColumns.Count * timeOffsetRange] = constant;
176      return matrix;
177    }
178
179    private double[] PrepareTargetVector(Dataset dataset, int targetVariable, List<int> allowedRows) {
180      int rowCount = allowedRows.Count;
181      double[] targetVector = new double[rowCount];
182      double[] samples = dataset.Samples;
183      for (int row = 0; row < rowCount; row++) {
184        targetVector[row] = dataset.GetValue(allowedRows[row], targetVariable);
185      }
186      return targetVector;
187    }
188  }
189}
Note: See TracBrowser for help on using the repository browser.