Repository: scikit-learn/scikit-learn
Pull request: [MRG+1] change threshold computation in trees
#### Reference Issues/PRs
Fixes #10267

#### What does this implement/fix? Explain your changes.
Change threshold computation in trees. The modification should not change anything in most cases, as the computation should be equal to the previous. It should avoid cases where the previous computation could give infinite value when value are near infinite.
This created empty leaf which caused the issue #10267.

Problem statement:
GradientBoostingRegressor with huber loss sometimes fails with an `IndexError: cannot do a non-empty take from an empty axes.`
If I use the first 63726 lines of my dataset for training everything works, but if I add one line more to the training set, then I get the error. My dataset contains no NaNs. The 63727th line has no obvious differences from the others, furthermore if I use lines 63000:64000 for training then I don't get the error, which suggests that the content of line 63727 isn't directly the problem.

I have tried and failed to make a small reproducible test, so I hope someone can make sense of what is happening here and why.

#### Versions

Windows-2008ServerR2-6.1.7601-SP1
Python 3.6.3 |Intel Corporation| (default, Oct 17 2017, 23:26:12) [MSC v.1900 64 bit (AMD64)]
NumPy 1.13.3
SciPy 0.19.1
Scikit-Learn 0.19.0

#### Error message which I got when using RandomSearchCV
```
Sub-process traceback:
---------------------------------------------------------------------------
IndexError                                         Thu Dec  7 19:08:00 2017
PID: 11688        Python 3.6.3: C:\ProgramData\Anaconda3\envs\i3\python.exe
...........................................................................
C:\ProgramData\Anaconda3\envs\i3\lib\site-packages\sklearn\externals\joblib\parallel.py in __call__(self=<sklearn.externals.joblib.parallel.BatchedCalls object>)
    126     def __init__(self, iterator_slice):
    127         self.items = list(iterator_slice)
    128         self._size = len(self.items)
    129
    130     def __call__(self):
--> 131         return [func(*args, **kwargs) for func, args, kwargs in self.items]
        self.items = [(<function _fit_and_score>, (GradientBoostingRegressor(alpha=0.9, criterion='...le=1.0, verbose=0, warm_start=False), memmap([[ -1.00000000e+00,  -2.00000003e-01,  -1 ....14299998e-02,   1.42734203e+01]], dtype=float32), memmap([ 27.,  35., -78., ..., -19.,  -9.,  -4.], type=float32), {'score': <function MaxWinRate>}, array([    0,     1,     2, ..., 12997, 12998, 12999]), memmap([ 13000,  13001,  13002, ..., 618893, 618894, 618895]), 10, {'learning_rate': 0.28522087352060566, 'max_depth': 16, 'max_features': 0.32058686573551248, 'min_samples_leaf': 1, 'n_estimators': 327}), {'error_score': 'raise', 'fit_params': {}, 'return_n_test_samples': True, 'return_parameters': False, 'return_times': True, 'return_train_score': True})]
    132
    133     def __len__(self):
    134         return self._size
    135

...........................................................................
C:\ProgramData\Anaconda3\envs\i3\lib\site-packages\sklearn\externals\joblib\parallel.py in <listcomp>(.0=<list_iterator object>)
    126     def __init__(self, iterator_slice):
    127         self.items = list(iterator_slice)
    128         self._size = len(self.items)
    129
    130     def __call__(self):
--> 131         return [func(*args, **kwargs) for func, args, kwargs in self.items]
        func = <function _fit_and_score>
        args = (GradientBoostingRegressor(alpha=0.9, criterion='...le=1.0, verbose=0, warm_start=False), memmap([[ -1.00000000e+00,  -2.00000003e-01,  -1 ....14299998e-02,   1.42734203e+01]], dtype=float32), memmap([ 27.,  35., -78., ..., -19.,  -9.,  -4.], dtype=float32), {'score': <function MaxWinRate>}, array([    0,     1,     2, ..., 12997, 12998, 12999]), memmap([ 13000,  13001,  13002, ..., 618893, 618894, 618895]), 10, {'learning_rate': 0.28522087352060566, 'max_depth': 16, 'max_features': 0.32058686573551248, 'min_samples_leaf': 1, 'n_estimators': 327})
        kwargs = {'error_score': 'raise', 'fit_params': {}, 'return_n_test_samples': True, 'return_parameters': False, 'return_times': True, 'return_train_score': True}
    132
    133     def __len__(self):
    134         return self._size
    135

...........................................................................
C:\ProgramData\Anaconda3\envs\i3\lib\site-packages\sklearn\model_selection\_validation.py in _fit_and_score(estimator=GradientBoostingRegressor(alpha=0.9, criterion='...le=1.0, verbose=0, warm_start=False), X=memmap([[ -1.00000000e+00,  -2.00000003e-01, -1....14299998e-02,   1.42734203e+01]], dtype=float32), y=memmap([ 27.,  35., -78., ..., -19.,  -9.,  -4.], dtype=float32), scorer={'score': <function MaxWinRate>}, train=array([    0,     1,     2, ..., 12997, 12998, 12999]), test=memmap([ 13000,  13001,  13002, ..., 618893, 618894, 618895]), verbose=10, parameters={'learning_rate': 0.28522087352060566, 'max_depth': 16, 'max_features': 0.32058686573551248, 'min_samples_leaf': 1, 'n_estimators': 327}, fit_params={}, return_train_score=True, return_parameters=False, return_n_test_samples=True, return_times=True, error_score='raise')
    432
    433     try:
    434         if y_train is None:
    435             estimator.fit(X_train, **fit_params)
    436         else:
--> 437             estimator.fit(X_train, y_train, **fit_params)
        estimator.fit = <bound method BaseGradientBoosting.fit of Gradie...e=1.0, verbose=0, warm_start=False)>
        X_train = memmap([[ -1.00000000e+00,  -2.00000003e-01,  -1....04051296e+09,   1.89901295e+01]], dtype=float32)
        y_train = memmap([ 27.,  35., -78., ...,   9.,  21.,  20.], dtype=float32)
        fit_params = {}
    438
    439     except Exception as e:
    440         # Note fit time as time until error
    441         fit_time = time.time() - start_time

...........................................................................
C:\ProgramData\Anaconda3\envs\i3\lib\site-packages\sklearn\ensemble\gradient_boosting.py in fit(self=GradientBoostingRegressor(alpha=0.9, criterion='...le=1.0, verbose=0, warm_start=False), X=array([[ -1.00000000e+00,  -2.00000003e-01,  -1.....04051296e+09,   1.89901295e+01]], dtype=float32), y=memmap([ 27.,  35., -78., ...,   9.,  21.,  20.], dtype=float32), sample_weight=array([ 1.,  1.,  1., ...,  1.,  1.,  1.], dtype=float32), monitor=None)
   1029                 X_idx_sorted = np.asfortranarray(np.argsort(X, axis=0),
   1030                                                  dtype=np.int32)
   1031
   1032         # fit the boosting stages
   1033         n_stages = self._fit_stages(X, y, y_pred, sample_weight, random_state,
-> 1034                                     begin_at_stage, monitor, X_idx_sorted)
        begin_at_stage = 0
        monitor = None
        X_idx_sorted = array([[    0,  8619,     0, ...,   859,   869, ...[ 6499,  7945,  6499, ..., 11039, 10053,  8487]])
   1035         # change shape of arrays after fit (early-stopping or additional ests)
   1036         if n_stages != self.estimators_.shape[0]:
   1037             self.estimators_ = self.estimators_[:n_stages]
   1038             self.train_score_ = self.train_score_[:n_stages]

...........................................................................
C:\ProgramData\Anaconda3\envs\i3\lib\site-packages\sklearn\ensemble\gradient_boosting.py in _fit_stages(self=GradientBoostingRegressor(alpha=0.9, criterion='...le=1.0, verbose=0, warm_start=False), X=array([[ -1.00000000e+00,  -2.00000003e-01,  -1.....04051296e+09,   1.89901295e+01]], dtype=float32), y=memmap([ 27.,  35., -78., ...,   9.,  21.,  20.], dtype=float32), y_pred=array([[ 26.99980487], [ 34.99961777], ...], [ 21.00002664], [ 19.99974011]]), sample_weight=array([ 1.,  1.,  1., ...,  1.,  1.,  1.], dtype=float32), random_state=<mtrand.RandomState object>, begin_at_stage=0, mo
nitor=None, X_idx_sorted=array([[    0,  8619,     0, ...,   859,   869, ...[ 6499,  7945,  6499, ..., 11039, 10053,  8487]]))
   1084                                       sample_weight[~sample_mask])
   1085
   1086             # fit next stage of trees
   1087             y_pred = self._fit_stage(i, X, y, y_pred, sample_weight,
   1088                                      sample_mask, random_state, X_idx_sorted,
-> 1089                                      X_csc, X_csr)
        X_csc = None
        X_csr = None
   1090
   1091             # track deviance (= loss)
   1092             if do_oob:
   1093                 self.train_score_[i] = loss_(y[sample_mask],

...........................................................................
C:\ProgramData\Anaconda3\envs\i3\lib\site-packages\sklearn\ensemble\gradient_boosting.py in _fit_stage(self=GradientBoostingRegressor(alpha=0.9, criterion='...le=1.0, verbose=0, warm_start=False), i=122, X=array([[ -1.00000000e+00,  -2.00000003e-01,  -1.....04051296e+09,   1.89901295e+01]], dtype=float32), y=memmap([ 27., 35., -78., ...,   9.,  21.,  20.], dtype=float32), y_pred=array([[ 26.99980487], [ 34.99961777], ...], [ 21.00002664], [ 19.99974011]]), sample_weight=array([ 1.,  1.,  1., ...,  1.,  1.,  1.], dtype=float32), sample_mask=array([ True,  True,  True, ...,  True,  True,  True], dtype=bool), random_state=<mtrand.RandomState object>, X_idx_sorted=array([[    0,  8619,     0, ...,   859,   869, ...[ 6499,  7945,  6499, ..., 11039, 10
053,  8487]]), X_csc=None, X_csr=None)
    793                                              sample_weight, sample_mask,

    794                                              self.learning_rate, k=k)
    795             else:
    796                 loss.update_terminal_regions(tree.tree_, X, y, residual, y_pred,
    797                                              sample_weight, sample_mask,
--> 798                                              self.learning_rate, k=k)
        self.learning_rate = 0.28522087352060566
        k = 0
    799
    800             # add tree to ensemble
    801             self.estimators_[i, k] = tree
    802

...........................................................................
C:\ProgramData\Anaconda3\envs\i3\lib\site-packages\sklearn\ensemble\gradient_boosting.py in update_terminal_regions(self=<sklearn.ensemble.gradient_boosting.HuberLossFunction object>, tree=<sklearn.tree._tree.Tree object>, X=array([[ -1.00000000e+00,  -2.00000003e-01,  -1.....04051296e+09,   1.89901295e+01]], dtype=float32), y=memmap([ 27.,  35., -78., ...,   9.,  21.,  20.], dtype=float32), residual=array([  1.95126553e-04,   3.82230297e-04,   9.5...1895189e-06,  -2.66392852e-05,   2.59891505e-04]), y_pred=array([[ 26.99980487], [ 34.99961777], ...], [ 21.00002664], [ 19.99974011]]), sample_weight=array([ 1.,  1.,  1., ...,  1.,  1.,  1.], dtype=float32), sample_mask=array([ True,  True,  True, ...,  True,  True,  True], dtype=bool), learning_rate=0.28522087352060566, k=0)
    244
    245         # update each leaf (= perform line search)
    246         for leaf in np.where(tree.children_left == TREE_LEAF)[0]:
    247             self._update_terminal_region(tree, masked_terminal_regions,
    248                                          leaf, X, y, residual,
--> 249                                          y_pred[:, k], sample_weight)
        y_pred = array([[ 26.99980487], [ 34.99961777], ...], [ 21.00002664], [ 19.99974011]])
        k = 0
        sample_weight = array([ 1.,  1.,  1., ...,  1.,  1.,  1.], dtype=float32
)
    250
    251         # update predictions (both in-bag and out-of-bag)
    252         y_pred[:, k] += (learning_rate
    253                          * tree.value[:, 0, 0].take(terminal_regions, ax
is=0))

...........................................................................
C:\ProgramData\Anaconda3\envs\i3\lib\site-packages\sklearn\ensemble\gradient_boosting.py in _update_terminal_region(self=<sklearn.ensemble.gradient_boosting.HuberLossFunction object>, tree=<sklearn.tree._tree.Tree object>, terminal_regions=array([237, 237,  86, ..., 237, 237, 237], dtype=int64), leaf=260, X=array([[ -1.00000000e+00,  -2.00000003e-01,  -1.....04051296e+09,   1.89901295e+01]], dtype=float32), y=memmap([ 27.,  35., -78., ...,   9.,  21.,  20.], dtype=float32), residual=array([  1.95126553e-04,   3.82230297e-04,   9.5...1895189e-06,  -2.66392852e-05,   2.59891505e-04]), pred=array([ 26.99980487,  34.99961777, -78.03339049,...  9.00000162, 21.00002664,  19.99974011]), sample_weight=array([], dtype=float32))
    385         terminal_region = np.where(terminal_regions == leaf)[0]
    386         sample_weight = sample_weight.take(terminal_region, axis=0)
    387         gamma = self.gamma
    388         diff = (y.take(terminal_region, axis=0)
    389                 - pred.take(terminal_region, axis=0))
--> 390         median = _weighted_percentile(diff, sample_weight, percentile=50)
        median = undefined
        diff = array([], dtype=float64)
        sample_weight = array([], dtype=float32)
    391         diff_minus_median = diff - median
    392         tree.value[leaf, 0] = median + np.mean(
    393             np.sign(diff_minus_median) *
    394             np.minimum(np.abs(diff_minus_median), gamma))

...........................................................................
C:\ProgramData\Anaconda3\envs\i3\lib\site-packages\sklearn\utils\stats.py in _weighted_percentile(array=array([], dtype=float64), sample_weight=array([], dtype=float32), percentile=50)
     17     Compute the weighted ``percentile`` of ``array`` with ``sample_weight``.
     18     """
     19     sorted_idx = np.argsort(array)
     20
     21     # Find index of median prediction for each sample
---> 22     weight_cdf = stable_cumsum(sample_weight[sorted_idx])
        weight_cdf = undefined
        sample_weight = array([], dtype=float32)
        sorted_idx = array([], dtype=int64)
     23     percentile_idx = np.searchsorted(
     24         weight_cdf, (percentile / 100.) * weight_cdf[-1])
     25     return array[sorted_idx[percentile_idx]]

...........................................................................
C:\ProgramData\Anaconda3\envs\i3\lib\site-packages\sklearn\utils\extmath.py in stable_cumsum(arr=array([], dtype=float32), axis=None, rtol=1e-05, atol=1e-08)
    757     if np_version < (1, 9):
    758         return np.cumsum(arr, axis=axis, dtype=np.float64)
    759
    760     out = np.cumsum(arr, axis=axis, dtype=np.float64)
    761     expected = np.sum(arr, axis=axis, dtype=np.float64)
--> 762     if not np.all(np.isclose(out.take(-1, axis=axis), expected, rtol=rto
l,
        out.take = <built-in method take of numpy.ndarray object>
        axis = None
        expected = 0.0
        rtol = 1e-05
        atol = 1e-08
    763                              atol=atol, equal_nan=True)):
    764         warnings.warn('cumsum was found to be unstable: '
    765                       'its last element does not correspond to sum',
    766                       RuntimeWarning)

IndexError: cannot do a non-empty take from an empty axes.
```

