[MRG+1] fused types in isotonic_regression (#9106)

vene · agramfort · commit 4e56b293d66f · 2019-02-26T15:17:14.000+01:00
* ENH fused types in isotonic_regression

* make X drive computation dtype in IsotonicRegression

* preserve current behaviour if transform w/o fit

* thoroughly check sample weights; avoid storing dtype explicitly

* consistent testing and behavior

* misc

* update what's new

* fix for interp1d upcast on old scipy

* flake8 remove blank line

sorry, i'm rusty

* add failing test

* FIX dtype bug in _make_unique
diff --git a/doc/whats_new/v0.21.rst b/doc/whats_new/v0.21.rst
@@ -328,6 +328,13 @@ Support for Python 3.4 and below has been officially dropped.
   :user:`Jérémie du Boisberranger <jeremiedbb>`.
 
 
+:mod:`sklearn.isotonic`
+.......................
+
+- |Feature| Allow different dtypes (such as float32) in
+  :class:`isotonic.IsotonicRegression` :issue:`8769` by :user:`Vlad Niculae <vene>`
+
+
 Multiple modules
 ................
 
diff --git a/sklearn/_isotonic.pyx b/sklearn/_isotonic.pyx
@@ -7,17 +7,16 @@
 import numpy as np
 cimport numpy as np
 cimport cython
-
-ctypedef np.float64_t DOUBLE
+from cython cimport floating
 
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
 @cython.cdivision(True)
-def _inplace_contiguous_isotonic_regression(DOUBLE[::1] y, DOUBLE[::1] w):
+def _inplace_contiguous_isotonic_regression(floating[::1] y, floating[::1] w):
     cdef:
         Py_ssize_t n = y.shape[0], i, k
-        DOUBLE prev_y, sum_wy, sum_w
+        floating prev_y, sum_wy, sum_w
         Py_ssize_t[::1] target = np.arange(n, dtype=np.intp)
 
     # target describes a list of blocks.  At any time, if [i..j] (inclusive) is
@@ -68,9 +67,9 @@ def _inplace_contiguous_isotonic_regression(DOUBLE[::1] y, DOUBLE[::1] w):
 @cython.boundscheck(False)
 @cython.wraparound(False)
 @cython.cdivision(True)
-def _make_unique(np.ndarray[dtype=np.float64_t] X,
-                 np.ndarray[dtype=np.float64_t] y,
-                 np.ndarray[dtype=np.float64_t] sample_weights):
+def _make_unique(np.ndarray[dtype=floating] X,
+                 np.ndarray[dtype=floating] y,
+                 np.ndarray[dtype=floating] sample_weights):
     """Average targets for duplicate X, drop duplicates.
 
     Aggregates duplicate X values into a single X value where
@@ -82,18 +81,20 @@ def _make_unique(np.ndarray[dtype=np.float64_t] X,
     unique_values = len(np.unique(X))
     if unique_values == len(X):
         return X, y, sample_weights
-    cdef np.ndarray[dtype=np.float64_t] y_out = np.empty(unique_values)
-    cdef np.ndarray[dtype=np.float64_t] x_out = np.empty(unique_values)
-    cdef np.ndarray[dtype=np.float64_t] weights_out = np.empty(unique_values)
 
-    cdef np.float64_t current_x = X[0]
-    cdef np.float64_t current_y = 0
-    cdef np.float64_t current_weight = 0
-    cdef np.float64_t y_old = 0
+    cdef np.ndarray[dtype=floating] y_out = np.empty(unique_values,
+                                                     dtype=X.dtype)
+    cdef np.ndarray[dtype=floating] x_out = np.empty_like(y_out)
+    cdef np.ndarray[dtype=floating] weights_out = np.empty_like(y_out)
+
+    cdef floating current_x = X[0]
+    cdef floating current_y = 0
+    cdef floating current_weight = 0
+    cdef floating y_old = 0
     cdef int i = 0
     cdef int current_count = 0
     cdef int j
-    cdef np.float64_t x
+    cdef floating x
     cdef int n_samples = len(X)
     for j in range(n_samples):
         x = X[j]
diff --git a/sklearn/isotonic.py b/sklearn/isotonic.py
@@ -119,11 +119,12 @@ def isotonic_regression(y, sample_weight=None, y_min=None, y_max=None,
     by Michael J. Best and Nilotpal Chakravarti, section 3.
     """
     order = np.s_[:] if increasing else np.s_[::-1]
-    y = np.array(y[order], dtype=np.float64)
+    y = as_float_array(y)
+    y = np.array(y[order], dtype=y.dtype)
     if sample_weight is None:
-        sample_weight = np.ones(len(y), dtype=np.float64)
+        sample_weight = np.ones(len(y), dtype=y.dtype)
     else:
-        sample_weight = np.array(sample_weight[order], dtype=np.float64)
+        sample_weight = np.array(sample_weight[order], dtype=y.dtype)
 
     _inplace_contiguous_isotonic_regression(y, sample_weight)
     if y_min is not None or y_max is not None:
@@ -239,10 +240,9 @@ def _build_f(self, X, y):
 
     def _build_y(self, X, y, sample_weight, trim_duplicates=True):
         """Build the y_ IsotonicRegression."""
+        X = as_float_array(X)
+        y = check_array(y, dtype=X.dtype, ensure_2d=False)
         check_consistent_length(X, y, sample_weight)
-        X, y = [check_array(x, ensure_2d=False) for x in [X, y]]
-
-        y = as_float_array(y)
         self._check_fit_data(X, y, sample_weight)
 
         # Determine increasing if auto-determination requested
@@ -254,15 +254,15 @@ def _build_y(self, X, y, sample_weight, trim_duplicates=True):
         # If sample_weights is passed, removed zero-weight values and clean
         # order
         if sample_weight is not None:
-            sample_weight = check_array(sample_weight, ensure_2d=False)
+            sample_weight = check_array(sample_weight, ensure_2d=False,
+                                        dtype=X.dtype)
             mask = sample_weight > 0
             X, y, sample_weight = X[mask], y[mask], sample_weight[mask]
         else:
-            sample_weight = np.ones(len(y))
+            sample_weight = np.ones(len(y), dtype=X.dtype)
 
         order = np.lexsort((y, X))
-        X, y, sample_weight = [array[order].astype(np.float64, copy=False)
-                               for array in [X, y, sample_weight]]
+        X, y, sample_weight = [array[order] for array in [X, y, sample_weight]]
         unique_X, unique_y, unique_sample_weight = _make_unique(
             X, y, sample_weight)
 
@@ -345,7 +345,14 @@ def transform(self, T):
         T_ : array, shape=(n_samples,)
             The transformed data
         """
-        T = as_float_array(T)
+
+        if hasattr(self, '_necessary_X_'):
+            dtype = self._necessary_X_.dtype
+        else:
+            dtype = np.float64
+
+        T = check_array(T, dtype=dtype, ensure_2d=False)
+
         if len(T.shape) != 1:
             raise ValueError("Isotonic regression input should be a 1d array")
 
@@ -357,7 +364,13 @@ def transform(self, T):
 
         if self.out_of_bounds == "clip":
             T = np.clip(T, self.X_min_, self.X_max_)
-        return self.f_(T)
+
+        res = self.f_(T)
+
+        # on scipy 0.17, interp1d up-casts to float64, so we cast back
+        res = res.astype(T.dtype)
+
+        return res
 
     def predict(self, T):
         """Predict new data by linear interpolation.
diff --git a/sklearn/tests/test_isotonic.py b/sklearn/tests/test_isotonic.py
@@ -4,8 +4,9 @@
 import copy
 
 from sklearn.isotonic import (check_increasing, isotonic_regression,
-                              IsotonicRegression)
+                              IsotonicRegression, _make_unique)
 
