Source code for equadratures.polytree

import numpy as np
from copy import deepcopy
from equadratures.parameter import Parameter
from equadratures import Weight
from equadratures.poly import Poly
from equadratures.basis import Basis
import equadratures.plot as plot
from urllib.parse import quote

[docs]class PolyTree(object): """ Definition of a polynomial tree object. Parameters ---------- splitting_criterion : str, optional The type of splitting_criterion to use in the fit function. Options include ``model_aware`` which fits polynomials for each candidate split, ``model_agnostic`` which uses a standard deviation based model-agnostic split criterion [1], and ``loss_gradient`` which uses a gradient based splitting criterion similar to that in [2]. max_depth : int, optional The maximum depth which the tree will grow to. min_samples_leaf : int, optional The minimum number of samples per leaf node. order : int, optional The order of the generated orthogonal polynomials. basis : str, optional The type of index set used for the basis. Options include: ``univariate``, ``total-order``, ``tensor-grid``, ``sparse-grid`` and ``hyperbolic-basis``. search : str, optional The method of search to be used. Options are ``grid`` or ``exhaustive``. samples : int, optional The interval between splits if ``grid`` search is chosen. verbose : bool, optional For debugging. all_data : bool, optional Store data at all nodes in :class:`~PolyTree` (instead of only leaf nodes). split_dims : list, optional List of dimensions along which to make splits. k : float, optional The smoothing parameter. Range from 0.0 to 1.0, with 0 giving no smoothing, and 1 giving maximum smoothing. distribution : str, optional The type of input parameter distributions. Either ``uniform`` or ``data``. Example ------- >>> tree = polytree.PolyTree() >>> X = np.loadtxt('inputs.txt') >>> Xtest = np.loadtxt('inputs_test.txt') >>> y = np.loadtxt('outputs.txt') >>>,y) >>> y_test = tree.predict(X_test) References ---------- 1. Wang, Y., Witten, I. H., (1997) Inducing Model Trees for Continuous Classes. In Proc. of the 9th European Conf. on Machine Learning Poster Papers. 128-137. `Paper <>`__ 2. Broelemann, K., Kasneci, G., (2019) A Gradient-Based Split Criterion for Highly Accurate and Transparent Model Trees. In Int. Joint Conf. on Artificial Intelligence (IJCAI). 2030-2037. `Paper <>`__ 3. Chan, T. F., Golub, G. H., LeVeque, R. J., (1983) Algorithms for computing the sample variance: Analysis and recommendations. The American Statistician. 37(3): 242–247. `Paper <>`__ """ def __init__(self, splitting_criterion='model_aware', max_depth=5, min_samples_leaf=None, order=1, basis='total-order', search='exhaustive', samples=50, verbose=False, poly_method="least-squares", poly_solver_args={},all_data=False,split_dims=None,k=0.05,distribution='uniform'): self.splitting_criterion = splitting_criterion self.max_depth = max_depth self.min_samples_leaf = min_samples_leaf self.order = order self.basis = basis self.tree = None = search self.samples = samples self.verbose = verbose self.cardinality = None self.poly_method = poly_method self.poly_solver_args = poly_solver_args self.actual_max_depth = 0 self.all_data = all_data self.k = k self.distribution = distribution if split_dims is not None: split_dims = [split_dims] if not isinstance(split_dims, list) else split_dims assert all(isinstance(dim, int) for dim in split_dims), "split_dims should be a list if ints" self.split_dims = split_dims assert max_depth >= 0, "max_depth must be >= 0" assert order > 0, "order must be a postive integer" assert samples > 0, "samples must be a postive integer" assert k > 0, "k must be a positive number"
[docs] def get_splits(self): """ Returns all of the data splits made. Returns ------- list A list of splits made in the format of a nested list: [[split, dimension], ...] """ def _search_tree(node, splits): if node["children"]["left"] != None: if [node["threshold"], node["j_feature"]] not in splits: splits.append([node["threshold"], node["j_feature"]]) splits = _search_tree(node["children"]["left"], splits) if node["children"]["right"] != None: if [node["threshold"], node["j_feature"]] not in splits: splits.append([node["threshold"], node["j_feature"]]) splits = _search_tree(node["children"]["right"], splits) return splits return _search_tree(self.tree, [])
def _split_data(self, j_feature, threshold, X, y): idx_left = np.where(X[:, j_feature] <= threshold)[0] idx_right = np.delete(np.arange(0, len(X)), idx_left) assert len(idx_left) + len(idx_right) == len(X) return (X[idx_left], y[idx_left]), (X[idx_right], y[idx_right])
[docs] def get_polys(self): """ Returns all of the polynomials fitted at each node in the tree. Returns ------- list A list of Poly objects. """ def _search_tree(node, polys): if node["children"]["left"] == None and node["children"]["right"] == None: polys.append(node["poly"]) if node["children"]["left"] != None: polys = _search_tree(node["children"]["left"], polys) if node["children"]["right"] != None: polys = _search_tree(node["children"]["right"], polys) return polys return _search_tree(self.tree, [])
[docs] def fit(self, X, y): """ Fits the PolyTree to the provided data. Parameters ---------- X : numpy.ndarray Training input data y : numpy.ndarray Training output data """ def _build_tree(): global index_node_global def _splitter(node): # Extract data X, y = node["data"] depth = node["depth"] N, d = X.shape # Dimensions to split along if self.split_dims is None: self.split_dims = range(d) # Find feature splits that might improve loss did_split = False if self.splitting_criterion == "model_aware": loss_best = node["loss"] elif self.splitting_criterion == "model_agnostic" or self.splitting_criterion=="loss_gradient": loss_best = np.inf else: raise Exception("invalid splitting_criterion") data_best = None polys_best = None j_feature_best = None threshold_best = None if self.verbose: polys_fit = 0 # Perform threshold split search only if node has not hit max depth if (depth >= 0) and (depth < self.max_depth): if self.splitting_criterion != "loss_gradient": for j_feature in range(d): last_threshold = np.inf if == 'exhaustive': threshold_search = X[:, j_feature] elif == 'grid': if self.samples > N: samples = N else: samples = self.samples threshold_search = np.linspace(np.min(X[:,j_feature]), np.max(X[:,j_feature]), num=samples) else: raise Exception('Incorrect search type! Must be \'exhaustive\' or \'grid\'') # Perform threshold split search on j_feature for threshold in np.unique(np.sort(threshold_search)): # Split data based on threshold (X_left, y_left), (X_right, y_right) = self._split_data(j_feature, threshold, X, y) #print(j_feature, threshold, X_left, X_right) N_left, N_right = len(X_left), len(X_right) # Do not attempt to split if split conditions not satisfied if not (N_left >= self.min_samples_leaf and N_right >= self.min_samples_leaf): continue # Compute weight loss function if self.splitting_criterion == "model_aware": loss_left, poly_left = _fit_poly(X_left, y_left) loss_right, poly_right = _fit_poly(X_right, y_right) loss_split = (N_left*loss_left + N_right*loss_right) / N if self.verbose: polys_fit += 2 elif self.splitting_criterion == "model_agnostic": loss_split = np.std(y) - (N_left*np.std(y_left) + N_right*np.std(y_right)) / N # Update best parameters if loss is lower if loss_split < loss_best: did_split = True loss_best = loss_split if self.splitting_criterion == "model_aware": polys_best = [poly_left, poly_right] data_best = [(X_left, y_left), (X_right, y_right)] j_feature_best = j_feature threshold_best = threshold # Gradient based splitting criterion from ref. [2] else: # Fit a single poly to parent node loss, poly = _fit_poly(X, y) # Now run the splitting algo using gradients from this poly did_split, j_feature_best, threshold_best = self._find_split_from_grad(poly, X, y.reshape(-1,1)) # If model_agnostic or gradient based, fit poly's to children now we have split if self.splitting_criterion != "model_aware" and did_split: (X_left, y_left), (X_right, y_right) = self._split_data(j_feature_best, threshold_best, X, y) loss_left, poly_left = _fit_poly(X_left, y_left) loss_right, poly_right = _fit_poly(X_right, y_right) N_left, N_right = len(X_left), len(X_right) loss_best = (N_left*loss_left + N_right*loss_right) / N polys_best = [poly_left, poly_right] if self.splitting_criterion == "loss_gradient": data_best = [(X_left, y_left), (X_right, y_right)] if self.verbose: polys_fit += 2 if self.verbose and did_split: print("Node (X.shape = {}) fitted with {} polynomials generated".format(X.shape, polys_fit)) elif self.verbose: print("Node (X.shape = {}) failed to fit after {} polynomials generated".format(X.shape, polys_fit)) if did_split and depth > self.actual_max_depth: self.actual_max_depth = depth # Return the best result result = {"did_split": did_split, "loss": loss_best, "polys": polys_best, "data": data_best, "j_feature": j_feature_best, "threshold": threshold_best, "N": N} return result def _fit_poly(X, y): # try: N, d = X.shape myParameters = [] for dimension in range(d): values = X[:,dimension] values_min = np.amin(values) values_max = np.amax(values) if (values_min - values_max) ** 2 < 0.01: values_min -= 0.01 values_max += 0.01 myParameters.append(Parameter(distribution='Uniform', lower=values_min, upper=values_max, order=self.order)) else: if self.distribution == 'uniform': myParameters.append(Parameter(distribution='Uniform', lower=values_min, upper=values_max, order=self.order)) elif self.distribution == 'data': input_dist = Weight(values, support=[values_min, values_max], pdf=False) myParameters.append(Parameter(distribution='data',weight_function=input_dist,order=self.order)) if self.basis == "hyperbolic-basis": myBasis = Basis(self.basis, orders=[self.order for _ in range(d)], q=0.5) else: myBasis = Basis(self.basis, orders=[self.order for _ in range(d)]) container["index_node_global"] += 1 poly = Poly(myParameters, myBasis, method=self.poly_method, sampling_args={'sample-points':X, 'sample-outputs':y}, solver_args=self.poly_solver_args) poly.set_model() mse = np.linalg.norm(y - poly.get_polyfit(X).reshape(-1)) ** 2 / N # except Exception as e: # print("Warning fitting of Poly failed:", e) # print(d, values_min, values_max) # mse, poly = np.inf, None return mse, poly def _create_node(X, y, depth, container): poly_loss, poly = _fit_poly(X, y) node = {"name": "node", "index": container["index_node_global"], "loss": poly_loss, "poly": poly, "data": (X, y), "n_samples": len(X), "j_feature": None, "threshold": None, "children": {"left": None, "right": None}, "depth": depth, "flag": False} container["index_node_global"] += 1 return node def _split_traverse_node(node, container): result = _splitter(node) if not result["did_split"]: return node["j_feature"] = result["j_feature"] node["threshold"] = result["threshold"] if not self.all_data: del node["data"] (X_left, y_left), (X_right, y_right) = result["data"] poly_left, poly_right = result["polys"] node["children"]["left"] = _create_node(X_left, y_left, node["depth"]+1, container) node["children"]["right"] = _create_node(X_right, y_right, node["depth"]+1, container) node["children"]["left"]["poly"] = poly_left node["children"]["right"]["poly"] = poly_right # Split nodes _split_traverse_node(node["children"]["left"], container) _split_traverse_node(node["children"]["right"], container) container = {"index_node_global": 0} root = _create_node(X, y, 0, container) _split_traverse_node(root, container) return root N, d = X.shape if self.basis == "hyperbolic-basis": self.cardinality = Basis(self.basis, orders=[self.order for _ in range(d)], q=0.5).get_cardinality() else: self.cardinality = Basis(self.basis, orders=[self.order for _ in range(d)]).get_cardinality() if self.min_samples_leaf == None or self.min_samples_leaf == self.cardinality: self.min_samples_leaf = int(np.ceil(self.cardinality * 1.25)) elif self.cardinality > self.min_samples_leaf: print("WARNING: Basis cardinality ({}) greater than the minimum samples per leaf ({}). This may cause reduced performance.".format(self.cardinality, self.min_samples_leaf)) self.k *= self.min_samples_leaf self.tree = _build_tree()
[docs] def prune(self, X, y, tol=0.0, percent=False): """ Prunes the tree that you have fitted. Parameters ---------- X : numpy.ndarray Training input data y : numpy.ndarray Training output data tol : float, optional Pruning tolerance (%). Prune nodes if they only improve loss by less than this tolerance. percent : bool, optional If true, tol is taken as a percentage of the parent node's error. Otherwise, tol is taken to be an absolute value. """ if percent: tol /= 100.0 def pruner(node, X_subset, y_subset): if X_subset.shape[0] < 1: node["test_loss"] = 0 node["n_samples"] = 0 return node node["test_loss"] = np.linalg.norm(y_subset - node["poly"].get_polyfit(X_subset).reshape(-1)) ** 2 / X_subset.shape[0] is_left = node["children"]["left"] != None is_right = node["children"]["right"] != None if is_left and is_right: (X_left, y_left), (X_right, y_right) = self._split_data(node["j_feature"], node["threshold"], X_subset, y_subset) node["children"]["left"] = pruner(node["children"]["left"], X_left, y_left) node["children"]["right"] = pruner(node["children"]["right"], X_right, y_right) lower_loss = ( node["children"]["left"]["test_loss"] * node["children"]["left"]["n_samples"] + node["children"]["right"]["test_loss"] * node["children"]["right"]["n_samples"] ) / ( node["children"]["left"]["n_samples"] + node["children"]["right"]["n_samples"] ) node["lower_loss"] = lower_loss if percent: loss_eps = tol * node["test_loss"] else: loss_eps = tol print(lower_loss + loss_eps, node["test_loss"]) if lower_loss + loss_eps > node["test_loss"]: if self.verbose: print("prune",lower_loss, node["test_loss"], node["children"]["left"]["test_loss"], node["children"]["left"]["n_samples"], node["children"]["right"]["test_loss"], node["children"]["right"]["n_samples"]) node["children"]["left"] = None node["children"]["right"] = None return node assert self.tree is not None, "Run fit() before prune()" (X_left, y_left), (X_right, y_right) = self._split_data(self.tree["j_feature"], self.tree["threshold"], X, y) self.tree["children"]["left"] = pruner(self.tree["children"]["left"], X_left, y_left) self.tree["children"]["right"] = pruner(self.tree["children"]["right"], X_right, y_right)
[docs] def predict(self, X): """ Evaluates the the polynomial tree approximation of the data. Parameters ---------- X : numpy.ndarray An ndarray with shape (number_of_observations, dimensions) at which the tree fit must be evaluated at. Returns ------- numpy.ndarray Array with shape (1, number_of_observations) corresponding to the polynomial approximations of the tree. """ def _predict(node, indexes): y_pred[indexes, node["depth"], 0] = node["poly"].get_polyfit(X[indexes]).reshape(-1) y_pred[indexes, node["depth"], 1] = np.full(fill_value=node["n_samples"], shape=len(indexes)) no_children = node["children"]["left"] is None and \ node["children"]["right"] is None if no_children: return idx_left = np.where(X[indexes, node["j_feature"]] <= node["threshold"])[0] idx_right = np.where(X[indexes, node["j_feature"]] > node["threshold"])[0] _predict(node["children"]["left"], indexes[idx_left]) _predict(node["children"]["right"], indexes[idx_right]) assert self.tree is not None y_pred = np.empty(shape=(X.shape[0], self.actual_max_depth + 2, 2)) * np.nan _predict(self.tree, np.arange(0, X.shape[0])) smoothed_y_pred = np.zeros(shape=(X.shape[0])) for y in range(0,X.shape[0]): i = self.actual_max_depth + 1 while np.isnan(y_pred[y][i][0]) and i > 0: i-=1 smoothed_y = y_pred[y][i][0] #print(y_pred[i]) while i > 0: n_i = y_pred[y][i][1] if n_i == 0: break #print(smoothed_y) smoothed_y = (smoothed_y * n_i + y_pred[y][i][0] * self.k) / (self.k + n_i) i-=1 #print("\n") smoothed_y_pred[y] = smoothed_y return smoothed_y_pred
[docs] def apply(self,X): """ Returns the leaf node index for each observation in the data. Parameters ---------- X : numpy.ndarray Array with shape (number_of_observations, dimensions) at which the tree fit must be evaluated at. Returns ------- numpy.ndarray A numpy.ndarray of shape (number_of_observations,1) corresponding to the node indices for each observation in X. """ def _apply(node, indexes): no_children = node["children"]["left"] is None and \ node["children"]["right"] is None if no_children: inode[indexes] = node["index"] return idx_left = np.where(X[indexes, node["j_feature"]] <= node["threshold"])[0] idx_right = np.where(X[indexes, node["j_feature"]] > node["threshold"])[0] _apply(node["children"]["left"], indexes[idx_left]) _apply(node["children"]["right"], indexes[idx_right]) if X.ndim == 1: X = X.reshape(1,-1) inode = np.zeros(shape=X.shape[0],dtype=int) _apply(self.tree, np.arange(0, X.shape[0])) return inode
[docs] def get_leaves(self): """ Returns the node indices for all leaf nodes. Returns ------- list Contains the node indices of all leaf nodes. """ def _recurse(node,leaf_list): no_children = node["children"]["left"] is None and \ node["children"]["right"] is None if no_children: leaf_list.append(node["index"]) return _recurse(node["children"]["left"],leaf_list) _recurse(node["children"]["right"],leaf_list) leaf_list = [] _recurse(self.tree,leaf_list) return leaf_list
[docs] def get_mean_and_variance(self): """ Computes the mean and variance of the polynomial tree model. Returns ------- tuple Tuple (mean,variance) containing two floats; the approximated mean and variance from the fitted PolyTree. """ # Get volume of polytree domain root_poly = self.tree["poly"] root_vol = self._calc_domain_vol(root_poly) # Get leaf nodes leaves = self.get_leaves() # Summation over all leaf nodes in the tree mean = 0. var = 0. for leaf in leaves: leaf_poly = self.get_node(leaf)["poly"] leaf_vol = self._calc_domain_vol(leaf_poly) coeffs = leaf_poly.coefficients # Compute mean mean += (leaf_vol/root_vol) * float(coeffs[0]) # Compute variance tmp = 0. for i in range(0, len(coeffs)): tmp += float(coeffs[i]**2) var += (leaf_vol/root_vol) * tmp var -= mean**2 return mean, var
[docs] def get_graphviz(self, X=None, feature_names=None, file_name=None): """ Generates a graphviz visualisation of the PolyTree. Parameters ---------- X : numpy.ndarray, optional An ndarray with shape (dimensions) containing an input vector for a given sample, to highlight in the tree. feature_names : list, optional A list of the names of the features used in the training data. filename : str, optional Filename to write graphviz data to. If ``None`` (default) then rendered in-place, if ``'source'``, the raw graphviz string is returned. """ from graphviz import Digraph g = Digraph('g', node_attr={'shape': 'record', 'height': '.1'}) if feature_names is None: dim = self.tree["poly"].dimensions feature_names = ['x_%d'%i for i in range(dim)] def _build_graphviz_recurse(node, parent_node_index=0, parent_depth=0, edge_label="",labelangle=0): # Empty node if node is None: return # Create node node_index = node["index"] if node["children"]["left"] is None and node["children"]["right"] is None: threshold_str = "" leaf = True else: threshold_str = "{} <= {:.3f}\\n".format(feature_names[node['j_feature']], node["threshold"]) leaf = False if "lower_loss" in node: label_str = "node {} \\n {} n_samples = {}\\n loss = {:.6f}\\n lower_loss = {}".format(node_index,threshold_str, node["n_samples"], node["test_loss"], node["lower_loss"]) elif "test_loss" in node: label_str = "node {} \\n {} n_samples = {}\\n loss = {:.6f}".format(node_index,threshold_str, node["n_samples"], node["test_loss"]) else: label_str = "node {} \\n {} n_samples = {}\\n loss = {:.6f}".format(node_index,threshold_str, node["n_samples"], node["loss"]) # Create node if leaf: nodeshape = "rectangle" style = ["rounded"] fillcolor = "#E4fEE4" else: nodeshape = "rectangle" style = ["filled"] fillcolor = "#EBFAFF" if node["flag"]: style.append('bold') bordercolor = "black" fontcolor = "black" g.attr('node', label=label_str, shape=nodeshape) g.node('{}'.format(node_index), color=bordercolor, style=', '.join(style), fillcolor=fillcolor, fontcolor=fontcolor) # Create edge if parent_depth > 0: if node["flag"]: edgecolor = 'orange' style = 'bold' else: edgecolor = 'black' style = 'solid' if parent_depth > 1: edge_label = '' # Only label True/False for root node g.edge('{}'.format(parent_node_index), '{}'.format(node_index), headlabel=edge_label, color=edgecolor,style=style,labeldistance="2.5",labelangle=labelangle) # Traverse child or append leaf value _build_graphviz_recurse(node["children"]["left"], parent_node_index=node_index, parent_depth=parent_depth + 1, edge_label="True",labelangle="45") _build_graphviz_recurse(node["children"]["right"], parent_node_index=node_index, parent_depth=parent_depth + 1, edge_label="False",labelangle="-45") def _flag_tree_walk(node,X): node["flag"] = True if node["children"]["left"] is None and \ node["children"]["right"] is None: return if X[node["j_feature"]] <= node["threshold"]: return _flag_tree_walk(node["children"]["left"],X) if X[node["j_feature"]] > node["threshold"]: return _flag_tree_walk(node["children"]["right"],X) # Flag the node path to highlight later if X is not None: _flag_tree_walk(self.tree,X) # Build graph _build_graphviz_recurse(self.tree, parent_node_index=0, parent_depth=0, edge_label="") if file_name == 'source': return g.source elif file_name is None: try: g.render(view=True) except: file_name = '' print("GraphViz source file written to " + file_name + " and can be viewed using an online renderer. Alternatively you can install graphviz on your system to render locally") if file_name is not None: # not elif here as file_name might be updated in try-except above with open(file_name, "w") as file: file.write(str(g.source))
[docs] def get_node(self, inode): """ Returns the node corresponding to a given node number. Parameters ---------- inode : int The node number. Returns ------- dict Dictionary containing the data for the requested node. """ # Find node with given index inode. Traverse all children until correct node found. def _get_node_from_n(node): if node is not None: # Need to check if node is None here as below _get_node_from_n() calls on children will result in None if leaf node if node["index"] == inode: return node else: result = _get_node_from_n(node["children"]["right"]) if result is None: result = _get_node_from_n(node["children"]["left"]) return result else: return None return _get_node_from_n(self.tree)
[docs] def get_paths(self,X=None): """ Returns the tree paths for the leaf nodes in the tree. Parameters ---------- X : numpy.ndarray, optional Array with shape (number_of_observations, dimensions) to apply the tree to. If given, paths will only be returned for leaves which contain observations. Returns ------- dict Dictionary containing a dict for each leaf node. Indexed by the node indices for the leaf nodes. """ def _find_path(node, path, i): """ Private recursive function to find path through a tree for a given leaf node. """ node_index = node["index"] info = {'node':node_index,'j':node["j_feature"],'thresh':node["threshold"]} path.append(info) if node_index == i: return True left = False right = False if node["children"]["left"] is not None: left = _find_path(node["children"]["left"], path, i) if node["children"]["right"] is not None: right = _find_path(node["children"]["right"], path, i) if left or right : return True path.remove(info) return False # Get leaf node id's if X is None: leave_id = self.get_leaves() else: # Get leaf nodes leave_id = self.apply(X) # Loop through leaves and find path for each. paths ={} for leaf in np.unique(leave_id): path_leaf = [] _find_path(self.tree, path_leaf, leaf) # Set split info to None for leaf node path_leaf[-1]["j"] = None path_leaf[-1]["thresh"] = None # Save in dict paths[leaf] = path_leaf return paths
[docs] def plot_decision_surface(self,ij,ax=None,X=None,y=None,max_depth=None,label=True, color='data',colorbar=True,show=True,kwargs={}): """ Plots the decision boundaries of the PolyTree over a 2D surface. See :meth:`~equadratures.plot.plot_decision_surface` for full description. """ return plot.plot_decision_surface(self,ij,ax,X,y,max_depth,label,color,colorbar,show,kwargs)
def _find_split_from_grad(self,model, X, y): """ Private method to find the optimal split point for a tree node based on the training data in that node. Parameters ---------- model : Poly An instance of the Poly class, corresponding to the Poly belonging to the tree node. X : numpy.ndarray An ndarray with shape (number_of_observations, dimensions) containing the input data belonging to the tree node. y : numpy.ndarray An ndarray with shape (number_of_observations, 1) containing the response data belonging to the tree node. Returns ------- tuple Tuple (did_split, split_dim, split_val), where: did_split (bool): True if a split was found, otherwise False. split_dim (int): The dimension in X within which the best split was found. split_val (float): The location of the best split. """ renorm = True N,D = np.shape(X) # Gradient of loss wrt model coefficients P = model.get_poly(X).T r = y-model.get_polyfit(X) g = r*P # Sum of gradients gsum = g.sum(axis=0) # Loop through all dimensions in X split_dim = None split_val = None gain_max = -np.inf for d in self.split_dims: # Sort along feature i sort = np.argsort(X[:,d]) Xd = X[sort,d] # Find unique values along one column. #TODO - grid search option _,splits = np.unique(Xd,return_index=True) splits = splits[1:] # Number of samples on left and right split N_l = splits N_r = N - N_l # Only take splits where both children have more than `min_samples_leaf` samples idx = np.minimum(N_l, N_r) >= self.min_samples_leaf splits = splits[idx] N_l = N_l[idx].reshape(-1,1) N_r = N_r[idx].reshape(-1,1) # If we've run out of candidate spilts, skip if len(splits) <= 1: continue # Sums of gradients for left and right gsum_left = g[sort,:].cumsum(axis=0) gsum_left = gsum_left[splits-1,:] gsum_right = gsum - gsum_left # Renorm. gradients to zero mean and unit std if renorm: mu_l, mu_r, sigma_l, sigma_r = self._get_mean_and_sigma(P[:,1:],splits,N_l,N_r,sort) gsum_left = self._renormalise( gsum_left, 1/sigma_l, -mu_l/sigma_l) gsum_right = self._renormalise(gsum_right, 1/sigma_r, -mu_r/sigma_r) # Compute the Gain (see Eq. (6) in [1]) gain = (gsum_left**2).sum(axis=1)/N_l.reshape(-1) + (gsum_right**2).sum(axis=1)/N_r.reshape(-1) # Find best gain and compare with previous best best_idx = np.argmax(gain) gain = gain[best_idx] if gain > gain_max: gain_max = gain best_split_dim = d best_split_val = 0.5*(Xd[splits[best_idx] - 1] + Xd[splits[best_idx]]) # If gain_max stilll == -np.inf, we must have passed through all features w/o finding a split # so return False. Otherwise return True and the spilt details. if gain_max == -np.inf: return False, None, None else: return True, best_split_dim, best_split_val @staticmethod def _get_mean_and_sigma(X,splits,N_l,N_r,sort): """ Computes mean and standard deviation of the data in array X, when it is split in two by the threshold values in the splits array. The data is offset by its mean to avoid catastrophic cancellation when computing the variance (see ref. [3]). Parameters ---------- X : numpy.ndarray Arrray with dimensions (N,ndim) containing the orthogonal polynomials P. splits : numpy.ndarray Array of split locations. N_l : numpy.ndarray Array containing info on number of samples to left of splits. N_r : numpy.ndarray Array containing info on number of samples to right of splits. sort : numpy.ndarray Index array to reorder X. """ # Min value of sigma (for stability later) epsilon = 0.001 # Reorder, and shift X by mean mu = np.reshape(np.mean(X, axis=0), (1, -1)) Xshift = X[sort] - mu # Cumulative sums (and sums of squares) for left and right splits Xsum_l = Xshift.cumsum(axis=0) Xsum_r = Xsum_l[-1:,:] - Xsum_l X2sum_l = (Xshift**2).cumsum(axis=0) X2sum_r = X2sum_l[-1:,:] - X2sum_l # Compute mean of left and right side for all splits mu_l = Xsum_l[splits-1,:] / N_l mu_r = Xsum_r[splits-1,:] / N_r # Compute standard deviation of left and right side for all splits sigma_l = np.sqrt(np.maximum(X2sum_l[splits-1,:]/(N_l-1)-mu_l**2, epsilon**2)) sigma_r = np.sqrt(np.maximum(X2sum_r[splits-1,:]/(N_r-1)-mu_r**2, epsilon**2)) # Correct for previous shift mu_l = mu_l + mu mu_r = mu_r + mu return mu_l, mu_r, sigma_l, sigma_r @staticmethod def _renormalise(gradients, a, c): """ Renormalises gradients according to according to eq. (14) of [1]. Parameters ---------- gradients : numpy.ndarray Array with shape (n_samples, n_params), containing the gradients. a : numpy.ndarray Array with shape (n_samples, n_params-1) containing the normalisation factors. c: numpy.ndarray Array with shape (n_samples, n_params-1) containing the normalisation offsets. Returns ------- gradients : numpy.ndarray Array with shape (n_samples, n_params) containing the renormalised gradients. """ c = c*gradients[:,0].reshape(-1,1) gradients[:,1:] = gradients[:,1:] * a + c return gradients @staticmethod def _calc_domain_vol(Polynomial): params = Polynomial.parameters vol = 1. for param in params: vol *= param.upper - param.lower return vol