Hints from issue discussion:
From that error message I can deduce that `_update_terminal_region` is being called with `terminal_regions` not containing `leaf`. I'm not familiar enough with this technique to know whether this is a state that should never happen, or whether this is a case that is unhandled but should be. (@pprett?) Either way, thanks for the bug report.
Maybe @glemaitre has a suggestion on this one?

At a first look, I have no clue. This is not really a part of the code that I know.
@jpeg729 do could you provide the data such that we can reproduce the bug and try to understand where this is coming from?
My dataset is private. Besides, it weighs several gigabytes. Furthermore, I don't think examining the dataset is going to enlighten you all that much, especially as I have recently noticed that changing the hyperparameters can sometimes cause the error to appear on a dataset that previously allowed training with no problem.

Based on @jnothman's observation I have added `if len(terminal_region) == 0: return`
after line 385: `terminal_region = np.where(terminal_regions == leaf)[0]` of [`...\sklearn\ensemble\gradient_boosting.py#L385`](https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/ensemble/gradient_boosting.py#L385)

I confirm that with this change I can train the model on my dataset without any errors. If I undo the modification then, as before, model.fit throws the above error.
                                                                                      I am not interesting by the dataset but by whatever can reproduce the bug. The easy thing is to have what you tried first. But this true that if you have a synthetic minimal example it is always better. 
If you have ideas for fixes I'll happily test them out on my dataset.
@jpeg729 we insist on having a non-regression test for bug fixes. That means that we will need to find a way to generate a small dataset that triggers this bug. It would be really valuable if you could post a stand-alone snippet reproducing the problem. You really are in the best position to do that because you are the one that managed to bump into the bug in the first place.
As I said, my dataset is several Gb and I have had no success finding a small subset that triggers the bug. Given the difficulty my dataset presents my guess is that if you try random datasets you will sooner or later come across a dataset that triggers the problem. I am running a script to search for random datasets that trigger the problem.

It may be more efficient to try to understand why empty leaves can occur and to engineer a dataset that triggers that eventuality.
They say that RandomForests and friends don't need data normalisation, but they forget to mention that you have to be careful with near infinite values.

This will provoke the error 9 times out of 10.
```python
import numpy as np
from sklearn.ensemble import GradientBoostingRegressor

data = np.random.randn(100, 100)*1e38
data = np.nan_to_num(data.astype('float32'))
X = data[:, :-1]
y = data[:, -1]
GradientBoostingRegressor(loss="huber").fit(X, y)
```

As far as I can see, the only way a GradientBoostingRegressor can have an empty leaf is if the DecisionTreeRegressor returns a tree with an empty leaf, and sure enough running the following code often reveals the presence of several empty leaves.
```python
import numpy as np
from sklearn.tree import DecisionTreeRegressor
from sklearn.tree._tree import TREE_LEAF

data = np.random.randn(100, 100)*1e38
data = np.nan_to_num(data.astype('float32'))
X = data[:, :-1]
y = data[:, -1]
tree = DecisionTreeRegressor().fit(X,y)
terminal_regions = tree.apply(X)
count_empty = 0
for leaf in np.where(tree.tree_.children_left == TREE_LEAF)[0]:
    count_empty += len(np.where(terminal_regions == leaf)[0]) == 0

print(count_empty)
```

### Possible solutions
* Fix DecisionTreeRegressor so that it never returns a tree with an empty leaf even if some of its input values are near infinite.
* Add a warning to the docs telling people to make sure their dataset contains no near infinite values.
OK thanks a lot @jpeg729! Tweaking your snippet just a bit, here is a stand-alone deterministic snippet reproducing the problem:

```py
import numpy as np
from sklearn.ensemble import GradientBoostingRegressor

rng = np.random.RandomState(42)
data = rng.randn(100, 100)*1e38
data = np.nan_to_num(data.astype('float32'))
X = data[:, :-1]
y = data[:, -1]
GradientBoostingRegressor(loss="huber", random_state=rng).fit(X, y)
```

Commit message: add line in whats new

Patch to review:
diff --git a/doc/whats_new/v0.20.rst b/doc/whats_new/v0.20.rst
index 7e7d39dbf1759..72294bb51bdc7 100644
--- a/doc/whats_new/v0.20.rst
+++ b/doc/whats_new/v0.20.rst
@@ -215,6 +215,11 @@ Classifiers and regressors
   unexpected TypeError given dense inputs.
   :issue:`10412` by :user:`Jiongyan Zhang <qmick>`.
 
+- Fixed a bug in :class:`tree._splitter.BestSplitter` and
+  :class:`tree._splitter.BestSparseSplitter` where threshold could become
+  infinite when values in X were near infinite.
+  :issue:`10536` by :user:`Jonathan Ohayon <Johayon>`.
+
 Decomposition, manifold learning and clustering
 
 - Fix for uninformative error in :class:`decomposition.IncrementalPCA`:
diff --git a/sklearn/tree/_splitter.pyx b/sklearn/tree/_splitter.pyx
index 4e9b619b11e47..6bbf6986c8b41 100644
--- a/sklearn/tree/_splitter.pyx
+++ b/sklearn/tree/_splitter.pyx
@@ -481,9 +481,12 @@ cdef class BestSplitter(BaseDenseSplitter):
 
                             if current_proxy_improvement > best_proxy_improvement:
                                 best_proxy_improvement = current_proxy_improvement
-                                current.threshold = (Xf[p - 1] + Xf[p]) / 2.0
+                                # sum of halves is used to avoid infinite value
+                                current.threshold = Xf[p - 1] / 2.0 + Xf[p] / 2.0
 
-                                if current.threshold == Xf[p]:
+                                if ((current.threshold == Xf[p]) or
+                                    (current.threshold == INFINITY) or
+                                    (current.threshold == -INFINITY)):
                                     current.threshold = Xf[p - 1]
 
                                 best = current  # copy
@@ -1382,9 +1385,12 @@ cdef class BestSparseSplitter(BaseSparseSplitter):
 
                             if current_proxy_improvement > best_proxy_improvement:
                                 best_proxy_improvement = current_proxy_improvement
+                                # sum of halves used to avoid infinite values
+                                current.threshold = Xf[p_prev] / 2.0 + Xf[p] / 2.0
 
-                                current.threshold = (Xf[p_prev] + Xf[p]) / 2.0
-                                if current.threshold == Xf[p]:
+                                if ((current.threshold == Xf[p]) or
+                                    (current.threshold == INFINITY) or
+                                    (current.threshold == -INFINITY)):
                                     current.threshold = Xf[p_prev]
 
                                 best = current
diff --git a/sklearn/tree/tests/test_tree.py b/sklearn/tree/tests/test_tree.py
index 53767a65ed891..8c93a3b1d8580 100644
--- a/sklearn/tree/tests/test_tree.py
+++ b/sklearn/tree/tests/test_tree.py
@@ -3,6 +3,7 @@
 """
 import copy
 import pickle
+import pytest
 from functools import partial
 from itertools import product
 import struct
@@ -1731,3 +1732,19 @@ def _pickle_copy(obj):
             assert_equal(typename, typename_)
             assert_equal(n_outputs, n_outputs_)
             assert_equal(n_samples, n_samples_)
+
+
+@pytest.mark.parametrize('data',
+                         np.random.RandomState(0).randn(50, 100, 11) * 2e38)
+def test_empty_leaf_infinite_threshold(data):
+    # try to make empty leaf by using near infinite value.
+    data = np.nan_to_num(data.astype('float32'))
+    X = data[:, :-1]
+    y = data[:, -1]
+    tree = DecisionTreeRegressor().fit(X, y)
+    terminal_regions = tree.apply(X)
+    empty_leaf = set(np.where(tree.tree_.children_left
+                              == TREE_LEAF)[0]).difference(terminal_regions)
+    infinite_threshold = np.where(~np.isfinite(tree.tree_.threshold))[0]
+    assert (len(infinite_threshold) == 0)
+    assert (len(empty_leaf) == 0)


Code review comments:
diff_hunk: @@ -215,6 +215,11 @@ Classifiers and regressors
   unexpected TypeError given dense inputs.
   :issue:`10412` by :user:`Jiongyan Zhang <qmick>`.
 
+- Fixed a bug in :class:`tree._splitter.BestSplitter` and
path: doc/whats_new/v0.20.rst
text: @user1:
Please refer to it by public API elements, i.e. without reference to `_splitter`

@author:
Should I refer to `tree.DecisionTreeClassifier` and `tree.DecisionTreeRegressor` ?

@user1:
You can just say decision trees, which covers forests, etc. Mention the splitter parameter