#region License Information
/* HeuristicLab
* Copyright (C) 2002-2015 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
* and the BEACON Center for the Study of Evolution in Action.
*
* 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 .
*/
#endregion
using System;
using System.Collections;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using HeuristicLab.Core;
using HeuristicLab.Problems.DataAnalysis;
namespace HeuristicLab.Algorithms.DataAnalysis {
// This class implements a greedy decision tree learner which selects splits with the maximum reduction in sum of squared errors.
// The tree builder also tracks variable relevance metrics based on the splits and improvement after the split.
// The implementation is tuned for gradient boosting where multiple trees have to be calculated for the same training data
// each time with a different target vector. Vectors of idx to allow iteration of intput variables in sorted order are
// pre-calculated so that optimal thresholds for splits can be calculated in O(n) for each input variable.
// After each split the row idx are partitioned in a left an right part.
public class RegressionTreeBuilder {
private readonly IRandom random;
private readonly IRegressionProblemData problemData;
private readonly int nCols;
private readonly double[][] x; // all training data (original order from problemData), x is constant
private double[] y; // training labels (original order from problemData), y can be changed
private Dictionary sumImprovements; // for variable relevance calculation
private readonly string[] allowedVariables; // all variables in shuffled order
private Dictionary varName2Index; // maps the variable names to column indexes
private int effectiveVars; // number of variables that are used from allowedVariables
private int effectiveRows; // number of rows that are used from
private readonly int[][] sortedIdxAll;
private readonly int[][] sortedIdx; // random selection from sortedIdxAll (for r < 1.0)
private int calls = 0;
// helper arrays which are allocated to maximal necessary size only once in the ctor
private readonly int[] internalIdx, which, leftTmp, rightTmp;
private readonly double[] outx;
private readonly int[] outSortedIdx;
private RegressionTreeModel.TreeNode[] tree; // tree is represented as a flat array of nodes
private int curTreeNodeIdx; // the index where the next tree node is stored
private class Partition {
public int ParentNodeIdx { get; set; }
public int Depth { get; set; }
public int StartIdx { get; set; }
public int EndIndex { get; set; }
public bool Left { get; set; }
}
private readonly SortedList queue;
// prepare and allocate buffer variables in ctor
public RegressionTreeBuilder(IRegressionProblemData problemData, IRandom random) {
this.problemData = problemData;
this.random = random;
var rows = problemData.TrainingIndices.Count();
this.nCols = problemData.AllowedInputVariables.Count();
allowedVariables = problemData.AllowedInputVariables.ToArray();
varName2Index = new Dictionary(allowedVariables.Length);
for (int i = 0; i < allowedVariables.Length; i++) varName2Index.Add(allowedVariables[i], i);
sortedIdxAll = new int[nCols][];
sortedIdx = new int[nCols][];
sumImprovements = new Dictionary();
internalIdx = new int[rows];
which = new int[rows];
leftTmp = new int[rows];
rightTmp = new int[rows];
outx = new double[rows];
outSortedIdx = new int[rows];
queue = new SortedList();
x = new double[nCols][];
y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TrainingIndices).ToArray();
int col = 0;
foreach (var inputVariable in problemData.AllowedInputVariables) {
x[col] = problemData.Dataset.GetDoubleValues(inputVariable, problemData.TrainingIndices).ToArray();
sortedIdxAll[col] = Enumerable.Range(0, rows).OrderBy(r => x[col][r]).ToArray();
sortedIdx[col] = new int[rows];
col++;
}
}
// simple API produces a single regression tree optimizing sum of squared errors
// this can be used if only a simple regression tree should be produced
// for a set of trees use the method CreateRegressionTreeForGradientBoosting below
//
// r and m work in the same way as for alglib random forest
// r is fraction of rows to use for training
// m is fraction of variables to use for training
public IRegressionModel CreateRegressionTree(int maxDepth, double r = 0.5, double m = 0.5) {
// subtract mean of y first
var yAvg = y.Average();
for (int i = 0; i < y.Length; i++) y[i] -= yAvg;
var seLoss = new SquaredErrorLoss();
var zeros = Enumerable.Repeat(0.0, y.Length);
var ones = Enumerable.Repeat(1.0, y.Length);
var model = CreateRegressionTreeForGradientBoosting(y, maxDepth, problemData.TrainingIndices.ToArray(), seLoss.GetLineSearchFunc(y, zeros, ones), r, m);
return new GradientBoostedTreesModel(new[] { new ConstantRegressionModel(yAvg), model }, new[] { 1.0, 1.0 });
}
// specific interface that allows to specify the target labels and the training rows which is necessary when for gradient boosted trees
public IRegressionModel CreateRegressionTreeForGradientBoosting(double[] y, int maxDepth, int[] idx, LineSearchFunc lineSearch, double r = 0.5, double m = 0.5) {
Debug.Assert(maxDepth > 0);
Debug.Assert(r > 0);
Debug.Assert(r <= 1.0);
Debug.Assert(y.Count() == this.y.Length);
Debug.Assert(m > 0);
Debug.Assert(m <= 1.0);
this.y = y; // y is changed in gradient boosting
// shuffle row idx
HeuristicLab.Random.ListExtensions.ShuffleInPlace(idx, random);
int nRows = idx.Count();
// shuffle variable idx
HeuristicLab.Random.ListExtensions.ShuffleInPlace(allowedVariables, random);
// only select a part of the rows and columns randomly
effectiveRows = (int)Math.Ceiling(nRows * r);
effectiveVars = (int)Math.Ceiling(nCols * m);
// the which array is used for partining row idxs
Array.Clear(which, 0, which.Length);
// mark selected rows
for (int row = 0; row < effectiveRows; row++) {
which[idx[row]] = 1; // we use the which vector as a temporary variable here
internalIdx[row] = idx[row];
}
for (int col = 0; col < nCols; col++) {
int i = 0;
for (int row = 0; row < nRows; row++) {
if (which[sortedIdxAll[col][row]] > 0) {
Debug.Assert(i < effectiveRows);
sortedIdx[col][i] = sortedIdxAll[col][row];
i++;
}
}
}
// prepare array for the tree nodes (a tree of maxDepth=1 has 1 node, a tree of maxDepth=d has 2^d - 1 nodes)
int numNodes = (int)Math.Pow(2, maxDepth) - 1;
this.tree = new RegressionTreeModel.TreeNode[numNodes];
this.curTreeNodeIdx = 0;
// start and end idx are inclusive
queue.Add(calls++, new Partition() { ParentNodeIdx = -1, Depth = maxDepth, StartIdx = 0, EndIndex = effectiveRows - 1 });
CreateRegressionTreeForIdx(lineSearch);
return new RegressionTreeModel(tree);
}
private void CreateRegressionTreeForIdx(LineSearchFunc lineSearch) {
while (queue.Any()) {
var f = queue.First().Value; // actually a stack
queue.RemoveAt(0);
var depth = f.Depth;
var startIdx = f.StartIdx;
var endIdx = f.EndIndex;
Debug.Assert(endIdx - startIdx >= 0);
Debug.Assert(startIdx >= 0);
Debug.Assert(endIdx < internalIdx.Length);
double threshold;
string bestVariableName;
// stop when only one row is left or no split is possible
if (depth <= 1 || endIdx - startIdx == 0 || !FindBestVariableAndThreshold(startIdx, endIdx, out threshold, out bestVariableName)) {
CreateLeafNode(startIdx, endIdx, lineSearch);
if (f.ParentNodeIdx >= 0) if (f.Left) {
tree[f.ParentNodeIdx].leftIdx = curTreeNodeIdx;
} else {
tree[f.ParentNodeIdx].rightIdx = curTreeNodeIdx;
}
curTreeNodeIdx++;
} else {
int splitIdx;
CreateInternalNode(f.StartIdx, f.EndIndex, bestVariableName, threshold, out splitIdx);
// connect to parent tree
if (f.ParentNodeIdx >= 0) if (f.Left) {
tree[f.ParentNodeIdx].leftIdx = curTreeNodeIdx;
} else {
tree[f.ParentNodeIdx].rightIdx = curTreeNodeIdx;
}
Debug.Assert(splitIdx + 1 <= endIdx);
Debug.Assert(startIdx <= splitIdx);
queue.Add(calls++, new Partition() { ParentNodeIdx = curTreeNodeIdx, Left = true, Depth = depth - 1, StartIdx = startIdx, EndIndex = splitIdx }); // left part before right part (stack organization)
queue.Add(calls++, new Partition() { ParentNodeIdx = curTreeNodeIdx, Left = false, Depth = depth - 1, StartIdx = splitIdx + 1, EndIndex = endIdx });
curTreeNodeIdx++;
}
}
}
private void CreateLeafNode(int startIdx, int endIdx, LineSearchFunc lineSearch) {
// max depth reached or only one element
tree[curTreeNodeIdx].varName = RegressionTreeModel.TreeNode.NO_VARIABLE;
tree[curTreeNodeIdx].val = lineSearch(internalIdx, startIdx, endIdx);
}
// routine for building the tree for the row idx stored in internalIdx between startIdx and endIdx
// the lineSearch function calculates the optimal prediction value for tree leaf nodes
// (in the case of squared errors it is the average of target values for the rows represented by the node)
// startIdx and endIdx are inclusive
private void CreateInternalNode(int startIdx, int endIdx, string splittingVar, double threshold, out int splitIdx) {
int bestVarIdx = varName2Index[splittingVar];
// split - two pass
// store which index goes where
for (int k = startIdx; k <= endIdx; k++) {
if (x[bestVarIdx][internalIdx[k]] <= threshold)
which[internalIdx[k]] = -1; // left partition
else
which[internalIdx[k]] = 1; // right partition
}
// partition sortedIdx for each variable
int i;
int j;
for (int col = 0; col < nCols; col++) {
i = 0;
j = 0;
int k;
for (k = startIdx; k <= endIdx; k++) {
Debug.Assert(Math.Abs(which[sortedIdx[col][k]]) == 1);
if (which[sortedIdx[col][k]] < 0) {
leftTmp[i++] = sortedIdx[col][k];
} else {
rightTmp[j++] = sortedIdx[col][k];
}
}
Debug.Assert(i > 0); // at least on element in the left partition
Debug.Assert(j > 0); // at least one element in the right partition
Debug.Assert(i + j == endIdx - startIdx + 1);
k = startIdx;
for (int l = 0; l < i; l++) sortedIdx[col][k++] = leftTmp[l];
for (int l = 0; l < j; l++) sortedIdx[col][k++] = rightTmp[l];
}
// partition row indices
i = startIdx;
j = endIdx;
while (i <= j) {
Debug.Assert(Math.Abs(which[internalIdx[i]]) == 1);
Debug.Assert(Math.Abs(which[internalIdx[j]]) == 1);
if (which[internalIdx[i]] < 0) i++;
else if (which[internalIdx[j]] > 0) j--;
else {
Debug.Assert(which[internalIdx[i]] > 0);
Debug.Assert(which[internalIdx[j]] < 0);
// swap
int tmp = internalIdx[i];
internalIdx[i] = internalIdx[j];
internalIdx[j] = tmp;
i++;
j--;
}
}
Debug.Assert(j + 1 == i);
Debug.Assert(i <= endIdx);
Debug.Assert(startIdx <= j);
tree[curTreeNodeIdx].varName = splittingVar;
tree[curTreeNodeIdx].val = threshold;
splitIdx = j;
}
private bool FindBestVariableAndThreshold(int startIdx, int endIdx, out double threshold, out string bestVar) {
Debug.Assert(startIdx < endIdx + 1); // at least 2 elements
int rows = endIdx - startIdx + 1;
Debug.Assert(rows >= 2);
double sumY = 0.0;
for (int i = startIdx; i <= endIdx; i++) {
sumY += y[internalIdx[i]];
}
double bestImprovement = 1.0 / rows * sumY * sumY;
double bestThreshold = double.PositiveInfinity;
bestVar = RegressionTreeModel.TreeNode.NO_VARIABLE;
for (int col = 0; col < effectiveVars; col++) {
// sort values for variable to prepare for threshold selection
var curVariable = allowedVariables[col];
var curVariableIdx = varName2Index[curVariable];
for (int i = startIdx; i <= endIdx; i++) {
var sortedI = sortedIdx[curVariableIdx][i];
outSortedIdx[i - startIdx] = sortedI;
outx[i - startIdx] = x[curVariableIdx][sortedI];
}
double curImprovement;
double curThreshold;
FindBestThreshold(outx, outSortedIdx, rows, y, sumY, out curThreshold, out curImprovement);
if (curImprovement > bestImprovement) {
bestImprovement = curImprovement;
bestThreshold = curThreshold;
bestVar = allowedVariables[col];
}
}
if (bestVar == RegressionTreeModel.TreeNode.NO_VARIABLE) {
threshold = bestThreshold;
return false;
} else {
UpdateVariableRelevance(bestVar, sumY, bestImprovement, rows);
threshold = bestThreshold;
return true;
}
}
// TODO: assumption is that the Average(y) = 0
private void UpdateVariableRelevance(string bestVar, double sumY, double bestImprovement, int rows) {
if (string.IsNullOrEmpty(bestVar)) return;
// update variable relevance
double err = sumY * sumY / rows;
double errAfterSplit = bestImprovement;
double delta = (errAfterSplit - err); // relative reduction in squared error
double v;
if (!sumImprovements.TryGetValue(bestVar, out v)) {
sumImprovements[bestVar] = delta;
}
sumImprovements[bestVar] = v + delta;
}
// x [0..N-1] contains rows sorted values in the range from [0..rows-1]
// sortedIdx [0..N-1] contains the idx of the values in x in the original dataset in the range from [0..rows-1]
// rows specifies the number of valid entries in x and sortedIdx
// y [0..N-1] contains the target values in original sorting order
// sumY is y.Sum()
//
// the routine returns the best threshold (x[i] + x[i+1]) / 2 for i = [0 .. rows-2] by calculating the reduction in squared error
// additionally the reduction in squared error is returned in bestImprovement
// if all elements of x are equal the routing fails to produce a threshold
private static void FindBestThreshold(double[] x, int[] sortedIdx, int rows, double[] y, double sumY, out double bestThreshold, out double bestImprovement) {
Debug.Assert(rows >= 2);
double sl = 0.0;
double sr = sumY;
double nl = 0.0;
double nr = rows;
bestImprovement = 1.0 / rows * sumY * sumY;
bestThreshold = double.NegativeInfinity;
// for all thresholds
// if we have n rows there are n-1 possible splits
for (int i = 0; i < rows - 1; i++) {
sl += y[sortedIdx[i]];
sr -= y[sortedIdx[i]];
nl++;
nr--;
Debug.Assert(nl > 0);
Debug.Assert(nr > 0);
if (x[i] < x[i + 1]) { // don't try to split when two elements are equal
// goal is to find the split with leading to minimal total variance of left and right parts
// without partitioning the variance is var(y) = E(y²) - E(y)²
// = 1/n * sum(y²) - (1/n * sum(y))²
// -------------
// if we split into right and left part the overall variance is the weigthed combination nl/n * var(y_l) + nr/n * var(y_r)
// = nl/n * (1/nl * sum(y_l²) - (1/nl * sum(y_l))²) + nr/n * (1/nr * sum(y_r²) - (1/nr * sum(y_r))²)
// = 1/n * sum(y_l²) - 1/nl * 1/n * sum(y_l)² + 1/n * sum(y_r²) - 1/nr * 1/n * sum(y_r)²
// = 1/n * (sum(y_l²) + sum(y_r²)) - 1/n * (sum(y_l)² / nl + sum(y_r)² / nr)
// = 1/n * sum(y²) - 1/n * (sum(y_l)² / nl + sum(y_r)² / nr)
// -------------
// not changed by split (and the same for total variance without partitioning)
//
// therefore we need to find the maximum value (sum(y_l)² / nl + sum(y_r)² / nr) (ignoring the factor 1/n)
// and this value must be larger than 1/n * sum(y)² to be an improvement over no split
double curQuality = sl * sl / nl + sr * sr / nr;
if (curQuality > bestImprovement) {
bestThreshold = (x[i] + x[i + 1]) / 2.0;
bestImprovement = curQuality;
}
}
}
// if all elements where the same then no split can be found
}
public IEnumerable> GetVariableRelevance() {
// values are scaled: the most important variable has relevance = 100
double scaling = 100 / sumImprovements.Max(t => t.Value);
return
sumImprovements
.Select(t => new KeyValuePair(t.Key, t.Value * scaling))
.OrderByDescending(t => t.Value);
}
}
}