| #------------------------------------------------------------- |
| # |
| # Licensed to the Apache Software Foundation (ASF) under one |
| # or more contributor license agreements. See the NOTICE file |
| # distributed with this work for additional information |
| # regarding copyright ownership. The ASF licenses this file |
| # to you under the Apache License, Version 2.0 (the |
| # "License"); you may not use this file except in compliance |
| # with the License. You may obtain a copy of the License at |
| # |
| # http://www.apache.org/licenses/LICENSE-2.0 |
| # |
| # Unless required by applicable law or agreed to in writing, |
| # software distributed under the License is distributed on an |
| # "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY |
| # KIND, either express or implied. See the License for the |
| # specific language governing permissions and limitations |
| # under the License. |
| # |
| #------------------------------------------------------------- |
| |
| # This script implements decision trees for recoded and binned categorical and |
| # numerical input features. We train a single CART (classification and |
| # regression tree) decision trees depending on the provided labels y, either |
| # classification (majority vote per leaf) or regression (average per leaf). |
| # |
| # .. code-block:: |
| # |
| # For example, give a feature matrix with features [a,b,c,d] |
| # and the following trees, M would look as follows: |
| # |
| # (L1) |d<5| |
| # / \ |
| # (L2) P1:2 |a<7| |
| # / \ |
| # (L3) P2:2 P3:1 |
| # |
| # --> M := |
| # [[4, 5, 0, 2, 1, 7, 0, 0, 0, 0, 0, 2, 0, 1]] |
| # |(L1)| | (L2) | | (L3) | |
| # |
| # |
| # |
| # INPUT: |
| # ------------------------------------------------------------------------------ |
| # X Feature matrix in recoded/binned representation |
| # y Label matrix in recoded/binned representation |
| # ctypes Row-Vector of column types [1 scale/ordinal, 2 categorical] |
| # of shape 1-by-(ncol(X)+1), where the last entry is the y type |
| # max_depth Maximum depth of the learned tree (stopping criterion) |
| # min_leaf Minimum number of samples in leaf nodes (stopping criterion), |
| # odd number recommended to avoid 50/50 leaf label decisions |
| # min_split Minimum number of samples in leaf for attempting a split |
| # max_features Parameter controlling the number of features used as split |
| # candidates at tree nodes: m = ceil(num_features^max_features) |
| # max_values Parameter controlling the number of values per feature used |
| # as split candidates: nb = ceil(num_values^max_values) |
| # max_dataratio Parameter in [0,1] controlling when to materialize data |
| # subsets of X and y on node splits. When set to 0, we always |
| # scan the original X and y, which has the benefit of avoiding |
| # the allocation and maintenance of data for all active nodes. |
| # When set to 0.01 we rematerialize whenever the sub-tree data |
| # would be less than 1% of last the parent materialize data size. |
| # impurity Impurity measure: entropy, gini (default), rss (regression) |
| # seed Fixed seed for randomization of samples and split candidates |
| # verbose Flag indicating verbose debug output |
| # ------------------------------------------------------------------------------ |
| # |
| # OUTPUT: |
| # ------------------------------------------------------------------------------ |
| # M Matrix M containing the learned trees, in linearized form |
| # ------------------------------------------------------------------------------ |
| |
| m_decisionTree = function(Matrix[Double] X, Matrix[Double] y, Matrix[Double] ctypes, |
| Int max_depth = 10, Int min_leaf = 20, Int min_split = 50, |
| Double max_features = 0.5, Double max_values = 1.0, Double max_dataratio = 0.25, |
| String impurity = "gini", Int seed = -1, Boolean verbose = FALSE) |
| return(Matrix[Double] M) |
| { |
| t1 = time(); |
| |
| # validation checks |
| if( max_depth > 32 ) |
| stop("decisionTree: invalid max_depth > 32: "+max_depth); |
| if( sum(X<=0) != 0 ) |
| stop("decisionTree: feature matrix X is not properly recoded/binned (values <= 0): "+sum(X<=0)); |
| if( sum(abs(X-round(X))>1e-14) != 0 ) |
| stop("decisionTree: feature matrix X is not properly recoded/binned (non-integer): "+sum(abs(X-round(X))>1e-14)); |
| if( sum(y<=0) != 0 ) |
| stop("decisionTree: label vector y is not properly recoded/binned: "+sum(y<=0)); |
| |
| # initialize input data and basic statistics |
| # (we keep y2 and the indicates I in transposed form for sparsity exploitation) |
| m = nrow(X); n = ncol(X); |
| classify = (as.scalar(ctypes[1,n+1]) == 2); |
| |
| fdom = max(colMaxs(X),2); # num distinct per feature |
| foffb = t(cumsum(t(fdom))) - fdom; # feature begin |
| foffe = t(cumsum(t(fdom))) # feature end |
| rix = matrix(seq(1,m)%*%matrix(1,1,n), m*n, 1) |
| cix = matrix(X + foffb, m*n, 1); |
| X2 = table(rix, cix, 1, m, as.scalar(foffe[,n]), FALSE); #one-hot encoded |
| y2 = table(seq(1,m), y); |
| cnt = colSums(X2); |
| I = matrix(1, rows=1, cols=nrow(X)); |
| |
| if( verbose ) { |
| print("decisionTree: initialize with max_depth=" + max_depth + ", max_features=" |
| + max_features +", max_dataratio" + max_dataratio + ", impurity=" |
| + impurity + ", seed=" + seed + "."); |
| print("decisionTree: basic statistics:"); |
| print("-- impurity: " + as.scalar(computeImpurity(y2, I, impurity))); |
| print("-- minFeatureCount: " + min(cnt)); |
| print("-- maxFeatureCount: " + max(cnt)); |
| } |
| |
| # queue-based node splitting |
| M = matrix(0, rows=1, cols=2*(2^max_depth-1)) |
| queue = list(list(1,I,X2,y2)); # node IDs / data indicators |
| maxPath = 1; |
| while( length(queue) > 0 ) { |
| # pop next node from queue for splitting |
| [queue, node0] = remove(queue, 1); |
| node = as.list(node0); |
| nID = as.scalar(node[1]); |
| nI = as.matrix(node[2]); |
| X2 = as.matrix(node[3]); |
| y2 = as.matrix(node[4]); |
| if(verbose) |
| print("decisionTree: attempting split of node "+nID+" ("+sum(nI)+" rows)"); |
| |
| # optional rematerialization of data per node |
| if( sum(nI) < max_dataratio*ncol(nI) ) { |
| if(verbose) |
| print("-- compacting data: "+ncol(nI)+" --> "+sum(nI)); |
| X2 = removeEmpty(target=X2, margin="rows", select=t(nI)); |
| y2 = removeEmpty(target=y2, margin="rows", select=t(nI)); |
| nI = matrix(1, rows=1, cols=nrow(X2)); |
| } |
| |
| # find best split attribute |
| nSeed = ifelse(seed==-1, seed, seed*nID); |
| [f, v, IDleft, Ileft, IDright, Iright] = findBestSplit( |
| X2, y2, foffb, foffe, nID, nI, min_leaf, max_features, max_values, impurity, nSeed); |
| validSplit = sum(Ileft) >= min_leaf & sum(Iright) >= min_leaf; |
| if(verbose) |
| print("-- best split: f"+f+" <= "+v+" --> valid="+validSplit); |
| if( validSplit ) |
| M[, 2*nID-1:2*nID] = t(as.matrix(list(f,v))); |
| else |
| M[, 2*nID] = computeLeafLabel(y2, nI, classify, verbose); |
| maxPath = max(maxPath, floor(log(nID,2)+1)); |
| |
| # split data, finalize or recurse |
| if( validSplit ) { |
| if( sum(Ileft) >= min_split & floor(log(IDleft,2))+2 < max_depth ) |
| queue = append(queue, list(IDleft,Ileft,X2,y2)); |
| else |
| M[,2*IDleft] = computeLeafLabel(y2, Ileft, classify, verbose) |
| if( sum(Iright) >= min_split & floor(log(IDright,2))+2 < max_depth ) |
| queue = append(queue, list(IDright,Iright,X2,y2)); |
| else |
| M[,2*IDright] = computeLeafLabel(y2, Iright, classify, verbose) |
| maxPath = max(maxPath, floor(log(IDleft,2)+1)); |
| } |
| } |
| |
| # summary and encoding |
| M = M[1, 1:2*(2^maxPath-1)]; |
| |
| if(verbose) { |
| print("decisionTree: final constructed tree (linearized):"); |
| print("--" + toString(M)); |
| } |
| } |
| |
| findBestSplit = function(Matrix[Double] X2, Matrix[Double] y2, Matrix[Double] foffb, Matrix[Double] foffe, |
| Int ID, Matrix[Double] I, Int min_leaf, Double max_features, Double max_values, String impurity, Int seed) |
| return(Int f, Int v, Int IDleft, Matrix[Double] Ileft, Int IDright, Matrix[Double] Iright) |
| { |
| # sample features iff max_features < 1 |
| n = ncol(foffb); |
| numI = sum(I); |
| feat = seq(1,n); |
| if( max_features < 1.0 ) { |
| rI = rand(rows=n, cols=1, seed=seed) <= (n^max_features/n); |
| feat = removeEmpty(target=feat, margin="rows", select=rI); |
| if( sum(feat) == 0 ) #sample at least one |
| feat[1,1] = round(rand(rows=1, cols=1, min=1, max=n)); |
| } |
| |
| # evaluate features and feature splits |
| # (both categorical and numerical are treated similarly by |
| # finding a cutoff point in the recoded/binned representation) |
| R = matrix(0, rows=3, cols=nrow(feat)); |
| parfor( i in 1:nrow(feat) ) { |
| f = as.scalar(feat[i]); # feature |
| beg = as.scalar(foffb[1,f])+1; # feature start in X2 |
| end = as.scalar(foffe[1,f]); # feature end in X2 |
| belen = end-beg; #numFeat - 1 |
| while(FALSE){} # make beg/end known |
| |
| # construct 0/1 predicate vectors with <= semantics |
| # find rows that match at least one value and appear in I |
| # (vectorized evaluation, each column in P is a split candidate) |
| fP = upper.tri(target=matrix(1,belen,belen), diag=TRUE); |
| vI = seq(1,belen); |
| if( max_values < 1.0 & ncol(fP)>10 ) { |
| rI2 = rand(rows=ncol(fP),cols=1,seed=seed) <= (ncol(fP)^max_values/ncol(fP)); |
| fP = removeEmpty(target=fP, margin="cols", select=t(rI2)); |
| vI = removeEmpty(target=vI, margin="rows", select=rI2); |
| } |
| P = matrix(0, ncol(X2), ncol(fP)); |
| P[beg:end-1,1:ncol(fP)] = fP; |
| Ileft = (t(X2 %*% P) * I) != 0; |
| Iright = (Ileft==0) * I; |
| |
| # compute information gain for all split candidates |
| ig = as.scalar(computeImpurity(y2, I, impurity)) |
| - rowSums(Ileft)/numI * computeImpurity(y2, Ileft, impurity) |
| - rowSums(Iright)/numI * computeImpurity(y2, Iright, impurity); |
| ig = replace(target=ig, pattern=NaN, replacement=0); |
| |
| # track best split value and index, incl validity |
| valid = (rowSums(Ileft)>=min_leaf) & (rowSums(Iright)>=min_leaf); |
| bestig = max(valid*ig); |
| bestv = ifelse(bestig>0, nrow(valid)-as.scalar(rowIndexMax(t(rev(valid*ig))))+beg, -1); |
| if( bestv >= 0 ) |
| bestv = as.scalar(vI[bestv-beg+1,1])+beg-1; |
| R[,i] = as.matrix(list(f, bestig, bestv)); |
| } |
| ix = as.scalar(rowIndexMax(R[2,])); |
| |
| # extract indicators and IDs |
| IDleft = 2 * ID; |
| IDright= 2 * ID + 1; |
| f = as.integer(as.scalar(feat[ix,1])); |
| beg = as.scalar(foffb[1,f]); |
| v = as.integer(as.scalar(R[3,ix])-beg); |
| if( max(R[2,]) > 0 ) { |
| p = table(seq(beg+1, beg+v), 1, ncol(X2), 1); |
| Ileft = (t(X2 %*% p) * I) != 0; |
| Iright = I * (Ileft==0); |
| } |
| else { # no information gain |
| Ileft = as.matrix(0); |
| Iright = as.matrix(0); |
| } |
| } |
| |
| computeImpurity = function(Matrix[Double] y2, Matrix[Double] I, String impurity) |
| return(Matrix[Double] score) |
| { |
| f = (I %*% y2) / rowSums(I); # rel. freq. per category/bin |
| score = matrix(0, nrow(I), 1); |
| if( impurity == "gini" ) |
| score = 1 - rowSums(f^2); # sum(f*(1-f)); |
| else if( impurity == "entropy" ) |
| score = rowSums(-f * log(f)); |
| else if( impurity == "rss" ) { # residual sum of squares |
| yhat = f %*% seq(1,ncol(f)); # yhat |
| res = outer(yhat, t(rowIndexMax(y2)), "-"); # yhat-y |
| score = rowSums((I * res)^2); # sum((yhat-y)^2) |
| } |
| else |
| stop("decisionTree: unsupported impurity measure: "+impurity); |
| } |
| |
| computeLeafLabel = function(Matrix[Double] y2, Matrix[Double] I, Boolean classify, Boolean verbose) |
| return(Double label) |
| { |
| f = (I %*% y2) / sum(I); |
| label = as.scalar(ifelse(classify, |
| rowIndexMax(f), f %*% seq(1,ncol(f)))); |
| if(verbose) |
| print("-- leaf node label: " + label +" ("+sum(I)*max(f)+"/"+sum(I)+")"); |
| } |
