Repository: civisanalytics/civisml-extensions

Issue:
`_base_est_fit_predict` looks wrong
I believe cv.split(X, y) might give different results when called at different times. There might be randomness involved like shuffling, etc.

This is problematic because `ybase`, which is based on the old `train` & `test` but not coincide with with the new `test` and therefore `y[test]`.

old `train` & `test`:

```
for train, test in cv.split(X, y):
    for name, est in self.estimator_list[:-1]:
        # adapted from sklearn.model_selection._fit_and_predict
        # Adjust length of sample weights
        fit_params_est_adjusted = dict([
            (k, _index_param_value(X, v, train))
            for k, v in fit_params_ests[name].items()])

        # Fit estimator on training set and score out-of-sample
        _jobs.append(delayed(_fit_predict)(
            clone(est),
            X[train],
            y[train],
            X[test],
            **fit_params_est_adjusted))
```

new `train` & `test`:

```
# Extract the results from joblib
Xmeta, ymeta = None, None
for train, test in cv.split(X, y):
    ybase = np.empty((y[test].shape[0], 0))
    for name, est in self.estimator_list[:-1]:
        # Build design matrix out of out-of-sample predictions
        ybase = np.hstack((ybase, _out.pop(0)))

    # Append the test outputs to what will eventually be the features
    # for the meta-estimator.
    if Xmeta is not None:
        ymeta = np.concatenate((ymeta, y[test]))
        Xmeta = np.vstack((Xmeta, ybase))
    else:
        Xmeta = ybase
        ymeta = y[test]
```

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6344c7c..c527285 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,11 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 
 ## Unreleased
 
+### Fixed
+- Now caching CV indices. When CV generators are passed with `shuffle=True` and
+  no `random_state` is set, they produce different CV folds on each call to
+  `split` (#22).
+
 ### Fixed
 - Updated `scipy` dependency in `requirements.txt` file to `scipy>=0.14,<2.0`
 - ``DataFrameETL`` now correctly handles all ``Categorial``-type columns
diff --git a/civismlext/stacking.py b/civismlext/stacking.py
index c06a986..8c04a47 100644
--- a/civismlext/stacking.py
+++ b/civismlext/stacking.py
@@ -233,20 +233,25 @@ class BaseStackedModel(BaseEstimator):
         out-of-sample predictions on the test folds as features for the
         meta-estimator. Also return the fit_params for the meta-estimator.
         """
-        # TO DO:
-        # -use joblib to split base estimator training up
-
         y = y.squeeze()
         # Construct CV iterator
         cv = self._check_cv(y=y)
+        # Extract CV indices since we need them twice, and un-seeded CV
+        # generators with `shuffle=True` split differently each time.
+        train_inds = []
+        test_inds = []
+        for train, test in cv.split(X, y):
+            train_inds.append(train)
+            test_inds.append(test)
 
         fit_params_ests = self._extract_fit_params(**fit_params)
         _fit_predict = self._get_fit_predict_function()
 
+        _jobs = []
+
         # Loop over CV folds to get out-of-sample predictions, which become the
         # features for the meta-estimator.
-        _jobs = []
-        for train, test in cv.split(X, y):
+        for train, test in zip(train_inds, test_inds):
             for name, est in self.estimator_list[:-1]:
                 # adapted from sklearn.model_selection._fit_and_predict
                 # Adjust length of sample weights
@@ -269,7 +274,7 @@ class BaseStackedModel(BaseEstimator):
 
         # Extract the results from joblib
         Xmeta, ymeta = None, None
-        for train, test in cv.split(X, y):
+        for test in test_inds:
             ybase = np.empty((y[test].shape[0], 0))
             for name, est in self.estimator_list[:-1]:
                 # Build design matrix out of out-of-sample predictions


Test patch:
diff --git a/civismlext/test/test_stacking.py b/civismlext/test/test_stacking.py
index eefd2dc..39741c6 100644
--- a/civismlext/test/test_stacking.py
+++ b/civismlext/test/test_stacking.py
@@ -8,7 +8,7 @@ from ..nonnegative import NonNegativeLinearRegression
 import pytest
 import numpy as np
 from sklearn.tree import DecisionTreeClassifier
-from sklearn.model_selection import train_test_split
+from sklearn.model_selection import KFold, train_test_split
 from sklearn.metrics import roc_curve, auc, mean_squared_error
 from sklearn.linear_model import LogisticRegression, LinearRegression
 from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
@@ -729,6 +729,32 @@ def test_fit_regression():
     assert sr.meta_estimator.fit_params == {'bar': 'b'}
 
 
+def test_cv_shuffle_indices():
+    """Make sure xmeta and ymeta retain the correct order, even when the CV
+    generator is shuffling. This is checking for the bug reported in issue #16.
+    """
+    estlist = [('be1', PassThruReg()),
+               ('be2', PassThruReg()),
+               ('meta', PassThruReg())]
+    sr = StackedRegressor(estlist, cv=KFold(n_splits=2, shuffle=True))
+    x = np.arange(6)
+    y = np.arange(6)
+
+    # Suppose the train indices of a 2-fold CV are:
+    #  [a, b, c]  and  [d, e, f]
+    # Then the test indices are:
+    #  [d, e, f]  and  [a, b, c]
+    # Since xmeta is just a pass-through of x[train] (horizontally stacked
+    # twice, due to the two base estimators) and ymeta is a pass-through of
+    # y[test], and since x = y, we should expect that:
+    #  xmeta[inds, 0] == ymeta
+    #
+    inds = np.array([3, 4, 5, 0, 1, 2])
+
+    xmeta, ymeta, _ = sr._base_est_fit_predict(x, y)
+    np.testing.assert_equal(xmeta[inds, 0], ymeta)
+
+
 @mock.patch('civismlext.stacking.clone', lambda x: x)
 def test_base_est_fit_predict_multiout_regression():
     estlist = [('be1', PassThruReg()),