+from sklearn.utils.validation import as_float_array
 from sklearn.utils.testing import (assert_raises, assert_array_equal,
                                    assert_equal,
                                    assert_array_almost_equal,
@@ -459,3 +460,32 @@ def test_isotonic_copy_before_fit():
     # https://github.com/scikit-learn/scikit-learn/issues/6628
     ir = IsotonicRegression()
     copy.copy(ir)
+
+
+def test_isotonic_dtype():
+    y = [2, 1, 4, 3, 5]
+    weights = np.array([.9, .9, .9, .9, .9], dtype=np.float64)
+    reg = IsotonicRegression()
+
+    for dtype in (np.int32, np.int64, np.float32, np.float64):
+        for sample_weight in (None, weights.astype(np.float32), weights):
+            y_np = np.array(y, dtype=dtype)
+            expected_dtype = as_float_array(y_np).dtype
+
+            res = isotonic_regression(y_np, sample_weight=sample_weight)
+            assert_equal(res.dtype, expected_dtype)
+
+            X = np.arange(len(y)).astype(dtype)
+            reg.fit(X, y_np, sample_weight=sample_weight)
+            res = reg.predict(X)
+            assert_equal(res.dtype, expected_dtype)
+
+
+def test_make_unique_dtype():
+    x_list = [2, 2, 2, 3, 5]
+    for dtype in (np.float32, np.float64):
+        x = np.array(x_list, dtype=dtype)
+        y = x.copy()
+        w = np.ones_like(x)
+        x, y, w = _make_unique(x, y, w)
+        assert_array_equal(x, [2, 3, 5])
