diff --git a/.binder/requirements.txt b/.binder/requirements.txt
index 9ecc5c6fba79c..1e46398a5b7b0 100644
--- a/.binder/requirements.txt
+++ b/.binder/requirements.txt
@@ -1,8 +1,5 @@
---extra-index https://pypi.anaconda.org/scipy-wheels-nightly/simple scikit-learn
---pre
 matplotlib
 scikit-image
 pandas
 sphinx-gallery
-scikit-learn
-
+scikit-learn>=0.23,<0.24
diff --git a/build_tools/azure/test_script.cmd b/build_tools/azure/test_script.cmd
index f1e516d81fd99..2078692b5dc23 100644
--- a/build_tools/azure/test_script.cmd
+++ b/build_tools/azure/test_script.cmd
@@ -10,7 +10,8 @@ mkdir %TMP_FOLDER%
 cd %TMP_FOLDER%
 
 if "%CHECK_WARNINGS%" == "true" (
-    set PYTEST_ARGS=%PYTEST_ARGS% -Werror::DeprecationWarning -Werror::FutureWarning
+    REM numpy's 1.19.0's tostring() deprecation is ignored until scipy and joblib removes its usage
+    set PYTEST_ARGS=%PYTEST_ARGS% -Werror::DeprecationWarning -Werror::FutureWarning -Wignore:tostring:DeprecationWarning
 )
 
 if "%COVERAGE%" == "true" (
diff --git a/build_tools/azure/test_script.sh b/build_tools/azure/test_script.sh
index 77a950d86549c..adef606df1d94 100755
--- a/build_tools/azure/test_script.sh
+++ b/build_tools/azure/test_script.sh
@@ -29,7 +29,8 @@ if [[ "$COVERAGE" == "true" ]]; then
 fi
 
 if [[ -n "$CHECK_WARNINGS" ]]; then
-    TEST_CMD="$TEST_CMD -Werror::DeprecationWarning -Werror::FutureWarning"
+    # numpy's 1.19.0's tostring() deprecation is ignored until scipy and joblib removes its usage
+    TEST_CMD="$TEST_CMD -Werror::DeprecationWarning -Werror::FutureWarning -Wignore:tostring:DeprecationWarning"
 fi
 
 if [[ "$PYTHON_VERSION" == "*" ]]; then
diff --git a/conftest.py b/conftest.py
index 874931341e195..f52ee5f19634d 100644
--- a/conftest.py
+++ b/conftest.py
@@ -7,7 +7,6 @@
 
 import platform
 import sys
-from distutils.version import LooseVersion
 import os
 
 import pytest
@@ -17,10 +16,11 @@
 from sklearn.utils import _IS_32BIT
 from sklearn.externals import _pilutil
 from sklearn._build_utils.deprecated_modules import _DEPRECATED_MODULES
+from sklearn.utils.fixes import np_version, parse_version
 
 PYTEST_MIN_VERSION = '3.3.0'
 
-if LooseVersion(pytest.__version__) < PYTEST_MIN_VERSION:
+if parse_version(pytest.__version__) < parse_version(PYTEST_MIN_VERSION):
     raise ImportError('Your version of pytest is too old, you should have '
                       'at least pytest >= {} installed.'
                       .format(PYTEST_MIN_VERSION))
@@ -54,8 +54,7 @@ def pytest_collection_modifyitems(config, items):
     # run doctests only for numpy >= 1.14.
     skip_doctests = False
     try:
-        import numpy as np
-        if LooseVersion(np.__version__) < LooseVersion('1.14'):
+        if np_version < parse_version('1.14'):
             reason = 'doctests are only run for numpy >= 1.14'
             skip_doctests = True
         elif _IS_32BIT:
diff --git a/doc/conf.py b/doc/conf.py
index 74e37d01307be..539229a605330 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -80,7 +80,7 @@
 
 # General information about the project.
 project = 'scikit-learn'
-copyright = '2007 - 2019, scikit-learn developers (BSD License)'
+copyright = '2007 - 2020, scikit-learn developers (BSD License)'
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
diff --git a/doc/modules/linear_model.rst b/doc/modules/linear_model.rst
index 477baca9c4de3..c4767d0cb2d64 100644
--- a/doc/modules/linear_model.rst
+++ b/doc/modules/linear_model.rst
@@ -997,8 +997,8 @@ Examples of use cases include:
 * Risk modeling / insurance policy pricing:  number of claim events /
   policyholder per year (Poisson), cost per event (Gamma), total cost per
   policyholder per year (Tweedie / Compound Poisson Gamma).
-* Predictive maintenance: number of production interruption events per year:
-  Poisson, duration of interruption: Gamma, total interruption time per year
+* Predictive maintenance: number of production interruption events per year
+  (Poisson), duration of interruption (Gamma), total interruption time per year
   (Tweedie / Compound Poisson Gamma).
 
 
diff --git a/doc/modules/model_persistence.rst b/doc/modules/model_persistence.rst
index 81b3dfce63b99..3c59b7a313a17 100644
--- a/doc/modules/model_persistence.rst
+++ b/doc/modules/model_persistence.rst
@@ -20,7 +20,7 @@ Persistence example
 -------------------
 
 It is possible to save a model in scikit-learn by using Python's built-in
-persistence model, namely `pickle <https://docs.python.org/2/library/pickle.html>`_::
+persistence model, namely `pickle <https://docs.python.org/3/library/pickle.html>`_::
 
   >>> from sklearn import svm
   >>> from sklearn import datasets
diff --git a/doc/templates/index.html b/doc/templates/index.html
index f49fbc2f4c540..58ba61a2c0ca3 100644
--- a/doc/templates/index.html
+++ b/doc/templates/index.html
@@ -155,6 +155,10 @@ <h4 class="sk-landing-call-header">News</h4>
         <ul class="sk-landing-call-list list-unstyled">
         <li><strong>On-going development:</strong>
         <a href="https://scikit-learn.org/dev/whats_new.html"><strong>What's new</strong> (Changelog)</a>
+        <li><strong>August 2020.</strong> scikit-learn 0.23.2 is available for download (<a href="whats_new/v0.23.html#version-0-23-2">Changelog</a>).
+        </li>
+        <li><strong>May 2020.</strong> scikit-learn 0.23.1 is available for download (<a href="whats_new/v0.23.html#version-0-23-1">Changelog</a>).
+        </li>
         <li><strong>May 2020.</strong> scikit-learn 0.23.0 is available for download (<a href="whats_new/v0.23.html#version-0-23-0">Changelog</a>).
         </li>
         <li><strong>Scikit-learn from 0.23 requires Python 3.6 or greater.</strong>
diff --git a/doc/themes/scikit-learn-modern/static/css/theme.css b/doc/themes/scikit-learn-modern/static/css/theme.css
index ceda27c6de093..b22b700736f56 100644
--- a/doc/themes/scikit-learn-modern/static/css/theme.css
+++ b/doc/themes/scikit-learn-modern/static/css/theme.css
@@ -83,12 +83,12 @@ span.highlighted {
 }
 
 div.highlight {
-  padding: 0.2rem 0.5rem;
   border: 1px solid #ddd;
   margin-bottom: 1rem;
 }
 
 div.highlight pre {
+  padding: 0.2rem 0.5rem;
   margin-bottom: 0;
   line-height: 1.2rem;
 }
diff --git a/doc/whats_new/v0.23.rst b/doc/whats_new/v0.23.rst
index b089b4d5a1d78..2ba8eb2797e5e 100644
--- a/doc/whats_new/v0.23.rst
+++ b/doc/whats_new/v0.23.rst
@@ -2,6 +2,122 @@
 
 .. currentmodule:: sklearn
 
+.. _changes_0_23_2:
+
+Version 0.23.2
+==============
+
+**August 3 2020**
+
+Changed models
+--------------
+
+The following estimators and functions, when fit with the same data and
+parameters, may produce different models from the previous version. This often
+occurs due to changes in the modelling logic (bug fixes or enhancements), or in
+random sampling procedures.
+
+- |Fix| ``inertia_`` attribute of :class:`cluster.KMeans` and
+  :class:`cluster.MiniBatchKMeans`.
+
+Details are listed in the changelog below.
+
+(While we are trying to better inform users by providing this information, we
+cannot assure that this list is complete.)
+
+Changelog
+---------
+
+:mod:`sklearn.cluster`
+......................
+
+- |Fix| Fixed a bug in :class:`cluster.KMeans` where rounding errors could
+  prevent convergence to be declared when `tol=0`. :pr:`17959` by
+  :user:`Jérémie du Boisberranger <jeremiedbb>`.
+
+- |Fix| Fixed a bug in :class:`cluster.KMeans` and
+  :class:`cluster.MiniBatchKMeans` where the reported inertia was incorrectly
+  weighted by the sample weights. :pr:`17848` by
+  :user:`Jérémie du Boisberranger <jeremiedbb>`.
+
+- |Fix| Fixed a bug in :class:`cluster.MeanShift` with `bin_seeding=True`. When
+  the estimated bandwidth is 0, the behavior is equivalent to
+  `bin_seeding=False`.
+  :pr:`17742` by :user:`Jeremie du Boisberranger <jeremiedbb>`.
+
+- |Fix| Fixed a bug in :class:`cluster.AffinityPropagation`, that
+  gives incorrect clusters when the array dtype is float32.
+  :pr:`17995` by :user:`Thomaz Santana  <Wikilicious>` and
+  :user:`Amanda Dsouza <amy12xx>`.
+
+:mod:`sklearn.decomposition`
+............................
+
+- |Fix| Fixed a bug in
+  :func:`decomposition.MiniBatchDictionaryLearning.partial_fit` which should
+  update the dictionary by iterating only once over a mini-batch.
+  :pr:`17433` by :user:`Chiara Marmo <cmarmo>`.
+
+- |Fix| Avoid overflows on Windows in
+  :func:`decomposition.IncrementalPCA.partial_fit` for large ``batch_size`` and
+  ``n_samples`` values.
+  :pr:`17985` by :user:`Alan Butler <aldee153>` and
+  :user:`Amanda Dsouza <amy12xx>`.
+
+:mod:`sklearn.ensemble`
+.......................
+
+- |Fix| Fixed bug in :class:`ensemble.MultinomialDeviance` where the
+  average of logloss was incorrectly calculated as sum of logloss.
+  :pr:`17694` by :user:`Markus Rempfler <rempfler>` and
+  :user:`Tsutomu Kusanagi <t-kusanagi2>`.
+
+- |Fix| Fixes :class:`ensemble.StackingClassifier` and
+  :class:`ensemble.StackingRegressor` compatibility with estimators that
+  do not define `n_features_in_`. :pr:`17357` by `Thomas Fan`_.
+
+:mod:`sklearn.feature_extraction`
+.................................
+
+- |Fix| Fixes bug in :class:`feature_extraction.text.CountVectorizer` where
+  sample order invariance was broken when `max_features` was set and features
+  had the same count. :pr:`18016` by `Thomas Fan`_, `Roman Yurchak`_, and
+  `Joel Nothman`_.
+
+:mod:`sklearn.linear_model`
+...........................
+
+- |Fix| :func:`linear_model.lars_path` does not overwrite `X` when
+  `X_copy=True` and `Gram='auto'`. :pr:`17914` by `Thomas Fan`_.
+
+:mod:`sklearn.manifold`
+.......................
+
+- |Fix| Fixed a bug where :func:`metrics.pairwise_distances` would raise an
+  error if ``metric='seuclidean'`` and ``X`` is not type ``np.float64``.
+  :pr:`15730` by :user:`Forrest Koch <ForrestCKoch>`.
+
+:mod:`sklearn.metrics`
+......................
+
+- |Fix| Fixed a bug in :func:`metrics.mean_squared_error` where the
+  average of multiple RMSE values was incorrectly calculated as the root of the
+  average of multiple MSE values.
+  :pr:`17309` by :user:`Swier Heeres <swierh>`.
+
+:mod:`sklearn.pipeline`
+.......................
+
+- |Fix| :class:`pipeline.FeatureUnion` raises a deprecation warning when
+  `None` is included in `transformer_list`. :pr:`17360` by `Thomas Fan`_.
+
+:mod:`sklearn.utils`
+....................
+
+- |Fix| Fix :func:`utils.estimator_checks.check_estimator` so that all test
+  cases support the `binary_only` estimator tag.
+  :pr:`17812` by :user:`Bruno Charron <brcharron>`.
+
 .. _changes_0_23_1:
 
 Version 0.23.1
@@ -23,6 +139,7 @@ Changelog
   provided by the user were modified in place. :pr:`17204` by
   :user:`Jeremie du Boisberranger <jeremiedbb>`.
 
+
 Miscellaneous
 .............
 
@@ -44,8 +161,6 @@ refer to
 
 .. include:: changelog_legend.inc
 
-Put the changes in their relevant module.
-
 Enforcing keyword-only arguments
 --------------------------------
 
@@ -171,7 +286,7 @@ Changelog
   deprecated. It has no effect. :pr:`11950` by
   :user:`Jeremie du Boisberranger <jeremiedbb>`.
 
-- |API| The ``random_state`` parameter has been added to 
+- |API| The ``random_state`` parameter has been added to
   :class:`cluster.AffinityPropagation`. :pr:`16801` by :user:`rcwoolston`
   and :user:`Chiara Marmo <cmarmo>`.
 
@@ -364,7 +479,7 @@ Changelog
   for each feature. :pr:`16403` by :user:`Narendra Mukherjee <narendramukherjee>`.
 
 - |Enhancement| :class:`impute.SimpleImputer`, :class:`impute.KNNImputer`, and
-  :class:`impute.SimpleImputer` accepts pandas' nullable integer dtype with
+  :class:`impute.IterativeImputer` accepts pandas' nullable integer dtype with
   missing values. :pr:`16508` by `Thomas Fan`_.
 
 :mod:`sklearn.inspection`
@@ -467,7 +582,7 @@ Changelog
   an error when `y_true` and `y_pred` were length zero and `labels` was
   not `None`. In addition, we raise an error when an empty list is given to
   the `labels` parameter.
-  :pr:`16442` by `Kyle Parsons <parsons-kyle-89>`.
+  :pr:`16442` by :user:`Kyle Parsons <parsons-kyle-89>`.
 
 - |API| Changed the formatting of values in
   :meth:`metrics.ConfusionMatrixDisplay.plot` and
@@ -491,7 +606,7 @@ Changelog
   :pr:`15622` by :user:`Gregory Morse <GregoryMorse>`.
 
 - |Fix| :func:`model_selection.cross_val_predict` supports
-  `method="predict_proba"` when `y=None`.:pr:`15918` by
+  `method="predict_proba"` when `y=None`. :pr:`15918` by
   :user:`Luca Kubin <lkubin>`.
 
 - |Fix| :func:`model_selection.fit_grid_point` is deprecated in 0.23 and will
diff --git a/examples/cluster/plot_coin_segmentation.py b/examples/cluster/plot_coin_segmentation.py
index bda1d717b2479..c662276459f92 100644
--- a/examples/cluster/plot_coin_segmentation.py
+++ b/examples/cluster/plot_coin_segmentation.py
@@ -25,7 +25,6 @@
 import time
 
 import numpy as np
-from distutils.version import LooseVersion
 from scipy.ndimage.filters import gaussian_filter
 import matplotlib.pyplot as plt
 import skimage
@@ -34,9 +33,10 @@
 
 from sklearn.feature_extraction import image
 from sklearn.cluster import spectral_clustering
+from sklearn.utils.fixes import parse_version
 
 # these were introduced in skimage-0.14
-if LooseVersion(skimage.__version__) >= '0.14':
+if parse_version(skimage.__version__) >= parse_version('0.14'):
     rescale_params = {'anti_aliasing': False, 'multichannel': False}
 else:
     rescale_params = {}
diff --git a/examples/cluster/plot_coin_ward_segmentation.py b/examples/cluster/plot_coin_ward_segmentation.py
index b45881e08c426..218e6fb4c8d86 100644
--- a/examples/cluster/plot_coin_ward_segmentation.py
+++ b/examples/cluster/plot_coin_ward_segmentation.py
@@ -17,7 +17,6 @@
 import time as time
 
 import numpy as np
-from distutils.version import LooseVersion
 from scipy.ndimage.filters import gaussian_filter
 
 import matplotlib.pyplot as plt
@@ -28,9 +27,10 @@
 
 from sklearn.feature_extraction.image import grid_to_graph
 from sklearn.cluster import AgglomerativeClustering
+from sklearn.utils.fixes import parse_version
 
 # these were introduced in skimage-0.14
-if LooseVersion(skimage.__version__) >= '0.14':
+if parse_version(skimage.__version__) >= parse_version('0.14'):
     rescale_params = {'anti_aliasing': False, 'multichannel': False}
 else:
     rescale_params = {}
diff --git a/examples/compose/plot_transformed_target.py b/examples/compose/plot_transformed_target.py
index 120948e9d685c..0d33f8009ad42 100755
--- a/examples/compose/plot_transformed_target.py
+++ b/examples/compose/plot_transformed_target.py
@@ -21,7 +21,6 @@
 import numpy as np
 import matplotlib
 import matplotlib.pyplot as plt
-from distutils.version import LooseVersion
 
 print(__doc__)
 
@@ -34,10 +33,11 @@
 from sklearn.linear_model import RidgeCV
 from sklearn.compose import TransformedTargetRegressor
 from sklearn.metrics import median_absolute_error, r2_score
+from sklearn.utils.fixes import parse_version
 
 
 # `normed` is being deprecated in favor of `density` in histograms
-if LooseVersion(matplotlib.__version__) >= '2.1':
+if parse_version(matplotlib.__version__) >= parse_version('2.1'):
     density_param = {'density': True}
 else:
     density_param = {'normed': True}
diff --git a/examples/decomposition/plot_sparse_coding.py b/examples/decomposition/plot_sparse_coding.py
index 8a87727a7c34b..4f8cf4035d6b0 100644
--- a/examples/decomposition/plot_sparse_coding.py
+++ b/examples/decomposition/plot_sparse_coding.py
@@ -16,12 +16,11 @@
 """
 print(__doc__)
 
-from distutils.version import LooseVersion
-
 import numpy as np
 import matplotlib.pyplot as plt
 
 from sklearn.decomposition import SparseCoder
+from sklearn.utils.fixes import np_version, parse_version
 
 
 def ricker_function(resolution, center, width):
@@ -68,7 +67,7 @@ def ricker_matrix(width, resolution, n_components):
               ('Lasso', 'lasso_lars', 2, None, 'turquoise'), ]
 lw = 2
 # Avoid FutureWarning about default value change when numpy >= 1.14
-lstsq_rcond = None if LooseVersion(np.__version__) >= '1.14' else -1
+lstsq_rcond = None if np_version >= parse_version('1.14') else -1
 
 plt.figure(figsize=(13, 6))
 for subplot, (D, title) in enumerate(zip((D_fixed, D_multi),
diff --git a/examples/miscellaneous/plot_johnson_lindenstrauss_bound.py b/examples/miscellaneous/plot_johnson_lindenstrauss_bound.py
index 9d369c6c6d46d..6e33aeea918f5 100644
--- a/examples/miscellaneous/plot_johnson_lindenstrauss_bound.py
+++ b/examples/miscellaneous/plot_johnson_lindenstrauss_bound.py
@@ -19,15 +19,15 @@
 import numpy as np
 import matplotlib
 import matplotlib.pyplot as plt
-from distutils.version import LooseVersion
 from sklearn.random_projection import johnson_lindenstrauss_min_dim
 from sklearn.random_projection import SparseRandomProjection
 from sklearn.datasets import fetch_20newsgroups_vectorized
 from sklearn.datasets import load_digits
 from sklearn.metrics.pairwise import euclidean_distances
+from sklearn.utils.fixes import parse_version
 
 # `normed` is being deprecated in favor of `density` in histograms
-if LooseVersion(matplotlib.__version__) >= '2.1':
+if parse_version(matplotlib.__version__) >= parse_version('2.1'):
     density_param = {'density': True}
 else:
     density_param = {'normed': True}
diff --git a/examples/neighbors/plot_kde_1d.py b/examples/neighbors/plot_kde_1d.py
index 39b1c38c1f59d..262391e926666 100644
--- a/examples/neighbors/plot_kde_1d.py
+++ b/examples/neighbors/plot_kde_1d.py
@@ -31,12 +31,12 @@
 import numpy as np
 import matplotlib
 import matplotlib.pyplot as plt
-from distutils.version import LooseVersion
 from scipy.stats import norm
 from sklearn.neighbors import KernelDensity
+from sklearn.utils.fixes import parse_version
 
 # `normed` is being deprecated in favor of `density` in histograms
-if LooseVersion(matplotlib.__version__) >= '2.1':
+if parse_version(matplotlib.__version__) >= parse_version('2.1'):
     density_param = {'density': True}
 else:
     density_param = {'normed': True}
diff --git a/sklearn/__init__.py b/sklearn/__init__.py
index 912752a8c676a..2f39a9fcc95ac 100644
--- a/sklearn/__init__.py
+++ b/sklearn/__init__.py
@@ -37,7 +37,7 @@
 # Dev branch marker is: 'X.Y.dev' or 'X.Y.devN' where N is an integer.
 # 'X.Y.dev0' is the canonical version of 'X.Y.dev'
 #
-__version__ = '0.23.1'
+__version__ = '0.23.2'
 
 
 # On OSX, we can get a runtime error due to multiple OpenMP libraries loaded
diff --git a/sklearn/_config.py b/sklearn/_config.py
index f183203e13228..e333095e817d1 100644
--- a/sklearn/_config.py
+++ b/sklearn/_config.py
@@ -61,8 +61,8 @@ def set_config(assume_finite=None, working_memory=None,
         .. versionadded:: 0.21
 
     display : {'text', 'diagram'}, optional
-        If 'diagram', estimators will be displayed as text in a jupyter lab
-        of notebook context. If 'text', estimators will be displayed as
+        If 'diagram', estimators will be displayed as a diagram in a Jupyter
+        lab or notebook context. If 'text', estimators will be displayed as
         text. Default is 'text'.
 
         .. versionadded:: 0.23
@@ -111,8 +111,8 @@ def config_context(**new_config):
            Default changed from False to True.
 
     display : {'text', 'diagram'}, optional
-        If 'diagram', estimators will be displayed as text in a jupyter lab
-        of notebook context. If 'text', estimators will be displayed as
+        If 'diagram', estimators will be displayed as a diagram in a Jupyter
+        lab or notebook context. If 'text', estimators will be displayed as
         text. Default is 'text'.
 
         .. versionadded:: 0.23
diff --git a/sklearn/cluster/_affinity_propagation.py b/sklearn/cluster/_affinity_propagation.py
index 3f4447eb52c5a..f26efb8ddef09 100644
--- a/sklearn/cluster/_affinity_propagation.py
+++ b/sklearn/cluster/_affinity_propagation.py
@@ -162,7 +162,7 @@ def affinity_propagation(S, *, preference=None, convergence_iter=15,
     tmp = np.zeros((n_samples, n_samples))
 
     # Remove degeneracies
-    S += ((np.finfo(np.double).eps * S + np.finfo(np.double).tiny * 100) *
+    S += ((np.finfo(S.dtype).eps * S + np.finfo(S.dtype).tiny * 100) *
           random_state.randn(n_samples, n_samples))
 
     # Execute parallel affinity propagation updates
diff --git a/sklearn/cluster/_kmeans.py b/sklearn/cluster/_kmeans.py
index 38d8120bf452c..421d7e94d23e5 100644
--- a/sklearn/cluster/_kmeans.py
+++ b/sklearn/cluster/_kmeans.py
@@ -22,7 +22,7 @@
 from ..utils.extmath import row_norms, stable_cumsum
 from ..utils.sparsefuncs_fast import assign_rows_csr
 from ..utils.sparsefuncs import mean_variance_axis
-from ..utils.validation import _num_samples, _deprecate_positional_args
+from ..utils.validation import _deprecate_positional_args
 from ..utils import check_array
 from ..utils import gen_batches
 from ..utils import check_random_state
@@ -145,15 +145,14 @@ def _k_init(X, n_clusters, x_squared_norms, random_state, n_local_trials=None):
 
 def _validate_center_shape(X, n_centers, centers):
     """Check if centers is compatible with X and n_centers"""
-    if len(centers) != n_centers:
-        raise ValueError('The shape of the initial centers (%s) '
-                         'does not match the number of clusters %i'
-                         % (centers.shape, n_centers))
+    if centers.shape[0] != n_centers:
+        raise ValueError(
+            f"The shape of the initial centers {centers.shape} does not "
+            f"match the number of clusters {n_centers}.")
     if centers.shape[1] != X.shape[1]:
         raise ValueError(
-            "The number of features of the initial centers %s "
-            "does not match the number of features of the data %s."
-            % (centers.shape[1], X.shape[1]))
+            f"The shape of the initial centers {centers.shape} does not "
+            f"match the number of features of the data {X.shape[1]}.")
 
 
 def _tolerance(X, tol):
@@ -167,21 +166,6 @@ def _tolerance(X, tol):
     return np.mean(variances) * tol
 
 
-def _check_normalize_sample_weight(sample_weight, X):
-    """Set sample_weight if None, and check for correct dtype"""
-
-    sample_weight_was_none = sample_weight is None
-
-    sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
-    if not sample_weight_was_none:
-        # normalize the weights to sum up to n_samples
-        # an array of 1 (i.e. samples_weight is None) is already normalized
-        n_samples = len(sample_weight)
-        scale = n_samples / sample_weight.sum()
-        sample_weight = sample_weight * scale
-    return sample_weight
-
-
 @_deprecate_positional_args
 def k_means(X, n_clusters, *, sample_weight=None, init='k-means++',
             precompute_distances='deprecated', n_init=10, max_iter=300,
@@ -252,8 +236,6 @@ def k_means(X, n_clusters, *, sample_weight=None, init='k-means++',
         Relative tolerance with regards to Frobenius norm of the difference
         in the cluster centers of two consecutive iterations to declare
         convergence.
-        It's not advised to set `tol=0` since convergence might never be
-        declared due to rounding errors. Use a very small number instead.
 
     random_state : int, RandomState instance, default=None
         Determines random number generation for centroid initialization. Use
@@ -399,7 +381,7 @@ def _kmeans_single_elkan(X, sample_weight, n_clusters, max_iter=300,
         Number of iterations run.
     """
     random_state = check_random_state(random_state)
-    sample_weight = _check_normalize_sample_weight(sample_weight, X)
+    sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
 
     # init
     centers = _init_centroids(X, n_clusters, init, random_state=random_state,
@@ -413,6 +395,7 @@ def _kmeans_single_elkan(X, sample_weight, n_clusters, max_iter=300,
     centers_new = np.zeros_like(centers)
     weight_in_clusters = np.zeros(n_clusters, dtype=X.dtype)
     labels = np.full(n_samples, -1, dtype=np.int32)
+    labels_old = labels.copy()
     center_half_distances = euclidean_distances(centers) / 2
     distance_next_center = np.partition(np.asarray(center_half_distances),
                                         kth=1, axis=0)[1]
@@ -432,6 +415,8 @@ def _kmeans_single_elkan(X, sample_weight, n_clusters, max_iter=300,
     init_bounds(X, centers, center_half_distances,
                 labels, upper_bounds, lower_bounds)
 
+    strict_convergence = False
+
     for i in range(max_iter):
         elkan_iter(X, sample_weight, centers, centers_new,
                    weight_in_clusters, center_half_distances,
@@ -448,17 +433,25 @@ def _kmeans_single_elkan(X, sample_weight, n_clusters, max_iter=300,
             inertia = _inertia(X, sample_weight, centers, labels)
             print("Iteration {0}, inertia {1}" .format(i, inertia))
 
-        center_shift_tot = (center_shift**2).sum()
-        if center_shift_tot <= tol:
+        if np.array_equal(labels, labels_old):
+            # First check the labels for strict convergence.
             if verbose:
-                print("Converged at iteration {0}: "
-                      "center shift {1} within tolerance {2}"
-                      .format(i, center_shift_tot, tol))
+                print(f"Converged at iteration {i}: strict convergence.")
+            strict_convergence = True
             break
+        else:
+            # No strict convergence, check for tol based convergence.
+            center_shift_tot = (center_shift**2).sum()
+            if center_shift_tot <= tol:
+                if verbose:
+                    print(f"Converged at iteration {i}: center shift "
+                          f"{center_shift_tot} within tolerance {tol}.")
+                break
 
         centers, centers_new = centers_new, centers
+        labels_old[:] = labels
 
-    if center_shift_tot > 0:
+    if not strict_convergence:
         # rerun E-step so that predicted labels match cluster centers
         elkan_iter(X, sample_weight, centers, centers, weight_in_clusters,
                    center_half_distances, distance_next_center,
@@ -546,7 +539,7 @@ def _kmeans_single_lloyd(X, sample_weight, n_clusters, max_iter=300,
         Number of iterations run.
     """
     random_state = check_random_state(random_state)
-    sample_weight = _check_normalize_sample_weight(sample_weight, X)
+    sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
 
     # init
     centers = _init_centroids(X, n_clusters, init, random_state=random_state,
@@ -557,6 +550,7 @@ def _kmeans_single_lloyd(X, sample_weight, n_clusters, max_iter=300,
 
     centers_new = np.zeros_like(centers)
     labels = np.full(X.shape[0], -1, dtype=np.int32)
+    labels_old = labels.copy()
     weight_in_clusters = np.zeros(n_clusters, dtype=X.dtype)
     center_shift = np.zeros(n_clusters, dtype=X.dtype)
 
@@ -567,6 +561,8 @@ def _kmeans_single_lloyd(X, sample_weight, n_clusters, max_iter=300,
         lloyd_iter = lloyd_iter_chunked_dense
         _inertia = _inertia_dense
 
+    strict_convergence = False
+
     # Threadpoolctl context to limit the number of threads in second level of
     # nested parallelism (i.e. BLAS) to avoid oversubsciption.
     with threadpool_limits(limits=1, user_api="blas"):
@@ -578,17 +574,25 @@ def _kmeans_single_lloyd(X, sample_weight, n_clusters, max_iter=300,
                 inertia = _inertia(X, sample_weight, centers, labels)
                 print("Iteration {0}, inertia {1}" .format(i, inertia))
 
-            center_shift_tot = (center_shift**2).sum()
-            if center_shift_tot <= tol:
+            if np.array_equal(labels, labels_old):
+                # First check the labels for strict convergence.
                 if verbose:
-                    print("Converged at iteration {0}: "
-                          "center shift {1} within tolerance {2}"
-                          .format(i, center_shift_tot, tol))
+                    print(f"Converged at iteration {i}: strict convergence.")
+                strict_convergence = True
                 break
+            else:
+                # No strict convergence, check for tol based convergence.
+                center_shift_tot = (center_shift**2).sum()
+                if center_shift_tot <= tol:
+                    if verbose:
+                        print(f"Converged at iteration {i}: center shift "
+                              f"{center_shift_tot} within tolerance {tol}.")
+                    break
 
             centers, centers_new = centers_new, centers
+            labels_old[:] = labels
 
-        if center_shift_tot > 0:
+        if not strict_convergence:
             # rerun E-step so that predicted labels match cluster centers
             lloyd_iter(X, sample_weight, x_squared_norms, centers, centers,
                        weight_in_clusters, labels, center_shift, n_threads,
@@ -639,7 +643,7 @@ def _labels_inertia(X, sample_weight, x_squared_norms, centers,
 
     n_threads = _openmp_effective_n_threads(n_threads)
 
-    sample_weight = _check_normalize_sample_weight(sample_weight, X)
+    sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
     labels = np.full(n_samples, -1, dtype=np.int32)
     weight_in_clusters = np.zeros(n_clusters, dtype=centers.dtype)
     center_shift = np.zeros_like(weight_in_clusters)
@@ -783,8 +787,6 @@ class KMeans(TransformerMixin, ClusterMixin, BaseEstimator):
         Relative tolerance with regards to Frobenius norm of the difference
         in the cluster centers of two consecutive iterations to declare
         convergence.
-        It's not advised to set `tol=0` since convergence might never be
-        declared due to rounding errors. Use a very small number instead.
 
     precompute_distances : {'auto', True, False}, default='auto'
         Precompute distances (faster but takes more memory).
@@ -924,15 +926,78 @@ def __init__(self, n_clusters=8, *, init='k-means++', n_init=10,
         self.n_jobs = n_jobs
         self.algorithm = algorithm
 
+    def _check_params(self, X):
+        # precompute_distances
+        if self.precompute_distances != 'deprecated':
+            warnings.warn("'precompute_distances' was deprecated in version "
+                          "0.23 and will be removed in 0.25. It has no "
+                          "effect", FutureWarning)
+
+        # n_jobs
+        if self.n_jobs != 'deprecated':
+            warnings.warn("'n_jobs' was deprecated in version 0.23 and will be"
+                          " removed in 0.25.", FutureWarning)
+            self._n_threads = self.n_jobs
+        else:
+            self._n_threads = None
+        self._n_threads = _openmp_effective_n_threads(self._n_threads)
+
+        # n_init
+        if self.n_init <= 0:
+            raise ValueError(
+                f"n_init should be > 0, got {self.n_init} instead.")
+        self._n_init = self.n_init
+
+        # max_iter
+        if self.max_iter <= 0:
+            raise ValueError(
+                f"max_iter should be > 0, got {self.max_iter} instead.")
+
+        # n_clusters
+        if X.shape[0] < self.n_clusters:
+            raise ValueError(f"n_samples={X.shape[0]} should be >= "
+                             f"n_clusters={self.n_clusters}.")
+
+        # tol
+        self._tol = _tolerance(X, self.tol)
+
+        # algorithm
+        if self.algorithm not in ("auto", "full", "elkan"):
+            raise ValueError(f"Algorithm must be 'auto', 'full' or 'elkan', "
+                             f"got {self.algorithm} instead.")
+
+        self._algorithm = self.algorithm
+        if self._algorithm == "auto":
+            self._algorithm = "full" if self.n_clusters == 1 else "elkan"
+        if self._algorithm == "elkan" and self.n_clusters == 1:
+            warnings.warn("algorithm='elkan' doesn't make sense for a single "
+                          "cluster. Using 'full' instead.", RuntimeWarning)
+            self._algorithm = "full"
+
+        # init
+        if not (hasattr(self.init, '__array__') or callable(self.init)
+                or (isinstance(self.init, str)
+                    and self.init in ["k-means++", "random"])):
+            raise ValueError(
+                f"init should be either 'k-means++', 'random', a ndarray or a "
+                f"callable, got '{self.init}' instead.")
+
+        if hasattr(self.init, '__array__') and self._n_init != 1:
+            warnings.warn(
+                f"Explicit initial center position passed: performing only"
+                f" one init in {self.__class__.__name__} instead of "
+                f"n_init={self._n_init}.", RuntimeWarning, stacklevel=2)
+            self._n_init = 1
+
     def _check_test_data(self, X):
         X = check_array(X, accept_sparse='csr', dtype=[np.float64, np.float32],
                         order='C', accept_large_sparse=False)
         n_samples, n_features = X.shape
         expected_n_features = self.cluster_centers_.shape[1]
         if not n_features == expected_n_features:
-            raise ValueError("Incorrect number of features. "
-                             "Got %d features, expected %d" % (
-                                 n_features, expected_n_features))
+            raise ValueError(
+                f"Incorrect number of features. Got {n_features} features, "
+                f"expected {expected_n_features}.")
 
         return X
 
@@ -962,56 +1027,20 @@ def fit(self, X, y=None, sample_weight=None):
         self
             Fitted estimator.
         """
-        random_state = check_random_state(self.random_state)
-
-        if self.precompute_distances != 'deprecated':
-            warnings.warn("'precompute_distances' was deprecated in version "
-                          "0.23 and will be removed in 0.25. It has no "
-                          "effect", FutureWarning)
-
-        if self.n_jobs != 'deprecated':
-            warnings.warn("'n_jobs' was deprecated in version 0.23 and will be"
-                          " removed in 0.25.", FutureWarning)
-            self._n_threads = self.n_jobs
-        else:
-            self._n_threads = None
-        self._n_threads = _openmp_effective_n_threads(self._n_threads)
-
-        n_init = self.n_init
-        if n_init <= 0:
-            raise ValueError("Invalid number of initializations."
-                             " n_init=%d must be bigger than zero." % n_init)
-
-        if self.max_iter <= 0:
-            raise ValueError(
-                'Number of iterations should be a positive number,'
-                ' got %d instead' % self.max_iter
-            )
-
         X = self._validate_data(X, accept_sparse='csr',
                                 dtype=[np.float64, np.float32],
                                 order='C', copy=self.copy_x,
                                 accept_large_sparse=False)
-        # verify that the number of samples given is larger than k
-        if _num_samples(X) < self.n_clusters:
-            raise ValueError("n_samples=%d should be >= n_clusters=%d" % (
-                _num_samples(X), self.n_clusters))
 
-        tol = _tolerance(X, self.tol)
+        self._check_params(X)
+        random_state = check_random_state(self.random_state)
 
         # Validate init array
         init = self.init
         if hasattr(init, '__array__'):
-            init = check_array(init, dtype=X.dtype.type, copy=True, order='C')
+            init = check_array(init, dtype=X.dtype, copy=True, order='C')
             _validate_center_shape(X, self.n_clusters, init)
 
-            if n_init != 1:
-                warnings.warn(
-                    'Explicit initial center position passed: '
-                    'performing only one init in k-means instead of n_init=%d'
-                    % n_init, RuntimeWarning, stacklevel=2)
-                n_init = 1
-
         # subtract of mean of x for more accurate distance computations
         if not sp.issparse(X):
             X_mean = X.mean(axis=0)
@@ -1024,33 +1053,21 @@ def fit(self, X, y=None, sample_weight=None):
         # precompute squared norms of data points
         x_squared_norms = row_norms(X, squared=True)
 
-        best_labels, best_inertia, best_centers = None, None, None
-
-        algorithm = self.algorithm
-        if algorithm == "elkan" and self.n_clusters == 1:
-            warnings.warn("algorithm='elkan' doesn't make sense for a single "
-                          "cluster. Using 'full' instead.", RuntimeWarning)
-            algorithm = "full"
-
-        if algorithm == "auto":
-            algorithm = "full" if self.n_clusters == 1 else "elkan"
-
-        if algorithm == "full":
+        if self._algorithm == "full":
             kmeans_single = _kmeans_single_lloyd
-        elif algorithm == "elkan":
-            kmeans_single = _kmeans_single_elkan
         else:
-            raise ValueError("Algorithm must be 'auto', 'full' or 'elkan', got"
-                             " {}".format(str(algorithm)))
+            kmeans_single = _kmeans_single_elkan
+
+        best_labels, best_inertia, best_centers = None, None, None
 
         # seeds for the initializations of the kmeans runs.
-        seeds = random_state.randint(np.iinfo(np.int32).max, size=n_init)
+        seeds = random_state.randint(np.iinfo(np.int32).max, size=self._n_init)
 
         for seed in seeds:
             # run a k-means once
             labels, inertia, centers, n_iter_ = kmeans_single(
                 X, sample_weight, self.n_clusters, max_iter=self.max_iter,
-                init=init, verbose=self.verbose, tol=tol,
+                init=init, verbose=self.verbose, tol=self._tol,
                 x_squared_norms=x_squared_norms, random_state=seed,
                 n_threads=self._n_threads)
             # determine if these results are the best so far
@@ -1581,6 +1598,46 @@ def __init__(self, n_clusters=8, *, init='k-means++', max_iter=100,
         self.init_size = init_size
         self.reassignment_ratio = reassignment_ratio
 
+    def _check_params(self, X):
+        super()._check_params(X)
+
+        # max_no_improvement
+        if self.max_no_improvement is not None and self.max_no_improvement < 0:
+            raise ValueError(
+                f"max_no_improvement should be >= 0, got "
+                f"{self.max_no_improvement} instead.")
+
+        # batch_size
+        if self.batch_size <= 0:
+            raise ValueError(
+                f"batch_size should be > 0, got {self.batch_size} instead.")
+
+        # init_size
+        if self.init_size is not None and self.init_size <= 0:
+            raise ValueError(
+                f"init_size should be > 0, got {self.init_size} instead.")
+        self._init_size = self.init_size
+        if self._init_size is None:
+            self._init_size = 3 * self.batch_size
+            if self._init_size < self.n_clusters:
+                self._init_size = 3 * self.n_clusters
+        elif self._init_size < self.n_clusters:
+            warnings.warn(
+                f"init_size={self._init_size} should be larger than "
+                f"n_clusters={self.n_clusters}. Setting it to "
+                f"min(3*n_clusters, n_samples)",
+                RuntimeWarning, stacklevel=2)
+            self._init_size = 3 * self.n_clusters
+        self._init_size = min(self._init_size, X.shape[0])
+        # FIXME: init_size_ will be deprecated and this line will be removed
+        self.init_size_ = self._init_size
+
+        # reassignment_ratio
+        if self.reassignment_ratio < 0:
+            raise ValueError(
+                f"reassignment_ratio should be >= 0, got "
+                f"{self.reassignment_ratio} instead.")
+
     def fit(self, X, y=None, sample_weight=None):
         """Compute the centroids on X by chunking it into mini-batches.
 
@@ -1604,27 +1661,21 @@ def fit(self, X, y=None, sample_weight=None):
         -------
         self
         """
+        X = self._validate_data(X, accept_sparse='csr',
+                                dtype=[np.float64, np.float32],
+                                order='C', accept_large_sparse=False)
+
+        self._check_params(X)
         random_state = check_random_state(self.random_state)
-        X = self._validate_data(X, accept_sparse="csr", order='C',
-                                dtype=[np.float64, np.float32])
-        n_samples, n_features = X.shape
-        if n_samples < self.n_clusters:
-            raise ValueError("n_samples=%d should be >= n_clusters=%d"
-                             % (n_samples, self.n_clusters))
-
-        sample_weight = _check_normalize_sample_weight(sample_weight, X)
-
-        n_init = self.n_init
-        if hasattr(self.init, '__array__'):
-            self.init = np.ascontiguousarray(self.init, dtype=X.dtype)
-            if n_init != 1:
-                warnings.warn(
-                    'Explicit initial center position passed: '
-                    'performing only one init in MiniBatchKMeans instead of '
-                    'n_init=%d'
-                    % self.n_init, RuntimeWarning, stacklevel=2)
-                n_init = 1
+        sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
 
+        # Validate init array
+        init = self.init
+        if hasattr(init, '__array__'):
+            init = check_array(init, dtype=X.dtype, copy=True, order='C')
+            _validate_center_shape(X, self.n_clusters, init)
+
+        n_samples, n_features = X.shape
         x_squared_norms = row_norms(X, squared=True)
 
         if self.tol > 0.0:
@@ -1644,24 +1695,18 @@ def fit(self, X, y=None, sample_weight=None):
         n_batches = int(np.ceil(float(n_samples) / self.batch_size))
         n_iter = int(self.max_iter * n_batches)
 
-        init_size = self.init_size
-        if init_size is None:
-            init_size = 3 * self.batch_size
-        if init_size > n_samples:
-            init_size = n_samples
-        self.init_size_ = init_size
-
-        validation_indices = random_state.randint(0, n_samples, init_size)
+        validation_indices = random_state.randint(0, n_samples,
+                                                  self._init_size)
         X_valid = X[validation_indices]
         sample_weight_valid = sample_weight[validation_indices]
         x_squared_norms_valid = x_squared_norms[validation_indices]
 
         # perform several inits with random sub-sets
         best_inertia = None
-        for init_idx in range(n_init):
+        for init_idx in range(self._n_init):
             if self.verbose:
                 print("Init %d/%d with method: %s"
-                      % (init_idx + 1, n_init, self.init))
+                      % (init_idx + 1, self._n_init, init))
             weight_sums = np.zeros(self.n_clusters, dtype=sample_weight.dtype)
 
             # TODO: once the `k_means` function works with sparse input we
@@ -1670,10 +1715,10 @@ def fit(self, X, y=None, sample_weight=None):
             # Initialize the centers using only a fraction of the data as we
             # expect n_samples to be very large when using MiniBatchKMeans
             cluster_centers = _init_centroids(
-                X, self.n_clusters, self.init,
+                X, self.n_clusters, init,
                 random_state=random_state,
                 x_squared_norms=x_squared_norms,
-                init_size=init_size)
+                init_size=self._init_size)
 
             # Compute the label assignment on the init dataset
             _mini_batch_step(
@@ -1689,7 +1734,7 @@ def fit(self, X, y=None, sample_weight=None):
                                          cluster_centers)
             if self.verbose:
                 print("Inertia for init %d/%d: %f"
-                      % (init_idx + 1, n_init, inertia))
+                      % (init_idx + 1, self._n_init, inertia))
             if best_inertia is None or inertia < best_inertia:
                 self.cluster_centers_ = cluster_centers
                 self.counts_ = weight_sums
@@ -1761,7 +1806,7 @@ def _labels_inertia_minibatch(self, X, sample_weight):
         """
         if self.verbose:
             print('Computing label assignment and total inertia')
-        sample_weight = _check_normalize_sample_weight(sample_weight, X)
+        sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
         x_squared_norms = row_norms(X, squared=True)
         slices = gen_batches(X.shape[0], self.batch_size)
         results = [_labels_inertia(X[s], sample_weight[s], x_squared_norms[s],
@@ -1789,27 +1834,32 @@ def partial_fit(self, X, y=None, sample_weight=None):
         -------
         self
         """
+        is_first_call_to_partial_fit = not hasattr(self, 'cluster_centers_')
 
-        X = check_array(X, accept_sparse="csr", order="C",
-                        dtype=[np.float64, np.float32])
-        n_samples, n_features = X.shape
-        if hasattr(self.init, '__array__'):
-            self.init = np.ascontiguousarray(self.init, dtype=X.dtype)
-
-        if n_samples == 0:
-            return self
-
-        sample_weight = _check_normalize_sample_weight(sample_weight, X)
+        X = self._validate_data(X, accept_sparse='csr',
+                                dtype=[np.float64, np.float32],
+                                order='C', accept_large_sparse=False,
+                                reset=is_first_call_to_partial_fit)
 
-        x_squared_norms = row_norms(X, squared=True)
         self.random_state_ = getattr(self, "random_state_",
                                      check_random_state(self.random_state))
-        if (not hasattr(self, 'counts_')
-                or not hasattr(self, 'cluster_centers_')):
-            # this is the first call partial_fit on this object:
+        sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
+
+        x_squared_norms = row_norms(X, squared=True)
+
+        if is_first_call_to_partial_fit:
+            # this is the first call to partial_fit on this object
+            self._check_params(X)
+
+            # Validate init array
+            init = self.init
+            if hasattr(init, '__array__'):
+                init = check_array(init, dtype=X.dtype, copy=True, order='C')
+                _validate_center_shape(X, self.n_clusters, init)
+
             # initialize the cluster centers
             self.cluster_centers_ = _init_centroids(
-                X, self.n_clusters, self.init,
+                X, self.n_clusters, init,
                 random_state=self.random_state_,
                 x_squared_norms=x_squared_norms, init_size=self.init_size)
 
@@ -1825,13 +1875,6 @@ def partial_fit(self, X, y=None, sample_weight=None):
                 10 * (1 + self.counts_.min())) == 0
             distances = np.zeros(X.shape[0], dtype=X.dtype)
 
-            # Raise error if partial_fit called on data with different number
-            # of features.
-            if X.shape[1] != self.cluster_centers_.shape[1]:
-                raise ValueError(
-                    "Number of features %d does not match previous "
-                    "data %d." % (X.shape[1], self.cluster_centers_.shape[1]))
-
         _mini_batch_step(X, sample_weight, x_squared_norms,
                          self.cluster_centers_, self.counts_,
                          np.zeros(0, dtype=X.dtype), 0,
diff --git a/sklearn/cluster/_mean_shift.py b/sklearn/cluster/_mean_shift.py
index 3c0c34cdbdffa..c56853770240b 100644
--- a/sklearn/cluster/_mean_shift.py
+++ b/sklearn/cluster/_mean_shift.py
@@ -218,6 +218,8 @@ def get_bin_seeds(X, bin_size, min_bin_freq=1):
     bin_seeds : array-like of shape (n_samples, n_features)
         Points used as initial kernel positions in clustering.mean_shift.
     """
+    if bin_size == 0:
+        return X
 
     # Bin points
     bin_sizes = defaultdict(int)
diff --git a/sklearn/cluster/tests/test_affinity_propagation.py b/sklearn/cluster/tests/test_affinity_propagation.py
index 826878061291b..1369d566c16fe 100644
--- a/sklearn/cluster/tests/test_affinity_propagation.py
+++ b/sklearn/cluster/tests/test_affinity_propagation.py
@@ -231,3 +231,16 @@ def test_affinity_propagation_convergence_warning_dense_sparse(centers):
         assert_array_equal(ap.predict(X),
                            np.zeros(X.shape[0], dtype=int))
     assert len(record) == 0
+
+
+def test_affinity_propagation_float32():
+    # Test to fix incorrect clusters due to dtype change
+    # (non-regression test for issue #10832)
+    X = np.array([[1, 0, 0, 0],
+                  [0, 1, 1, 0],
+                  [0, 1, 1, 0],
+                  [0, 0, 0, 1]], dtype='float32')
+    afp = AffinityPropagation(preference=1, affinity='precomputed',
+                              random_state=0).fit(X)
+    expected = np.array([0, 1, 1, 2])
+    assert_array_equal(afp.labels_, expected)
diff --git a/sklearn/cluster/tests/test_k_means.py b/sklearn/cluster/tests/test_k_means.py
index c821e7e13d58b..b53fd5be22c7c 100644
--- a/sklearn/cluster/tests/test_k_means.py
+++ b/sklearn/cluster/tests/test_k_means.py
@@ -1,4 +1,5 @@
 """Testing for K-means"""
+import re
 import sys
 
 import numpy as np
@@ -14,7 +15,7 @@
 from sklearn.utils._testing import assert_warns
 from sklearn.utils._testing import assert_warns_message
 from sklearn.utils._testing import assert_raise_message
-from sklearn.utils.validation import _num_samples
+from sklearn.utils.fixes import _astype_copy_false
 from sklearn.base import clone
 from sklearn.exceptions import ConvergenceWarning
 
@@ -49,18 +50,19 @@
 X_csr = sp.csr_matrix(X)
 
 
-@pytest.mark.parametrize("representation", ["dense", "sparse"])
+@pytest.mark.parametrize("array_constr", [np.array, sp.csr_matrix],
+                         ids=["dense", "sparse"])
 @pytest.mark.parametrize("algo", ["full", "elkan"])
 @pytest.mark.parametrize("dtype", [np.float32, np.float64])
-def test_kmeans_results(representation, algo, dtype):
-    # cheks that kmeans works as intended
-    array_constr = {'dense': np.array, 'sparse': sp.csr_matrix}[representation]
+def test_kmeans_results(array_constr, algo, dtype):
+    # Checks that KMeans works as intended on toy dataset by comparing with
+    # expected results computed by hand.
     X = array_constr([[0, 0], [0.5, 0], [0.5, 1], [1, 1]], dtype=dtype)
-    sample_weight = [3, 1, 1, 3]  # will be rescaled to [1.5, 0.5, 0.5, 1.5]
+    sample_weight = [3, 1, 1, 3]
     init_centers = np.array([[0, 0], [1, 1]], dtype=dtype)
 
     expected_labels = [0, 0, 1, 1]
-    expected_inertia = 0.1875
+    expected_inertia = 0.375
     expected_centers = np.array([[0.125, 0], [0.875, 1]], dtype=dtype)
     expected_n_iter = 2
 
@@ -68,8 +70,8 @@ def test_kmeans_results(representation, algo, dtype):
     kmeans.fit(X, sample_weight=sample_weight)
 
     assert_array_equal(kmeans.labels_, expected_labels)
-    assert_almost_equal(kmeans.inertia_, expected_inertia)
-    assert_array_almost_equal(kmeans.cluster_centers_, expected_centers)
+    assert_allclose(kmeans.inertia_, expected_inertia)
+    assert_allclose(kmeans.cluster_centers_, expected_centers)
     assert kmeans.n_iter_ == expected_n_iter
 
 
@@ -135,10 +137,12 @@ def test_relocate_empty_clusters(representation):
     assert_allclose(centers_new, [[-36], [10], [9.5]])
 
 
-@pytest.mark.parametrize('distribution', ['normal', 'blobs'])
-@pytest.mark.parametrize('tol', [1e-2, 1e-4, 1e-8])
-def test_elkan_results(distribution, tol):
-    # check that results are identical between lloyd and elkan algorithms
+@pytest.mark.parametrize("distribution", ["normal", "blobs"])
+@pytest.mark.parametrize("array_constr", [np.array, sp.csr_matrix],
+                         ids=["dense", "sparse"])
+@pytest.mark.parametrize("tol", [1e-2, 1e-8, 1e-100, 0])
+def test_kmeans_elkan_results(distribution, array_constr, tol):
+    # Check that results are identical between lloyd and elkan algorithms
     rnd = np.random.RandomState(0)
     if distribution == 'normal':
         X = rnd.normal(size=(5000, 10))
@@ -164,11 +168,12 @@ def test_kmeans_convergence(algorithm):
     # Check that KMeans stops when convergence is reached when tol=0. (#16075)
     rnd = np.random.RandomState(0)
     X = rnd.normal(size=(5000, 10))
+    max_iter = 300
 
-    km = KMeans(algorithm=algorithm, n_clusters=5, random_state=0, n_init=1,
-                tol=0, max_iter=300).fit(X)
+    km = KMeans(algorithm=algorithm, n_clusters=5, random_state=0,
+                n_init=1, tol=0, max_iter=max_iter).fit(X)
 
-    assert km.n_iter_ < 300
+    assert km.n_iter_ < max_iter
 
 
 @pytest.mark.parametrize('distribution', ['normal', 'blobs'])
@@ -343,47 +348,16 @@ def test_k_means_init(data, init):
     _check_fitted_model(km)
 
 
-def test_k_means_n_init():
-    rnd = np.random.RandomState(0)
-    X = rnd.normal(size=(40, 2))
-
-    # two regression tests on bad n_init argument
-    # previous bug: n_init <= 0 threw non-informative TypeError (#3858)
-    with pytest.raises(ValueError, match="n_init"):
-        KMeans(n_init=0).fit(X)
-    with pytest.raises(ValueError, match="n_init"):
-        KMeans(n_init=-1).fit(X)
-
-
-@pytest.mark.parametrize('Class', [KMeans, MiniBatchKMeans])
-def test_k_means_explicit_init_shape(Class):
-    # test for sensible errors when giving explicit init
-    # with wrong number of features or clusters
-    rnd = np.random.RandomState(0)
-    X = rnd.normal(size=(40, 3))
-
-    # mismatch of number of features
-    km = Class(n_init=1, init=X[:, :2], n_clusters=len(X))
-    msg = "does not match the number of features of the data"
-    with pytest.raises(ValueError, match=msg):
-        km.fit(X)
-    # for callable init
-    km = Class(n_init=1,
-               init=lambda X_, k, random_state: X_[:, :2],
-               n_clusters=len(X))
-    with pytest.raises(ValueError, match=msg):
-        km.fit(X)
-    # mismatch of number of clusters
-    msg = "does not match the number of clusters"
-    km = Class(n_init=1, init=X[:2, :], n_clusters=3)
-    with pytest.raises(ValueError, match=msg):
-        km.fit(X)
-    # for callable init
-    km = Class(n_init=1,
-               init=lambda X_, k, random_state: X_[:2, :],
-               n_clusters=3)
-    with pytest.raises(ValueError, match=msg):
-        km.fit(X)
+@pytest.mark.parametrize("init", ["random", "k-means++", centers,
+                                  lambda X, k, random_state: centers],
+                         ids=["random", "k-means++", "ndarray", "callable"])
+def test_minibatch_kmeans_partial_fit_init(init):
+    # Check MiniBatchKMeans init with partial_fit
+    km = MiniBatchKMeans(init=init, n_clusters=n_clusters, random_state=0)
+    for i in range(100):
+        # "random" init requires many batches to recover the true labels.
+        km.partial_fit(X)
+    _check_fitted_model(km)
 
 
 def test_k_means_fortran_aligned_data():
@@ -439,22 +413,42 @@ def test_k_means_fit_predict(algo, dtype, constructor, seed, max_iter, tol):
     assert v_measure_score(labels_1, labels_2) == 1
 
 
-def test_mb_kmeans_verbose():
-    mb_k_means = MiniBatchKMeans(init="k-means++", n_clusters=n_clusters,
-                                 random_state=42, verbose=1)
+def test_minibatch_kmeans_verbose():
+    # Check verbose mode of MiniBatchKMeans for better coverage.
+    km = MiniBatchKMeans(n_clusters=n_clusters, random_state=42, verbose=1)
     old_stdout = sys.stdout
     sys.stdout = StringIO()
     try:
-        mb_k_means.fit(X)
+        km.fit(X)
     finally:
         sys.stdout = old_stdout
 
 
-def test_minibatch_init_with_large_k():
-    mb_k_means = MiniBatchKMeans(init='k-means++', init_size=10, n_clusters=20)
-    # Check that a warning is raised, as the number clusters is larger
-    # than the init_size
-    assert_warns(RuntimeWarning, mb_k_means.fit, X)
+@pytest.mark.parametrize("algorithm", ["full", "elkan"])
+@pytest.mark.parametrize("tol", [1e-2, 0])
+def test_kmeans_verbose(algorithm, tol, capsys):
+    # Check verbose mode of KMeans for better coverage.
+    X = np.random.RandomState(0).normal(size=(5000, 10))
+
+    KMeans(algorithm=algorithm, n_clusters=n_clusters, random_state=42,
+           init="random", n_init=1, tol=tol, verbose=1).fit(X)
+
+    captured = capsys.readouterr()
+
+    assert re.search(r"Initialization complete", captured.out)
+    assert re.search(r"Iteration [0-9]+, inertia", captured.out)
+
+    if tol == 0:
+        assert re.search(r"strict convergence", captured.out)
+    else:
+        assert re.search(r"center shift .* within tolerance", captured.out)
+
+
+def test_minibatch_kmeans_warning_init_size():
+    # Check that a warning is raised when init_size is smaller than n_clusters
+    with pytest.warns(RuntimeWarning,
+                      match=r"init_size.* should be larger than n_clusters"):
+        MiniBatchKMeans(init_size=10, n_clusters=20).fit(X)
 
 
 def test_minibatch_k_means_init_multiple_runs_with_explicit_centers():
@@ -571,13 +565,6 @@ def test_sparse_mb_k_means_callable_init():
     def test_init(X, k, random_state):
         return centers
 
-    # Small test to check that giving the wrong number of centers
-    # raises a meaningful error
-    msg = "does not match the number of clusters"
-    with pytest.raises(ValueError, match=msg):
-        MiniBatchKMeans(init=test_init, random_state=42).fit(X_csr)
-
-    # Now check that the fit actually works
     mb_k_means = MiniBatchKMeans(n_clusters=3, init=test_init,
                                  random_state=42).fit(X_csr)
     _check_fitted_model(mb_k_means)
@@ -595,12 +582,21 @@ def test_mini_batch_k_means_random_init_partial_fit():
     assert v_measure_score(true_labels, labels) == 1.0
 
 
-def test_minibatch_default_init_size():
-    mb_k_means = MiniBatchKMeans(init=centers.copy(), n_clusters=n_clusters,
-                                 batch_size=10, random_state=42,
-                                 n_init=1).fit(X)
-    assert mb_k_means.init_size_ == 3 * mb_k_means.batch_size
-    _check_fitted_model(mb_k_means)
+def test_minibatch_kmeans_default_init_size():
+    # Check the internal _init_size attribute of MiniBatchKMeans
+
+    # default init size should be 3 * batch_size
+    km = MiniBatchKMeans(n_clusters=10, batch_size=5, n_init=1).fit(X)
+    assert km._init_size == 15
+
+    # if 3 * batch size < n_clusters, it should then be 3 * n_clusters
+    km = MiniBatchKMeans(n_clusters=10, batch_size=1, n_init=1).fit(X)
+    assert km._init_size == 30
+
+    # it should not be larger than n_samples
+    km = MiniBatchKMeans(n_clusters=10, batch_size=5, n_init=1,
+                         init_size=n_samples + 1).fit(X)
+    assert km._init_size == n_samples
 
 
 def test_minibatch_tol():
@@ -618,13 +614,6 @@ def test_minibatch_set_init_size():
     _check_fitted_model(mb_k_means)
 
 
-@pytest.mark.parametrize("Estimator", [KMeans, MiniBatchKMeans])
-def test_k_means_invalid_init(Estimator):
-    km = Estimator(init="invalid", n_init=1, n_clusters=n_clusters)
-    with pytest.raises(ValueError):
-        km.fit(X)
-
-
 def test_k_means_copyx():
     # Check if copy_x=False returns nearly equal X after de-centering.
     my_X = X.copy()
@@ -757,7 +746,9 @@ def test_full_vs_elkan():
     km1 = KMeans(algorithm='full', random_state=13).fit(X)
     km2 = KMeans(algorithm='elkan', random_state=13).fit(X)
 
-    assert homogeneity_score(km1.predict(X), km2.predict(X)) == 1.0
+    assert homogeneity_score(
+        km1.predict(X), km2.predict(X)
+    ) == pytest.approx(1.0)
 
 
 def test_n_init():
@@ -803,10 +794,6 @@ def test_k_means_function():
     assert_warns(RuntimeWarning, k_means, X, n_clusters=n_clusters,
                  sample_weight=None, init=centers)
 
-    # to many clusters desired
-    with pytest.raises(ValueError):
-        k_means(X, n_clusters=X.shape[0] + 1, sample_weight=None)
-
 
 def test_x_squared_norms_init_centroids():
     # Test that x_squared_norms can be None in _init_centroids
@@ -820,51 +807,40 @@ def test_x_squared_norms_init_centroids():
         _init_centroids(X, 3, "k-means++", random_state=0))
 
 
-def test_max_iter_error():
-    km = KMeans(max_iter=-1)
-    assert_raise_message(ValueError, 'Number of iterations should be',
-                         km.fit, X)
-
-
-@pytest.mark.parametrize('Estimator', [KMeans, MiniBatchKMeans])
-@pytest.mark.parametrize('is_sparse', [False, True])
-def test_float_precision(Estimator, is_sparse):
-
-    estimator = Estimator(n_init=1, random_state=30)
+@pytest.mark.parametrize("data", [X, X_csr], ids=["dense", "sparse"])
+@pytest.mark.parametrize("Estimator", [KMeans, MiniBatchKMeans])
+def test_float_precision(Estimator, data):
+    # Check that the results are the same for single and double precision.
+    km = Estimator(n_init=1, random_state=0)
 
     inertia = {}
-    X_new = {}
+    Xt = {}
     centers = {}
+    labels = {}
 
     for dtype in [np.float64, np.float32]:
-        if is_sparse:
-            X_test = sp.csr_matrix(X_csr, dtype=dtype)
-        else:
-            X_test = X.astype(dtype)
-        estimator.fit(X_test)
-        # dtype of cluster centers has to be the dtype of the input
-        # data
-        assert estimator.cluster_centers_.dtype == dtype
-        inertia[dtype] = estimator.inertia_
-        X_new[dtype] = estimator.transform(X_test)
-        centers[dtype] = estimator.cluster_centers_
-        # ensure the extracted row is a 2d array
-        assert estimator.predict(X_test[:1]) == estimator.labels_[0]
-        if hasattr(estimator, 'partial_fit'):
-            estimator.partial_fit(X_test[0:3])
-            # dtype of cluster centers has to stay the same after
-            # partial_fit
-            assert estimator.cluster_centers_.dtype == dtype
-
-    # compare arrays with low precision since the difference between
-    # 32 and 64 bit sometimes makes a difference up to the 4th decimal
-    # place
-    assert_array_almost_equal(inertia[np.float32], inertia[np.float64],
-                              decimal=4)
-    assert_array_almost_equal(X_new[np.float32], X_new[np.float64],
-                              decimal=4)
-    assert_array_almost_equal(centers[np.float32], centers[np.float64],
-                              decimal=4)
+        X = data.astype(dtype, **_astype_copy_false(data))
+        km.fit(X)
+
+        inertia[dtype] = km.inertia_
+        Xt[dtype] = km.transform(X)
+        centers[dtype] = km.cluster_centers_
+        labels[dtype] = km.labels_
+
+        # dtype of cluster centers has to be the dtype of the input data
+        assert km.cluster_centers_.dtype == dtype
+
+        # same with partial_fit
+        if Estimator is MiniBatchKMeans:
+            km.partial_fit(X[0:3])
+            assert km.cluster_centers_.dtype == dtype
+
+    # compare arrays with low precision since the difference between 32 and
+    # 64 bit comes from an accumulation of rounding errors.
+    assert_allclose(inertia[np.float32], inertia[np.float64], rtol=1e-5)
+    assert_allclose(Xt[np.float32], Xt[np.float64], rtol=1e-5)
+    assert_allclose(centers[np.float32], centers[np.float64], rtol=1e-5)
+    assert_array_equal(labels[np.float32], labels[np.float64])
 
 
 def test_k_means_init_centers():
@@ -893,24 +869,6 @@ def test_k_means_init_fitted_centers(data):
     assert_array_almost_equal(centers, new_centers)
 
 
-def test_sparse_validate_centers():
-    from sklearn.datasets import load_iris
-
-    iris = load_iris()
-    X = iris.data
-
-    # Get a local optimum
-    centers = KMeans(n_clusters=4).fit(X).cluster_centers_
-
-    # Test that a ValueError is raised for validate_center_shape
-    classifier = KMeans(n_clusters=3, init=centers, n_init=1)
-
-    msg = r"The shape of the initial centers \(\(4L?, 4L?\)\) " \
-          "does not match the number of clusters 3"
-    with pytest.raises(ValueError, match=msg):
-        classifier.fit(X)
-
-
 def test_less_centers_than_unique_points():
     X = np.asarray([[0, 0],
                     [0, 1],
@@ -986,24 +944,6 @@ def test_scaled_weights():
                             _sort_centers(est_2.cluster_centers_))
 
 
-def test_sample_weight_length():
-    # check that an error is raised when passing sample weights
-    # with an incompatible shape
-    km = KMeans(n_clusters=n_clusters, random_state=42)
-    msg = r'sample_weight.shape == \(2,\), expected \(100,\)'
-    with pytest.raises(ValueError, match=msg):
-        km.fit(X, sample_weight=np.ones(2))
-
-
-def test_check_normalize_sample_weight():
-    from sklearn.cluster._kmeans import _check_normalize_sample_weight
-    sample_weight = None
-    checked_sample_weight = _check_normalize_sample_weight(sample_weight, X)
-    assert _num_samples(X) == _num_samples(checked_sample_weight)
-    assert_almost_equal(checked_sample_weight.sum(), _num_samples(X))
-    assert X.dtype == checked_sample_weight.dtype
-
-
 def test_iter_attribute():
     # Regression test on bad n_iter_ value. Previous bug n_iter_ was one off
     # it's right value (#11340).
@@ -1177,3 +1117,54 @@ def test_sample_weight_unchanged():
 
     # internally, sample_weight is rescale to sum up to n_samples = 3
     assert_array_equal(sample_weight, np.array([0.5, 0.2, 0.3]))
+
+
+@pytest.mark.parametrize("Estimator", [KMeans, MiniBatchKMeans])
+@pytest.mark.parametrize("param, match", [
+    ({"n_init": 0}, r"n_init should be > 0"),
+    ({"max_iter": 0}, r"max_iter should be > 0"),
+    ({"n_clusters": n_samples + 1}, r"n_samples.* should be >= n_clusters"),
+    ({"init": X[:2]},
+     r"The shape of the initial centers .* does not match "
+     r"the number of clusters"),
+    ({"init": lambda X_, k, random_state: X_[:2]},
+     r"The shape of the initial centers .* does not match "
+     r"the number of clusters"),
+    ({"init": X[:8, :2]},
+     r"The shape of the initial centers .* does not match "
+     r"the number of features of the data"),
+    ({"init": lambda X_, k, random_state: X_[:8, :2]},
+     r"The shape of the initial centers .* does not match "
+     r"the number of features of the data"),
+    ({"init": "wrong"},
+     r"init should be either 'k-means\+\+', 'random', "
+     r"a ndarray or a callable")]
+)
+def test_wrong_params(Estimator, param, match):
+    # Check that error are raised with clear error message when wrong values
+    # are passed for the parameters
+    with pytest.raises(ValueError, match=match):
+        Estimator(**param).fit(X)
+
+
+@pytest.mark.parametrize("param, match", [
+    ({"algorithm": "wrong"}, r"Algorithm must be 'auto', 'full' or 'elkan'")]
+)
+def test_kmeans_wrong_params(param, match):
+    # Check that error are raised with clear error message when wrong values
+    # are passed for the KMeans specific parameters
+    with pytest.raises(ValueError, match=match):
+        KMeans(**param).fit(X)
+
+
+@pytest.mark.parametrize("param, match", [
+    ({"max_no_improvement": -1}, r"max_no_improvement should be >= 0"),
+    ({"batch_size": -1}, r"batch_size should be > 0"),
+    ({"init_size": -1}, r"init_size should be > 0"),
+    ({"reassignment_ratio": -1}, r"reassignment_ratio should be >= 0")]
+)
+def test_minibatch_kmeans_wrong_params(param, match):
+    # Check that error are raised with clear error message when wrong values
+    # are passed for the MiniBatchKMeans specific parameters
+    with pytest.raises(ValueError, match=match):
+        MiniBatchKMeans(**param).fit(X)
diff --git a/sklearn/cluster/tests/test_mean_shift.py b/sklearn/cluster/tests/test_mean_shift.py
index 934f9fc8c3d62..7d2300711466a 100644
--- a/sklearn/cluster/tests/test_mean_shift.py
+++ b/sklearn/cluster/tests/test_mean_shift.py
@@ -12,12 +12,14 @@
 from sklearn.utils._testing import assert_array_equal
 from sklearn.utils._testing import assert_array_almost_equal
 from sklearn.utils._testing import assert_raise_message
+from sklearn.utils._testing import assert_allclose
 
 from sklearn.cluster import MeanShift
 from sklearn.cluster import mean_shift
 from sklearn.cluster import estimate_bandwidth
 from sklearn.cluster import get_bin_seeds
 from sklearn.datasets import make_blobs
+from sklearn.metrics import v_measure_score
 
 
 n_clusters = 3
@@ -168,3 +170,25 @@ def test_max_iter(max_iter):
 
     for c1, c2 in zip(clusters1, clusters2):
         assert np.allclose(c1, c2)
+
+
+def test_mean_shift_zero_bandwidth():
+    # Check that mean shift works when the estimated bandwidth is 0.
+    X = np.array([1, 1, 1, 2, 2, 2, 3, 3]).reshape(-1, 1)
+
+    # estimate_bandwidth with default args returns 0 on this dataset
+    bandwidth = estimate_bandwidth(X)
+    assert bandwidth == 0
+
+    # get_bin_seeds with a 0 bin_size should return the dataset itself
+    assert get_bin_seeds(X, bin_size=bandwidth) is X
+
+    # MeanShift with binning and a 0 estimated bandwidth should be equivalent
+    # to no binning.
+    ms_binning = MeanShift(bin_seeding=True, bandwidth=None).fit(X)
+    ms_nobinning = MeanShift(bin_seeding=False).fit(X)
+    expected_labels = np.array([0, 0, 0, 1, 1, 1, 2, 2])
+
+    assert v_measure_score(ms_binning.labels_, expected_labels) == 1
+    assert v_measure_score(ms_nobinning.labels_, expected_labels) == 1
+    assert_allclose(ms_binning.cluster_centers_, ms_nobinning.cluster_centers_)
diff --git a/sklearn/conftest.py b/sklearn/conftest.py
index d38e45f57b4f8..5dadb7e68469b 100644
--- a/sklearn/conftest.py
+++ b/sklearn/conftest.py
@@ -15,7 +15,7 @@ def pyplot():
         The ``matplotlib.pyplot`` module.
     """
     matplotlib = pytest.importorskip('matplotlib')
-    matplotlib.use('agg', warn=False, force=True)
+    matplotlib.use('agg')
     pyplot = pytest.importorskip('matplotlib.pyplot')
     yield pyplot
     pyplot.close('all')
diff --git a/sklearn/datasets/_lfw.py b/sklearn/datasets/_lfw.py
index b8db75010e8f2..33586fac94f44 100644
--- a/sklearn/datasets/_lfw.py
+++ b/sklearn/datasets/_lfw.py
@@ -12,7 +12,6 @@
 from os.path import dirname, join, exists, isdir
 
 import logging
-from distutils.version import LooseVersion
 
 import numpy as np
 import joblib
@@ -21,6 +20,7 @@
 from ._base import get_data_home, _fetch_remote, RemoteFileMetadata
 from ..utils import Bunch
 from ..utils.validation import _deprecate_positional_args
+from ..utils.fixes import parse_version
 
 logger = logging.getLogger(__name__)
 
@@ -304,7 +304,7 @@ def fetch_lfw_people(*, data_home=None, funneled=True, resize=0.5,
 
     # wrap the loader in a memoizing function that will return memmaped data
     # arrays for optimal memory usage
-    if LooseVersion(joblib.__version__) < LooseVersion('0.12'):
+    if parse_version(joblib.__version__) < parse_version('0.12'):
         # Deal with change of API in joblib
         m = Memory(cachedir=lfw_home, compress=6, verbose=0)
     else:
@@ -477,7 +477,7 @@ def fetch_lfw_pairs(*, subset='train', data_home=None, funneled=True,
 
     # wrap the loader in a memoizing function that will return memmaped data
     # arrays for optimal memory usage
-    if LooseVersion(joblib.__version__) < LooseVersion('0.12'):
+    if parse_version(joblib.__version__) < parse_version('0.12'):
         # Deal with change of API in joblib
         m = Memory(cachedir=lfw_home, compress=6, verbose=0)
     else:
diff --git a/sklearn/decomposition/_dict_learning.py b/sklearn/decomposition/_dict_learning.py
index 7a9650353408f..143d70bbaebd2 100644
--- a/sklearn/decomposition/_dict_learning.py
+++ b/sklearn/decomposition/_dict_learning.py
@@ -1490,7 +1490,7 @@ def partial_fit(self, X, y=None, iter_offset=None):
             iter_offset = getattr(self, 'iter_offset_', 0)
         U, (A, B) = dict_learning_online(
             X, self.n_components, alpha=self.alpha,
-            n_iter=self.n_iter, method=self.fit_algorithm,
+            n_iter=1, method=self.fit_algorithm,
             method_max_iter=self.transform_max_iter,
             n_jobs=self.n_jobs, dict_init=dict_init,
             batch_size=len(X), shuffle=False,
@@ -1504,5 +1504,5 @@ def partial_fit(self, X, y=None, iter_offset=None):
         # Keep track of the state of the algorithm to be able to do
         # some online fitting (partial_fit)
         self.inner_stats_ = (A, B)
-        self.iter_offset_ = iter_offset + self.n_iter
+        self.iter_offset_ = iter_offset + 1
         return self
diff --git a/sklearn/decomposition/_incremental_pca.py b/sklearn/decomposition/_incremental_pca.py
index bc34c17326f19..5437f2cd0b2a8 100644
--- a/sklearn/decomposition/_incremental_pca.py
+++ b/sklearn/decomposition/_incremental_pca.py
@@ -290,8 +290,8 @@ def partial_fit(self, X, y=None, check_input=True):
             X -= col_batch_mean
             # Build matrix of combined previous basis and new data
             mean_correction = \
-                np.sqrt((self.n_samples_seen_ * n_samples) /
-                        n_total_samples) * (self.mean_ - col_batch_mean)
+                np.sqrt((self.n_samples_seen_ / n_total_samples) *
+                        n_samples) * (self.mean_ - col_batch_mean)
             X = np.vstack((self.singular_values_.reshape((-1, 1)) *
                            self.components_, X, mean_correction))
 
diff --git a/sklearn/decomposition/tests/test_dict_learning.py b/sklearn/decomposition/tests/test_dict_learning.py
index f733ca75c9d6e..547dc175ee26d 100644
--- a/sklearn/decomposition/tests/test_dict_learning.py
+++ b/sklearn/decomposition/tests/test_dict_learning.py
@@ -391,6 +391,23 @@ def test_dict_learning_online_partial_fit():
                               decimal=2)
 
 
+def test_dict_learning_iter_offset():
+    n_components = 12
+    rng = np.random.RandomState(0)
+    V = rng.randn(n_components, n_features)
+    dict1 = MiniBatchDictionaryLearning(n_components, n_iter=10,
+                                        dict_init=V, random_state=0,
+                                        shuffle=False)
+    dict2 = MiniBatchDictionaryLearning(n_components, n_iter=10,
+                                        dict_init=V, random_state=0,
+                                        shuffle=False)
+    dict1.fit(X)
+    for sample in X:
+        dict2.partial_fit(sample[np.newaxis, :])
+
+    assert dict1.iter_offset_ == dict2.iter_offset_
+
+
 def test_sparse_encode_shapes():
     n_components = 12
     rng = np.random.RandomState(0)
diff --git a/sklearn/decomposition/tests/test_incremental_pca.py b/sklearn/decomposition/tests/test_incremental_pca.py
index 632c02b6c8618..d198b67c720c1 100644
--- a/sklearn/decomposition/tests/test_incremental_pca.py
+++ b/sklearn/decomposition/tests/test_incremental_pca.py
@@ -384,3 +384,18 @@ def test_incremental_pca_partial_fit_float_division():
 
     np.testing.assert_allclose(singular_vals_float_samples_seen,
                                singular_vals_int_samples_seen)
+
+
+def test_incremental_pca_fit_overflow_error():
+    # Test for overflow error on Windows OS
+    # (non-regression test for issue #17693)
+    rng = np.random.RandomState(0)
+    A = rng.rand(500000, 2)
+
+    ipca = IncrementalPCA(n_components=2, batch_size=10000)
+    ipca.fit(A)
+
+    pca = PCA(n_components=2)
+    pca.fit(A)
+
+    np.testing.assert_allclose(ipca.singular_values_, pca.singular_values_)
diff --git a/sklearn/ensemble/_gb_losses.py b/sklearn/ensemble/_gb_losses.py
index 7bd5faca1d7d9..28bab32ff1c31 100644
--- a/sklearn/ensemble/_gb_losses.py
+++ b/sklearn/ensemble/_gb_losses.py
@@ -710,13 +710,11 @@ def __call__(self, y, raw_predictions, sample_weight=None):
         for k in range(self.K):
             Y[:, k] = y == k
 
-        if sample_weight is None:
-            return np.sum(-1 * (Y * raw_predictions).sum(axis=1) +
-                          logsumexp(raw_predictions, axis=1))
-        else:
-            return np.sum(
-                -1 * sample_weight * (Y * raw_predictions).sum(axis=1) +
-                logsumexp(raw_predictions, axis=1))
+        return np.average(
+            -1 * (Y * raw_predictions).sum(axis=1) +
+            logsumexp(raw_predictions, axis=1),
+            weights=sample_weight
+        )
 
     def negative_gradient(self, y, raw_predictions, k=0, **kwargs):
         """Compute negative gradient for the ``k``-th class.
diff --git a/sklearn/ensemble/_hist_gradient_boosting/tests/test_loss.py b/sklearn/ensemble/_hist_gradient_boosting/tests/test_loss.py
index c3f6ded7be39a..adaac4946bb5b 100644
--- a/sklearn/ensemble/_hist_gradient_boosting/tests/test_loss.py
+++ b/sklearn/ensemble/_hist_gradient_boosting/tests/test_loss.py
@@ -3,7 +3,7 @@
 from numpy.testing import assert_allclose
 from scipy.optimize import newton
 from sklearn.utils import assert_all_finite
-from sklearn.utils.fixes import sp_version
+from sklearn.utils.fixes import sp_version, parse_version
 import pytest
 
 from sklearn.ensemble._hist_gradient_boosting.loss import _LOSSES
@@ -57,7 +57,7 @@ def get_hessians(y_true, raw_predictions):
     ('poisson', 0., 2.),
     ('poisson', -22., 10.),
 ])
-@pytest.mark.skipif(sp_version == (1, 2, 0),
+@pytest.mark.skipif(sp_version == parse_version('1.2.0'),
                     reason='bug in scipy 1.2.0, see scipy issue #9608')
 @skip_if_32bit
 def test_derivatives(loss, x0, y_true):
diff --git a/sklearn/ensemble/_stacking.py b/sklearn/ensemble/_stacking.py
index 73aa55c0575a7..3d08d1e3a9eaf 100644
--- a/sklearn/ensemble/_stacking.py
+++ b/sklearn/ensemble/_stacking.py
@@ -13,6 +13,7 @@
 from ..base import clone
 from ..base import ClassifierMixin, RegressorMixin, TransformerMixin
 from ..base import is_classifier, is_regressor
+from ..exceptions import NotFittedError
 from ..utils._estimator_html_repr import _VisualBlock
 
 from ._base import _fit_single_estimator
@@ -146,7 +147,6 @@ def fit(self, X, y, sample_weight=None):
             delayed(_fit_single_estimator)(clone(est), X, y, sample_weight)
             for est in all_estimators if est != 'drop'
         )
-        self.n_features_in_ = self.estimators_[0].n_features_in_
 
         self.named_estimators_ = Bunch()
         est_fitted_idx = 0
@@ -197,6 +197,17 @@ def fit(self, X, y, sample_weight=None):
 
         return self
 
+    @property
+    def n_features_in_(self):
+        """Number of features seen during :term:`fit`."""
+        try:
+            check_is_fitted(self)
+        except NotFittedError as nfe:
+            raise AttributeError(
+                f"{self.__class__.__name__} object has no attribute "
+                f"n_features_in_") from nfe
+        return self.estimators_[0].n_features_in_
+
     def _transform(self, X):
         """Concatenate and return the predictions of the estimators."""
         check_is_fitted(self)
diff --git a/sklearn/ensemble/tests/test_forest.py b/sklearn/ensemble/tests/test_forest.py
index 775ed851d5a53..17964da691cf7 100644
--- a/sklearn/ensemble/tests/test_forest.py
+++ b/sklearn/ensemble/tests/test_forest.py
@@ -11,7 +11,6 @@
 import pickle
 import math
 from collections import defaultdict
-from distutils.version import LooseVersion
 import itertools
 from itertools import combinations
 from itertools import product
@@ -35,6 +34,7 @@
 from sklearn.utils._testing import assert_warns_message
 from sklearn.utils._testing import ignore_warnings
 from sklearn.utils._testing import skip_if_no_parallel
+from sklearn.utils.fixes import parse_version
 
 from sklearn.exceptions import NotFittedError
 
@@ -1274,7 +1274,7 @@ def start_call(self):
 joblib.register_parallel_backend('testing', MyBackend)
 
 
-@pytest.mark.skipif(joblib.__version__ < LooseVersion('0.12'),
+@pytest.mark.skipif(parse_version(joblib.__version__) < parse_version('0.12'),
                     reason='tests not yet supported in joblib <0.12')
 @skip_if_no_parallel
 def test_backend_respected():
diff --git a/sklearn/ensemble/tests/test_gradient_boosting_loss_functions.py b/sklearn/ensemble/tests/test_gradient_boosting_loss_functions.py
index 6b24f90d0239d..dae708f297906 100644
--- a/sklearn/ensemble/tests/test_gradient_boosting_loss_functions.py
+++ b/sklearn/ensemble/tests/test_gradient_boosting_loss_functions.py
@@ -27,7 +27,7 @@ def test_binomial_deviance():
 
     # pred has the same BD for y in {0, 1}
     assert (bd(np.array([0.0]), np.array([0.0])) ==
-                 bd(np.array([1.0]), np.array([0.0])))
+            bd(np.array([1.0]), np.array([0.0])))
 
     assert_almost_equal(bd(np.array([1.0, 1.0, 1.0]),
                            np.array([100.0, 100.0, 100.0])),
@@ -36,8 +36,9 @@ def test_binomial_deviance():
                            np.array([100.0, -100.0, -100.0])), 0)
 
     # check if same results as alternative definition of deviance (from ESLII)
-    alt_dev = lambda y, pred: np.mean(np.logaddexp(0.0, -2.0 *
-                                                   (2.0 * y - 1) * pred))
+    def alt_dev(y, pred):
+        return np.mean(np.logaddexp(0.0, -2.0 * (2.0 * y - 1) * pred))
+
     test_data = [(np.array([1.0, 1.0, 1.0]), np.array([100.0, 100.0, 100.0])),
                  (np.array([0.0, 0.0, 0.0]), np.array([100.0, 100.0, 100.0])),
                  (np.array([0.0, 0.0, 0.0]),
@@ -49,7 +50,9 @@ def test_binomial_deviance():
         assert_almost_equal(bd(*datum), alt_dev(*datum))
 
     # check the gradient against the
-    alt_ng = lambda y, pred: (2 * y - 1) / (1 + np.exp(2 * (2 * y - 1) * pred))
+    def alt_ng(y, pred):
+        return (2 * y - 1) / (1 + np.exp(2 * (2 * y - 1) * pred))
+
     for datum in test_data:
         assert_almost_equal(bd.negative_gradient(*datum), alt_ng(*datum))
 
@@ -176,6 +179,50 @@ def test_sample_weight_deviance():
         assert deviance_wo_w == deviance_w_w
 
 
+@pytest.mark.parametrize(
+    'n_classes, n_samples', [(3, 100), (5, 57), (7, 13)]
+)
+def test_multinomial_deviance(n_classes, n_samples):
+    # Check multinomial deviance with and without sample weights.
+    rng = np.random.RandomState(13)
+    sample_weight = np.ones(n_samples)
+    y_true = rng.randint(0, n_classes, size=n_samples)
+    y_pred = np.zeros((n_samples, n_classes), dtype=np.float64)
+    for klass in range(y_pred.shape[1]):
+        y_pred[:, klass] = y_true == klass
+
+    loss = MultinomialDeviance(n_classes)
+    loss_wo_sw = loss(y_true, y_pred)
+    assert loss_wo_sw > 0
+    loss_w_sw = loss(y_true, y_pred, sample_weight=sample_weight)
+    assert loss_wo_sw == pytest.approx(loss_w_sw)
+
+    # Multinomial deviance uses weighted average loss rather than
+    # weighted sum loss, so we make sure that the value remains the same
+    # when we device the weight by 2.
+    loss_w_sw = loss(y_true, y_pred, sample_weight=0.5 * sample_weight)
+    assert loss_wo_sw == pytest.approx(loss_w_sw)
+
+
+def test_mdl_computation_weighted():
+    raw_predictions = np.array([[1., -1., -.1], [-2., 1., 2.]])
+    y_true = np.array([0, 1])
+    weights = np.array([1, 3])
+    expected_loss = 1.0909323
+    # MultinomialDeviance loss computation with weights.
+    loss = MultinomialDeviance(3)
+    assert (loss(y_true, raw_predictions, weights)
+            == pytest.approx(expected_loss))
+
+
+@pytest.mark.parametrize('n', [0, 1, 2])
+def test_mdl_exception(n):
+    # Check that MultinomialDeviance throws an exception when n_classes <= 2
+    err_msg = 'MultinomialDeviance requires more than 2 classes.'
+    with pytest.raises(ValueError, match=err_msg):
+        MultinomialDeviance(n)
+
+
 def test_init_raw_predictions_shapes():
     # Make sure get_init_raw_predictions returns float64 arrays with shape
     # (n_samples, K) where K is 1 for binary classification and regression, and
diff --git a/sklearn/ensemble/tests/test_stacking.py b/sklearn/ensemble/tests/test_stacking.py
index f8a3f290e96b5..a7b02a3571c11 100644
--- a/sklearn/ensemble/tests/test_stacking.py
+++ b/sklearn/ensemble/tests/test_stacking.py
@@ -17,6 +17,8 @@
 from sklearn.datasets import load_iris
 from sklearn.datasets import load_diabetes
 from sklearn.datasets import load_breast_cancer
+from sklearn.datasets import make_regression
+from sklearn.datasets import make_classification
 
 from sklearn.dummy import DummyClassifier
 from sklearn.dummy import DummyRegressor
@@ -491,3 +493,32 @@ def test_stacking_cv_influence(stacker, X, y):
     with pytest.raises(AssertionError, match='Not equal'):
         assert_allclose(stacker_cv_3.final_estimator_.coef_,
                         stacker_cv_5.final_estimator_.coef_)
+
+
+@pytest.mark.parametrize("make_dataset, Stacking, Estimator", [
+    (make_classification, StackingClassifier, LogisticRegression),
+    (make_regression, StackingRegressor, LinearRegression)
+])
+def test_stacking_without_n_features_in(make_dataset, Stacking, Estimator):
+    # Stacking supports estimators without `n_features_in_`. Regression test
+    # for #17353
+
+    class MyEstimator(Estimator):
+        """Estimator without n_features_in_"""
+        def fit(self, X, y):
+            super().fit(X, y)
+            del self.n_features_in_
+
+    X, y = make_dataset(random_state=0, n_samples=100)
+    stacker = Stacking(estimators=[('lr', MyEstimator())])
+
+    msg = f"{Stacking.__name__} object has no attribute n_features_in_"
+    with pytest.raises(AttributeError, match=msg):
+        stacker.n_features_in_
+
+    # Does not raise
+    stacker.fit(X, y)
+
+    msg = "'MyEstimator' object has no attribute 'n_features_in_'"
+    with pytest.raises(AttributeError, match=msg):
+        stacker.n_features_in_
diff --git a/sklearn/feature_extraction/tests/test_text.py b/sklearn/feature_extraction/tests/test_text.py
index 86ae2fd6c149e..477008041dd99 100644
--- a/sklearn/feature_extraction/tests/test_text.py
+++ b/sklearn/feature_extraction/tests/test_text.py
@@ -1363,3 +1363,12 @@ class MyVectorizer(VectorizerMixin):
            "in version 0.24.")
     with pytest.warns(FutureWarning, match=msg):
         MyVectorizer()
+
+
+def test_tie_breaking_sample_order_invariance():
+    # Checks the sample order invariance when setting max_features
+    # non-regression test for #17939
+    vec = CountVectorizer(max_features=1)
+    vocab1 = vec.fit(['hello', 'world']).vocabulary_
+    vocab2 = vec.fit(['world', 'hello']).vocabulary_
+    assert vocab1 == vocab2
diff --git a/sklearn/feature_extraction/text.py b/sklearn/feature_extraction/text.py
index 661f638b000fc..a98bcbc4303b2 100644
--- a/sklearn/feature_extraction/text.py
+++ b/sklearn/feature_extraction/text.py
@@ -1212,13 +1212,14 @@ def fit_transform(self, raw_documents, y=None):
             if max_doc_count < min_doc_count:
                 raise ValueError(
                     "max_df corresponds to < documents than min_df")
+            if max_features is not None:
+                X = self._sort_features(X, vocabulary)
             X, self.stop_words_ = self._limit_features(X, vocabulary,
                                                        max_doc_count,
                                                        min_doc_count,
                                                        max_features)
-
-            X = self._sort_features(X, vocabulary)
-
+            if max_features is None:
+                X = self._sort_features(X, vocabulary)
             self.vocabulary_ = vocabulary
 
         return X
diff --git a/sklearn/linear_model/_least_angle.py b/sklearn/linear_model/_least_angle.py
index 255baacea9a59..c1ee3a738e003 100644
--- a/sklearn/linear_model/_least_angle.py
+++ b/sklearn/linear_model/_least_angle.py
@@ -411,12 +411,6 @@ def _lars_path_solver(X, y, Xy=None, Gram=None, n_samples=None, max_iter=500,
         Gram = None
         if X is None:
             raise ValueError('X and Gram cannot both be unspecified.')
-        if copy_X:
-            # force copy. setting the array to be fortran-ordered
-            # speeds up the calculation of the (partial) Gram matrix
-            # and allows to easily swap columns
-            X = X.copy('F')
-
     elif isinstance(Gram, str) and Gram == 'auto' or Gram is True:
         if Gram is True or X.shape[0] > X.shape[1]:
             Gram = np.dot(X.T, X)
@@ -432,6 +426,13 @@ def _lars_path_solver(X, y, Xy=None, Gram=None, n_samples=None, max_iter=500,
         if Gram.shape != (n_features, n_features):
             raise ValueError('The shapes of the inputs Gram and Xy'
                              ' do not match.')
+
+    if copy_X and X is not None and Gram is None:
+        # force copy. setting the array to be fortran-ordered
+        # speeds up the calculation of the (partial) Gram matrix
+        # and allows to easily swap columns
+        X = X.copy('F')
+
     max_features = min(max_iter, n_features)
 
     if return_path:
diff --git a/sklearn/linear_model/tests/test_base.py b/sklearn/linear_model/tests/test_base.py
index 0d30c4dd13022..af76826715241 100644
--- a/sklearn/linear_model/tests/test_base.py
+++ b/sklearn/linear_model/tests/test_base.py
@@ -5,8 +5,6 @@
 
 import pytest
 
-from distutils.version import LooseVersion
-
 import numpy as np
 from scipy import sparse
 from scipy import linalg
@@ -15,6 +13,7 @@
 from sklearn.utils._testing import assert_array_equal
 from sklearn.utils._testing import assert_almost_equal
 from sklearn.utils._testing import assert_allclose
+from sklearn.utils.fixes import parse_version
 
 from sklearn.linear_model import LinearRegression
 from sklearn.linear_model._base import _preprocess_data
@@ -210,7 +209,7 @@ def test_linear_regression_sparse_multiple_outcome(random_state=0):
 def test_linear_regression_pd_sparse_dataframe_warning():
     pd = pytest.importorskip('pandas')
     # restrict the pd versions < '0.24.0' as they have a bug in is_sparse func
-    if LooseVersion(pd.__version__) < '0.24.0':
+    if parse_version(pd.__version__) < parse_version('0.24.0'):
         pytest.skip("pandas 0.24+ required.")
 
     # Warning is raised only when some of the columns is sparse
diff --git a/sklearn/linear_model/tests/test_coordinate_descent.py b/sklearn/linear_model/tests/test_coordinate_descent.py
index 1b2f7c656f015..9213c7d8e3193 100644
--- a/sklearn/linear_model/tests/test_coordinate_descent.py
+++ b/sklearn/linear_model/tests/test_coordinate_descent.py
@@ -8,7 +8,6 @@
 from scipy import interpolate, sparse
 from copy import deepcopy
 import joblib
-from distutils.version import LooseVersion
 
 from sklearn.datasets import load_boston
 from sklearn.datasets import make_regression
@@ -23,6 +22,7 @@
 from sklearn.utils._testing import ignore_warnings
 from sklearn.utils._testing import assert_array_equal
 from sklearn.utils._testing import TempMemmap
+from sklearn.utils.fixes import parse_version
 
 from sklearn.linear_model import Lasso, \
     LassoCV, ElasticNet, ElasticNetCV, MultiTaskLasso, MultiTaskElasticNet, \
@@ -1034,7 +1034,8 @@ def test_linear_models_cv_fit_for_all_backends(backend, estimator):
     # memmapped when using loky backend, causing an error due to unexpected
     # behavior of fancy indexing of read-only memmaps (cf. numpy#14132).
 
-    if joblib.__version__ < LooseVersion('0.12') and backend == 'loky':
+    if (parse_version(joblib.__version__) < parse_version('0.12')
+            and backend == 'loky'):
         pytest.skip('loky backend does not exist in joblib <0.12')
 
     # Create a problem sufficiently large to cause memmapping (1MB).
diff --git a/sklearn/linear_model/tests/test_least_angle.py b/sklearn/linear_model/tests/test_least_angle.py
index e198dfb15e323..fc4005b33952c 100644
--- a/sklearn/linear_model/tests/test_least_angle.py
+++ b/sklearn/linear_model/tests/test_least_angle.py
@@ -1,7 +1,5 @@
 import warnings
 
-from distutils.version import LooseVersion
-
 import numpy as np
 import pytest
 from scipy import linalg
@@ -14,6 +12,7 @@
 from sklearn.utils._testing import ignore_warnings
 from sklearn.utils._testing import assert_warns
 from sklearn.utils._testing import TempMemmap
+from sklearn.utils.fixes import np_version, parse_version
 from sklearn.exceptions import ConvergenceWarning
 from sklearn import linear_model, datasets
 from sklearn.linear_model._least_angle import _lars_path_residues
@@ -123,7 +122,7 @@ def test_lars_lstsq():
     clf = linear_model.LassoLars(alpha=0.)
     clf.fit(X1, y)
     # Avoid FutureWarning about default value change when numpy >= 1.14
-    rcond = None if LooseVersion(np.__version__) >= '1.14' else -1
+    rcond = None if np_version >= parse_version('1.14') else -1
     coef_lstsq = np.linalg.lstsq(X1, y, rcond=rcond)[0]
     assert_array_almost_equal(clf.coef_, coef_lstsq)
 
@@ -761,3 +760,16 @@ def test_X_none_gram_not_none():
     with pytest.raises(ValueError,
                        match="X cannot be None if Gram is not None"):
         lars_path(X=None, y=[1], Gram='not None')
+
+
+def test_copy_X_with_auto_gram():
+    # Non-regression test for #17789, `copy_X=True` and Gram='auto' does not
+    # overwrite X
+    rng = np.random.RandomState(42)
+    X = rng.rand(6, 6)
+    y = rng.rand(6)
+
+    X_before = X.copy()
+    linear_model.lars_path(X, y, Gram='auto', copy_X=True, method='lasso')
+    # X did not change
+    assert_allclose(X, X_before)
diff --git a/sklearn/linear_model/tests/test_sgd.py b/sklearn/linear_model/tests/test_sgd.py
index 9246413112558..b487d6b36c3c0 100644
--- a/sklearn/linear_model/tests/test_sgd.py
+++ b/sklearn/linear_model/tests/test_sgd.py
@@ -1,4 +1,3 @@
-from distutils.version import LooseVersion
 import pickle
 import pytest
 
@@ -13,6 +12,7 @@
 from sklearn.utils._testing import assert_raises_regexp
 from sklearn.utils._testing import assert_warns
 from sklearn.utils._testing import ignore_warnings
+from sklearn.utils.fixes import parse_version
 
 from sklearn import linear_model, datasets, metrics
 from sklearn.base import clone, is_classifier
@@ -1594,7 +1594,8 @@ def test_SGDClassifier_fit_for_all_backends(backend):
     # a segmentation fault when trying to write in a readonly memory mapped
     # buffer.
 
-    if joblib.__version__ < LooseVersion('0.12') and backend == 'loky':
+    if (parse_version(joblib.__version__) < parse_version('0.12')
+            and backend == 'loky'):
         pytest.skip('loky backend does not exist in joblib <0.12')
 
     random_state = np.random.RandomState(42)
diff --git a/sklearn/manifold/tests/test_spectral_embedding.py b/sklearn/manifold/tests/test_spectral_embedding.py
index 3d9d87af3a09a..acd1ae4729eb8 100644
--- a/sklearn/manifold/tests/test_spectral_embedding.py
+++ b/sklearn/manifold/tests/test_spectral_embedding.py
@@ -117,7 +117,8 @@ def test_spectral_embedding_two_components(seed=36):
         se_precomp.fit_transform(affinity.astype(np.float32))
     # thresholding on the first components using 0.
     label_ = np.array(embedded_coordinate.ravel() < 0, dtype="float")
-    assert normalized_mutual_info_score(true_label, label_) == 1.0
+    assert normalized_mutual_info_score(
+        true_label, label_) == pytest.approx(1.0)
 
 
 @pytest.mark.parametrize("X", [S, sparse.csr_matrix(S)],
diff --git a/sklearn/metrics/_regression.py b/sklearn/metrics/_regression.py
index fd89d32a07c29..9e3e0b012a82f 100644
--- a/sklearn/metrics/_regression.py
+++ b/sklearn/metrics/_regression.py
@@ -244,6 +244,8 @@ def mean_squared_error(y_true, y_pred, *,
     >>> y_pred = [[0, 2],[-1, 2],[8, -5]]
     >>> mean_squared_error(y_true, y_pred)
     0.708...
+    >>> mean_squared_error(y_true, y_pred, squared=False)
+    0.822...
     >>> mean_squared_error(y_true, y_pred, multioutput='raw_values')
     array([0.41666667, 1.        ])
     >>> mean_squared_error(y_true, y_pred, multioutput=[0.3, 0.7])
@@ -255,15 +257,18 @@ def mean_squared_error(y_true, y_pred, *,
     check_consistent_length(y_true, y_pred, sample_weight)
     output_errors = np.average((y_true - y_pred) ** 2, axis=0,
                                weights=sample_weight)
+
+    if not squared:
+        output_errors = np.sqrt(output_errors)
+
     if isinstance(multioutput, str):
         if multioutput == 'raw_values':
-            return output_errors if squared else np.sqrt(output_errors)
+            return output_errors
         elif multioutput == 'uniform_average':
             # pass None as weights to np.average: uniform mean
             multioutput = None
 
-    mse = np.average(output_errors, weights=multioutput)
-    return mse if squared else np.sqrt(mse)
+    return np.average(output_errors, weights=multioutput)
 
 
 @_deprecate_positional_args
diff --git a/sklearn/metrics/cluster/tests/test_supervised.py b/sklearn/metrics/cluster/tests/test_supervised.py
index 65b585894d075..5fd02e1e49770 100644
--- a/sklearn/metrics/cluster/tests/test_supervised.py
+++ b/sklearn/metrics/cluster/tests/test_supervised.py
@@ -61,13 +61,13 @@ def test_generalized_average():
 @ignore_warnings(category=FutureWarning)
 def test_perfect_matches():
     for score_func in score_funcs:
-        assert score_func([], []) == 1.0
-        assert score_func([0], [1]) == 1.0
-        assert score_func([0, 0, 0], [0, 0, 0]) == 1.0
-        assert score_func([0, 1, 0], [42, 7, 42]) == 1.0
-        assert score_func([0., 1., 0.], [42., 7., 42.]) == 1.0
-        assert score_func([0., 1., 2.], [42., 7., 2.]) == 1.0
-        assert score_func([0, 1, 2], [42, 7, 2]) == 1.0
+        assert score_func([], []) == pytest.approx(1.0)
+        assert score_func([0], [1]) == pytest.approx(1.0)
+        assert score_func([0, 0, 0], [0, 0, 0]) == pytest.approx(1.0)
+        assert score_func([0, 1, 0], [42, 7, 42]) == pytest.approx(1.0)
+        assert score_func([0., 1., 0.], [42., 7., 42.]) == pytest.approx(1.0)
+        assert score_func([0., 1., 2.], [42., 7., 2.]) == pytest.approx(1.0)
+        assert score_func([0, 1, 2], [42, 7, 2]) == pytest.approx(1.0)
     score_funcs_with_changing_means = [
         normalized_mutual_info_score,
         adjusted_mutual_info_score,
@@ -75,13 +75,17 @@ def test_perfect_matches():
     means = {"min", "geometric", "arithmetic", "max"}
     for score_func in score_funcs_with_changing_means:
         for mean in means:
-            assert score_func([], [], mean) == 1.0
-            assert score_func([0], [1], mean) == 1.0
-            assert score_func([0, 0, 0], [0, 0, 0], mean) == 1.0
-            assert score_func([0, 1, 0], [42, 7, 42], mean) == 1.0
-            assert score_func([0., 1., 0.], [42., 7., 42.], mean) == 1.0
-            assert score_func([0., 1., 2.], [42., 7., 2.], mean) == 1.0
-            assert score_func([0, 1, 2], [42, 7, 2], mean) == 1.0
+            assert score_func([], [], mean) == pytest.approx(1.0)
+            assert score_func([0], [1], mean) == pytest.approx(1.0)
+            assert score_func([0, 0, 0], [0, 0, 0], mean) == pytest.approx(1.0)
+            assert score_func(
+                [0, 1, 0], [42, 7, 42], mean) == pytest.approx(1.0)
+            assert score_func(
+                [0., 1., 0.], [42., 7., 42.], mean) == pytest.approx(1.0)
+            assert score_func(
+                [0., 1., 2.], [42., 7., 2.], mean) == pytest.approx(1.0)
+            assert score_func(
+                [0, 1, 2], [42, 7, 2], mean) == pytest.approx(1.0)
 
 
 def test_homogeneous_but_not_complete_labeling():
@@ -212,7 +216,7 @@ def test_adjusted_mutual_info_score():
     ami = adjusted_mutual_info_score(labels_a, labels_b)
     assert_almost_equal(ami, 0.27821, 5)
     ami = adjusted_mutual_info_score([1, 1, 2, 2], [2, 2, 3, 3])
-    assert ami == 1.0
+    assert ami == pytest.approx(1.0)
     # Test with a very large array
     a110 = np.array([list(labels_a) * 110]).flatten()
     b110 = np.array([list(labels_b) * 110]).flatten()
@@ -272,15 +276,19 @@ def test_exactly_zero_info_score():
     for i in np.logspace(1, 4, 4).astype(np.int):
         labels_a, labels_b = (np.ones(i, dtype=np.int),
                               np.arange(i, dtype=np.int))
-        assert normalized_mutual_info_score(labels_a, labels_b) == 0.0
-        assert v_measure_score(labels_a, labels_b) == 0.0
-        assert adjusted_mutual_info_score(labels_a, labels_b) == 0.0
-        assert normalized_mutual_info_score(labels_a, labels_b) == 0.0
+        assert normalized_mutual_info_score(
+            labels_a, labels_b) == pytest.approx(0.0)
+        assert v_measure_score(
+            labels_a, labels_b) == pytest.approx(0.0)
+        assert adjusted_mutual_info_score(
+            labels_a, labels_b) == pytest.approx(0.0)
+        assert normalized_mutual_info_score(
+            labels_a, labels_b) == pytest.approx(0.0)
         for method in ["min", "geometric", "arithmetic", "max"]:
-            assert adjusted_mutual_info_score(labels_a, labels_b,
-                                              method) == 0.0
-            assert normalized_mutual_info_score(labels_a, labels_b,
-                                                method) == 0.0
+            assert adjusted_mutual_info_score(
+                labels_a, labels_b,  method) == pytest.approx(0.0)
+            assert normalized_mutual_info_score(
+                labels_a, labels_b, method) == pytest.approx(0.0)
 
 
 def test_v_measure_and_mutual_information(seed=36):
diff --git a/sklearn/metrics/pairwise.py b/sklearn/metrics/pairwise.py
index 2e1332d18a20c..24c658f51b647 100644
--- a/sklearn/metrics/pairwise.py
+++ b/sklearn/metrics/pairwise.py
@@ -29,6 +29,7 @@
 from ..preprocessing import normalize
 from ..utils._mask import _get_mask
 from ..utils.validation import _deprecate_positional_args
+from ..utils.fixes import sp_version, parse_version
 
 from ._pairwise_fast import _chi2_kernel_fast, _sparse_manhattan
 from ..exceptions import DataConversionWarning
@@ -1441,13 +1442,16 @@ def _precompute_metric_params(X, Y, metric=None, **kwds):
     """Precompute data-derived metric parameters if not provided
     """
     if metric == "seuclidean" and 'V' not in kwds:
+        # There is a bug in scipy < 1.5 that will cause a crash if
+        # X.dtype != np.double (float64). See PR #15730
+        dtype = np.float64 if sp_version < parse_version('1.5') else None
         if X is Y:
-            V = np.var(X, axis=0, ddof=1)
+            V = np.var(X, axis=0, ddof=1, dtype=dtype)
         else:
             warnings.warn("from version 0.25, pairwise_distances for "
                           "metric='seuclidean' will require V to be "
                           "specified if Y is passed.", FutureWarning)
-            V = np.var(np.vstack([X, Y]), axis=0, ddof=1)
+            V = np.var(np.vstack([X, Y]), axis=0, ddof=1, dtype=dtype)
         return {'V': V}
     if metric == "mahalanobis" and 'VI' not in kwds:
         if X is Y:
diff --git a/sklearn/metrics/tests/test_pairwise.py b/sklearn/metrics/tests/test_pairwise.py
index d7a96de12c9e3..b03356cfcb8e9 100644
--- a/sklearn/metrics/tests/test_pairwise.py
+++ b/sklearn/metrics/tests/test_pairwise.py
@@ -681,7 +681,7 @@ def test_euclidean_distances(dtype, x_array_constr, y_array_constr):
     distances = euclidean_distances(X, Y)
 
     # the default rtol=1e-7 is too close to the float32 precision
-    # and fails due too rounding errors.
+    # and fails due to rounding errors.
     assert_allclose(distances, expected, rtol=1e-6)
     assert distances.dtype == dtype
 
@@ -702,7 +702,7 @@ def test_euclidean_distances_sym(dtype, x_array_constr):
     distances = euclidean_distances(X)
 
     # the default rtol=1e-7 is too close to the float32 precision
-    # and fails due too rounding errors.
+    # and fails due to rounding errors.
     assert_allclose(distances, expected, rtol=1e-6)
     assert distances.dtype == dtype
 
@@ -729,7 +729,7 @@ def test_euclidean_distances_upcast(batch_size, x_array_constr,
     distances = np.sqrt(np.maximum(distances, 0))
 
     # the default rtol=1e-7 is too close to the float32 precision
-    # and fails due too rounding errors.
+    # and fails due to rounding errors.
     assert_allclose(distances, expected, rtol=1e-6)
 
 
@@ -749,7 +749,7 @@ def test_euclidean_distances_upcast_sym(batch_size, x_array_constr):
     distances = np.sqrt(np.maximum(distances, 0))
 
     # the default rtol=1e-7 is too close to the float32 precision
-    # and fails due too rounding errors.
+    # and fails due to rounding errors.
     assert_allclose(distances, expected, rtol=1e-6)
 
 
@@ -1294,3 +1294,47 @@ def test_pairwise_distances_data_derived_params(n_jobs, metric, dist_function,
 
         assert_allclose(dist, expected_dist_explicit_params)
         assert_allclose(dist, expected_dist_default_params)
+
+
+@pytest.mark.parametrize(
+        'metric', [
+            'braycurtis', 'canberra', 'chebyshev',
+            'correlation', 'hamming', 'mahalanobis', 'minkowski', 'seuclidean',
+            'sqeuclidean', 'cityblock', 'cosine', 'euclidean'])
+@pytest.mark.parametrize(
+        "dtype",
+        [np.float32, np.float64])
+@pytest.mark.parametrize("y_is_x", [True, False], ids=["Y is X", "Y is not X"])
+def test_numeric_pairwise_distances_datatypes(metric, dtype, y_is_x):
+    # Check that pairwise distances gives the same result as pdist and cdist
+    # regardless of input datatype when using any scipy metric for comparing
+    # numeric vectors
+    #
+    # This test is necessary because pairwise_distances used to throw an
+    # error when using metric='seuclidean' and the input data was not
+    # of type np.float64 (#15730)
+
+    rng = np.random.RandomState(0)
+
+    X = rng.random_sample((5, 4)).astype(dtype)
+
+    params = {}
+    if y_is_x:
+        Y = X
+        expected_dist = squareform(pdist(X, metric=metric))
+    else:
+        Y = rng.random_sample((5, 4)).astype(dtype)
+        expected_dist = cdist(X, Y, metric=metric)
+        # precompute parameters for seuclidean & mahalanobis when x is not y
+        if metric == 'seuclidean':
+            params = {'V': np.var(np.vstack([X, Y]),
+                                  axis=0, ddof=1, dtype=np.float64)}
+        elif metric == 'mahalanobis':
+            params = {'VI': np.linalg.inv(np.cov(np.vstack([X, Y]).T)).T}
+
+    dist = pairwise_distances(X, Y, metric=metric, **params)
+
+    # the default rtol=1e-7 is too close to the float32 precision
+    # and fails due to rounding errors
+    rtol = 1e-5 if dtype is np.float32 else 1e-7
+    assert_allclose(dist, expected_dist, rtol=rtol)
diff --git a/sklearn/metrics/tests/test_regression.py b/sklearn/metrics/tests/test_regression.py
index 06c44b2b6f59e..c5e9743539612 100644
--- a/sklearn/metrics/tests/test_regression.py
+++ b/sklearn/metrics/tests/test_regression.py
@@ -76,7 +76,7 @@ def test_multioutput_regression():
     assert_almost_equal(error, (1. / 3 + 2. / 3 + 2. / 3) / 4.)
 
     error = mean_squared_error(y_true, y_pred, squared=False)
-    assert_almost_equal(error, 0.645, decimal=2)
+    assert_almost_equal(error, 0.454, decimal=2)
 
     error = mean_squared_log_error(y_true, y_pred)
     assert_almost_equal(error, 0.200, decimal=2)
@@ -258,7 +258,7 @@ def test_regression_custom_weights():
     evsw = explained_variance_score(y_true, y_pred, multioutput=[0.4, 0.6])
 
     assert_almost_equal(msew, 0.39, decimal=2)
-    assert_almost_equal(rmsew, 0.62, decimal=2)
+    assert_almost_equal(rmsew, 0.59, decimal=2)
     assert_almost_equal(maew, 0.475, decimal=3)
     assert_almost_equal(rw, 0.94, decimal=2)
     assert_almost_equal(evsw, 0.94, decimal=2)
diff --git a/sklearn/model_selection/tests/test_search.py b/sklearn/model_selection/tests/test_search.py
index 1673040f96bc6..75b48a52423a5 100644
--- a/sklearn/model_selection/tests/test_search.py
+++ b/sklearn/model_selection/tests/test_search.py
@@ -13,7 +13,7 @@
 import scipy.sparse as sp
 import pytest
 
-from sklearn.utils.fixes import sp_version
+from sklearn.utils.fixes import sp_version, parse_version
 from sklearn.utils._testing import assert_raises
 from sklearn.utils._testing import assert_warns
 from sklearn.utils._testing import assert_warns_message
@@ -813,7 +813,7 @@ def test_param_sampler():
                                n_iter=3, random_state=0)
     assert [x for x in sampler] == [x for x in sampler]
 
-    if sp_version >= (0, 16):
+    if sp_version >= parse_version("0.16"):
         param_distributions = {"C": uniform(0, 1)}
         sampler = ParameterSampler(param_distributions=param_distributions,
                                    n_iter=10, random_state=0)
diff --git a/sklearn/neighbors/_base.py b/sklearn/neighbors/_base.py
index a1eebdcf78648..d51aac550fdea 100644
--- a/sklearn/neighbors/_base.py
+++ b/sklearn/neighbors/_base.py
@@ -7,7 +7,6 @@
 #
 # License: BSD 3 clause (C) INRIA, University of Amsterdam
 from functools import partial
-from distutils.version import LooseVersion
 
 import warnings
 from abc import ABCMeta, abstractmethod
@@ -28,6 +27,7 @@
 from ..utils.multiclass import check_classification_targets
 from ..utils.validation import check_is_fitted
 from ..utils.validation import check_non_negative
+from ..utils.fixes import parse_version
 from ..exceptions import DataConversionWarning, EfficiencyWarning
 
 VALID_METRICS = dict(ball_tree=BallTree.valid_metrics,
@@ -649,7 +649,7 @@ class from an array representing our data set and ask who's
                     "%s does not work with sparse matrices. Densify the data, "
                     "or set algorithm='brute'" % self._fit_method)
             old_joblib = (
-                    LooseVersion(joblib.__version__) < LooseVersion('0.12'))
+                    parse_version(joblib.__version__) < parse_version('0.12'))
             if old_joblib:
                 # Deal with change of API in joblib
                 check_pickle = False if old_joblib else None
@@ -956,7 +956,7 @@ class from an array representing our data set and ask who's
                     "or set algorithm='brute'" % self._fit_method)
 
             n_jobs = effective_n_jobs(self.n_jobs)
-            if LooseVersion(joblib.__version__) < LooseVersion('0.12'):
+            if parse_version(joblib.__version__) < parse_version('0.12'):
                 # Deal with change of API in joblib
                 delayed_query = delayed(_tree_query_radius_parallel_helper,
                                         check_pickle=False)
diff --git a/sklearn/neighbors/tests/test_dist_metrics.py b/sklearn/neighbors/tests/test_dist_metrics.py
index 7a2ec2ce4cd73..60a4e97880e15 100644
--- a/sklearn/neighbors/tests/test_dist_metrics.py
+++ b/sklearn/neighbors/tests/test_dist_metrics.py
@@ -6,13 +6,12 @@
 
 import pytest
 
-from distutils.version import LooseVersion
-from scipy import __version__ as scipy_version
 from scipy.spatial.distance import cdist
 from sklearn.neighbors import DistanceMetric
 from sklearn.neighbors import BallTree
 from sklearn.utils import check_random_state
 from sklearn.utils._testing import assert_raises_regex
+from sklearn.utils.fixes import sp_version, parse_version
 
 
 def dist_func(x1, x2, p):
@@ -106,7 +105,7 @@ def check_pdist_bool(metric, D_true):
     # Based on https://github.com/scipy/scipy/pull/7373
     # When comparing two all-zero vectors, scipy>=1.2.0 jaccard metric
     # was changed to return 0, instead of nan.
-    if metric == 'jaccard' and LooseVersion(scipy_version) < '1.2.0':
+    if metric == 'jaccard' and sp_version < parse_version('1.2.0'):
         D_true[np.isnan(D_true)] = 0
     assert_array_almost_equal(D12, D_true)
 
diff --git a/sklearn/preprocessing/_data.py b/sklearn/preprocessing/_data.py
index cc8776951f114..a435bfaf208fe 100644
--- a/sklearn/preprocessing/_data.py
+++ b/sklearn/preprocessing/_data.py
@@ -2642,7 +2642,7 @@ def quantile_transform(X, *, axis=0, n_quantiles=1000,
         input is already a numpy array). If True, a copy of `X` is transformed,
         leaving the original `X` unchanged
 
-        ..versionchnanged:: 0.23
+        ..versionchanged:: 0.23
             The default value of `copy` changed from False to True in 0.23.
 
     Returns
diff --git a/sklearn/tests/test_pipeline.py b/sklearn/tests/test_pipeline.py
index b6a3db5329add..9d5281f40e569 100644
--- a/sklearn/tests/test_pipeline.py
+++ b/sklearn/tests/test_pipeline.py
@@ -1,7 +1,6 @@
 """
 Test the pipeline module.
 """
-from distutils.version import LooseVersion
 from tempfile import mkdtemp
 import shutil
 import time
@@ -20,6 +19,7 @@
 from sklearn.utils._testing import assert_array_equal
 from sklearn.utils._testing import assert_array_almost_equal
 from sklearn.utils._testing import assert_no_warnings
+from sklearn.utils.fixes import parse_version
 
 from sklearn.base import clone, BaseEstimator, TransformerMixin
 from sklearn.pipeline import Pipeline, FeatureUnion, make_pipeline, make_union
@@ -1027,7 +1027,7 @@ def test_pipeline_memory():
     y = iris.target
     cachedir = mkdtemp()
     try:
-        if LooseVersion(joblib.__version__) < LooseVersion('0.12'):
+        if parse_version(joblib.__version__) < parse_version('0.12'):
             # Deal with change of API in joblib
             memory = joblib.Memory(cachedir=cachedir, verbose=10)
         else:
@@ -1089,7 +1089,7 @@ def test_pipeline_memory():
 
 def test_make_pipeline_memory():
     cachedir = mkdtemp()
-    if LooseVersion(joblib.__version__) < LooseVersion('0.12'):
+    if parse_version(joblib.__version__) < parse_version('0.12'):
         # Deal with change of API in joblib
         memory = joblib.Memory(cachedir=cachedir, verbose=10)
     else:
diff --git a/sklearn/utils/__init__.py b/sklearn/utils/__init__.py
index a14f1b4a596ac..67ff91f5bc95c 100644
--- a/sklearn/utils/__init__.py
+++ b/sklearn/utils/__init__.py
@@ -24,7 +24,7 @@
 from . import _joblib
 from ..exceptions import DataConversionWarning
 from .deprecation import deprecated
-from .fixes import np_version
+from .fixes import np_version, parse_version
 from ._estimator_html_repr import estimator_html_repr
 from .validation import (as_float_array,
                          assert_all_finite,
@@ -171,7 +171,7 @@ def axis0_safe_slice(X, mask, len_mask):
 
 def _array_indexing(array, key, key_dtype, axis):
     """Index an array or scipy.sparse consistently across NumPy version."""
-    if np_version < (1, 12) or issparse(array):
+    if np_version < parse_version('1.12') or issparse(array):
         # FIXME: Remove the check for NumPy when using >= 1.12
         # check if we have an boolean array-likes to make the proper indexing
         if key_dtype == 'bool':
diff --git a/sklearn/utils/estimator_checks.py b/sklearn/utils/estimator_checks.py
index ec28cb22919f0..d1c45f9fe4f75 100644
--- a/sklearn/utils/estimator_checks.py
+++ b/sklearn/utils/estimator_checks.py
@@ -719,15 +719,12 @@ def check_estimator_sparse_data(name, estimator_orig):
     X[X < .8] = 0
     X = _pairwise_estimator_convert_X(X, estimator_orig)
     X_csr = sparse.csr_matrix(X)
-    tags = estimator_orig._get_tags()
-    if tags['binary_only']:
-        y = (2 * rng.rand(40)).astype(np.int)
-    else:
-        y = (4 * rng.rand(40)).astype(np.int)
+    y = (4 * rng.rand(40)).astype(int)
     # catch deprecation warnings
     with ignore_warnings(category=FutureWarning):
         estimator = clone(estimator_orig)
     y = _enforce_estimator_tags_y(estimator, y)
+    tags = estimator_orig._get_tags()
     for matrix_format, X in _generate_sparse_matrix(X_csr):
         # catch deprecation warnings
         with ignore_warnings(category=FutureWarning):
@@ -825,10 +822,7 @@ def check_sample_weights_list(name, estimator_orig):
         n_samples = 30
         X = _pairwise_estimator_convert_X(rnd.uniform(size=(n_samples, 3)),
                                           estimator_orig)
-        if estimator._get_tags()['binary_only']:
-            y = np.arange(n_samples) % 2
-        else:
-            y = np.arange(n_samples) % 3
+        y = np.arange(n_samples) % 3
         y = _enforce_estimator_tags_y(estimator, y)
         sample_weight = [3] * n_samples
         # Test that estimators don't raise any exception
@@ -905,10 +899,7 @@ def check_dtype_object(name, estimator_orig):
     X = _pairwise_estimator_convert_X(rng.rand(40, 10), estimator_orig)
     X = X.astype(object)
     tags = estimator_orig._get_tags()
-    if tags['binary_only']:
-        y = (X[:, 0] * 2).astype(np.int)
-    else:
-        y = (X[:, 0] * 4).astype(np.int)
+    y = (X[:, 0] * 4).astype(int)
     estimator = clone(estimator_orig)
     y = _enforce_estimator_tags_y(estimator, y)
 
@@ -1007,9 +998,7 @@ def check_dont_overwrite_parameters(name, estimator_orig):
     rnd = np.random.RandomState(0)
     X = 3 * rnd.uniform(size=(20, 3))
     X = _pairwise_estimator_convert_X(X, estimator_orig)
-    y = X[:, 0].astype(np.int)
-    if estimator._get_tags()['binary_only']:
-        y[y == 2] = 1
+    y = X[:, 0].astype(int)
     y = _enforce_estimator_tags_y(estimator, y)
 
     if hasattr(estimator, "n_components"):
@@ -1060,8 +1049,6 @@ def check_fit2d_predict1d(name, estimator_orig):
     X = _pairwise_estimator_convert_X(X, estimator_orig)
     y = X[:, 0].astype(np.int)
     tags = estimator_orig._get_tags()
-    if tags['binary_only']:
-        y[y == 2] = 1
     estimator = clone(estimator_orig)
     y = _enforce_estimator_tags_y(estimator, y)
 
@@ -1109,9 +1096,7 @@ def check_methods_subset_invariance(name, estimator_orig):
     rnd = np.random.RandomState(0)
     X = 3 * rnd.uniform(size=(20, 3))
     X = _pairwise_estimator_convert_X(X, estimator_orig)
-    y = X[:, 0].astype(np.int)
-    if estimator_orig._get_tags()['binary_only']:
-        y[y == 2] = 1
+    y = X[:, 0].astype(int)
     estimator = clone(estimator_orig)
     y = _enforce_estimator_tags_y(estimator, y)
 
@@ -1383,10 +1368,7 @@ def check_fit_score_takes_y(name, estimator_orig):
     n_samples = 30
     X = rnd.uniform(size=(n_samples, 3))
     X = _pairwise_estimator_convert_X(X, estimator_orig)
-    if estimator_orig._get_tags()['binary_only']:
-        y = np.arange(n_samples) % 2
-    else:
-        y = np.arange(n_samples) % 3
+    y = np.arange(n_samples) % 3
     estimator = clone(estimator_orig)
     y = _enforce_estimator_tags_y(estimator, y)
     set_random_state(estimator)
@@ -1416,8 +1398,6 @@ def check_estimators_dtypes(name, estimator_orig):
     X_train_int_64 = X_train_32.astype(np.int64)
     X_train_int_32 = X_train_32.astype(np.int32)
     y = X_train_int_64[:, 0]
-    if estimator_orig._get_tags()['binary_only']:
-        y[y == 2] = 1
     y = _enforce_estimator_tags_y(estimator_orig, y)
 
     methods = ["predict", "transform", "decision_function", "predict_proba"]
@@ -1596,6 +1576,7 @@ def check_estimators_partial_fit_n_features(name, estimator_orig):
     estimator = clone(estimator_orig)
     X, y = make_blobs(n_samples=50, random_state=1)
     X -= X.min()
+    y = _enforce_estimator_tags_y(estimator_orig, y)
 
     try:
         if is_classifier(estimator):
@@ -2062,11 +2043,7 @@ def check_classifiers_multilabel_representation_invariance(name,
 def check_estimators_fit_returns_self(name, estimator_orig,
                                       readonly_memmap=False):
     """Check if self is returned when calling fit"""
-    if estimator_orig._get_tags()['binary_only']:
-        n_centers = 2
-    else:
-        n_centers = 3
-    X, y = make_blobs(random_state=0, n_samples=21, centers=n_centers)
+    X, y = make_blobs(random_state=0, n_samples=21)
     # some want non-negative input
     X -= X.min()
     X = _pairwise_estimator_convert_X(X, estimator_orig)
@@ -2108,10 +2085,7 @@ def check_supervised_y_2d(name, estimator_orig):
     X = _pairwise_estimator_convert_X(
         rnd.uniform(size=(n_samples, 3)), estimator_orig
     )
-    if tags['binary_only']:
-        y = np.arange(n_samples) % 2
-    else:
-        y = np.arange(n_samples) % 3
+    y = np.arange(n_samples) % 3
     y = _enforce_estimator_tags_y(estimator_orig, y)
     estimator = clone(estimator_orig)
     set_random_state(estimator)
@@ -2436,11 +2410,7 @@ def check_class_weight_balanced_linear_classifier(name, Classifier):
 
 @ignore_warnings(category=FutureWarning)
 def check_estimators_overwrite_params(name, estimator_orig):
-    if estimator_orig._get_tags()['binary_only']:
-        n_centers = 2
-    else:
-        n_centers = 3
-    X, y = make_blobs(random_state=0, n_samples=21, centers=n_centers)
+    X, y = make_blobs(random_state=0, n_samples=21)
     # some want non-negative input
     X -= X.min()
     X = _pairwise_estimator_convert_X(X, estimator_orig, kernel=rbf_kernel)
@@ -2511,7 +2481,8 @@ def check_no_attributes_set_in_init(name, estimator_orig):
 def check_sparsify_coefficients(name, estimator_orig):
     X = np.array([[-2, -1], [-1, -1], [-1, -2], [1, 1], [1, 2], [2, 1],
                   [-1, -2], [2, 2], [-2, -2]])
-    y = [1, 1, 1, 2, 2, 2, 3, 3, 3]
+    y = np.array([1, 1, 1, 2, 2, 2, 3, 3, 3])
+    y = _enforce_estimator_tags_y(estimator_orig, y)
     est = clone(estimator_orig)
 
     est.fit(X, y)
@@ -2535,7 +2506,7 @@ def check_classifier_data_not_an_array(name, estimator_orig):
     X = np.array([[3, 0], [0, 1], [0, 2], [1, 1], [1, 2], [2, 1],
                   [0, 3], [1, 0], [2, 0], [4, 4], [2, 3], [3, 2]])
     X = _pairwise_estimator_convert_X(X, estimator_orig)
-    y = [1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2]
+    y = np.array([1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2])
     y = _enforce_estimator_tags_y(estimator_orig, y)
     for obj_type in ["NotAnArray", "PandasDataframe"]:
         check_estimators_data_not_an_array(name, estimator_orig, X, y,
@@ -2682,6 +2653,9 @@ def _enforce_estimator_tags_y(estimator, y):
         # Create strictly positive y. The minimal increment above 0 is 1, as
         # y could be of integer dtype.
         y += 1 + abs(y.min())
+    # Estimators with a `binary_only` tag only accept up to two unique y values
+    if estimator._get_tags()["binary_only"] and y.size > 0:
+        y = np.where(y == y.flat[0], y, y.flat[0] + 1)
     # Estimators in mono_output_task_error raise ValueError if y is of 1-D
     # Convert into a 2-D y for those estimators.
     if estimator._get_tags()["multioutput_only"]:
diff --git a/sklearn/utils/fixes.py b/sklearn/utils/fixes.py
index f912df19c323b..bc820cba364d0 100644
--- a/sklearn/utils/fixes.py
+++ b/sklearn/utils/fixes.py
@@ -21,27 +21,22 @@
 
 from .deprecation import deprecated
 
+try:
+    from pkg_resources import parse_version  # type: ignore
+except ImportError:
+    # setuptools not installed
+    parse_version = LooseVersion  # type: ignore
 
-def _parse_version(version_string):
-    version = []
-    for x in version_string.split('.'):
-        try:
-            version.append(int(x))
-        except ValueError:
-            # x may be of the form dev-1ea1592
-            version.append(x)
-    return tuple(version)
 
+np_version = parse_version(np.__version__)
+sp_version = parse_version(scipy.__version__)
 
-np_version = _parse_version(np.__version__)
-sp_version = _parse_version(scipy.__version__)
 
-
-if sp_version >= (1, 4):
+if sp_version >= parse_version('1.4'):
     from scipy.sparse.linalg import lobpcg
 else:
     # Backport of lobpcg functionality from scipy 1.4.0, can be removed
-    # once support for sp_version < (1, 4) is dropped
+    # once support for sp_version < parse_version('1.4') is dropped
     # mypy error: Name 'lobpcg' already defined (possibly by an import)
     from ..externals._lobpcg import lobpcg  # type: ignore  # noqa
 
@@ -56,7 +51,7 @@ def _astype_copy_false(X):
     {ndarray, csr_matrix, csc_matrix}.astype when possible,
     otherwise don't specify
     """
-    if sp_version >= (1, 1) or not sp.issparse(X):
+    if sp_version >= parse_version('1.1') or not sp.issparse(X):
         return {'copy': False}
     else:
         return {}
@@ -85,7 +80,7 @@ def _joblib_parallel_args(**kwargs):
     """
     import joblib
 
-    if joblib.__version__ >= LooseVersion('0.12'):
+    if parse_version(joblib.__version__) >= parse_version('0.12'):
         return kwargs
 
     extra_args = set(kwargs.keys()).difference({'prefer', 'require'})
diff --git a/sklearn/utils/tests/test_estimator_checks.py b/sklearn/utils/tests/test_estimator_checks.py
index ade575f6f3206..71e7f2e0f5f71 100644
--- a/sklearn/utils/tests/test_estimator_checks.py
+++ b/sklearn/utils/tests/test_estimator_checks.py
@@ -25,7 +25,7 @@
 from sklearn.utils.estimator_checks import check_regressor_data_not_an_array
 from sklearn.utils.validation import check_is_fitted
 from sklearn.utils.estimator_checks import check_outlier_corruption
-from sklearn.utils.fixes import _parse_version
+from sklearn.utils.fixes import np_version, parse_version
 from sklearn.ensemble import RandomForestClassifier
 from sklearn.linear_model import LinearRegression, SGDClassifier
 from sklearn.mixture import GaussianMixture
@@ -34,7 +34,6 @@
 from sklearn.linear_model import MultiTaskElasticNet, LogisticRegression
 from sklearn.svm import SVC
 from sklearn.neighbors import KNeighborsRegressor
-from sklearn.tree import DecisionTreeClassifier
 from sklearn.utils.validation import check_array
 from sklearn.utils import all_estimators
 
@@ -306,11 +305,19 @@ def predict(self, X):
         return np.array([self.value_] * X.shape[0])
 
 
-class UntaggedBinaryClassifier(DecisionTreeClassifier):
+class UntaggedBinaryClassifier(SGDClassifier):
     # Toy classifier that only supports binary classification, will fail tests.
-    def fit(self, X, y, sample_weight=None):
-        super().fit(X, y, sample_weight)
-        if np.all(self.n_classes_ > 2):
+    def fit(self, X, y, coef_init=None, intercept_init=None,
+            sample_weight=None):
+        super().fit(X, y, coef_init, intercept_init, sample_weight)
+        if len(self.classes_) > 2:
+            raise ValueError('Only 2 classes are supported')
+        return self
+
+    def partial_fit(self, X, y, classes=None, sample_weight=None):
+        super().partial_fit(X=X, y=y, classes=classes,
+                            sample_weight=sample_weight)
+        if len(self.classes_) > 2:
             raise ValueError('Only 2 classes are supported')
         return self
 
@@ -334,8 +341,7 @@ def _more_tags(self):
 
 
 def test_not_an_array_array_function():
-    np_version = _parse_version(np.__version__)
-    if np_version < (1, 17):
+    if np_version < parse_version('1.17'):
         raise SkipTest("array_function protocol not supported in numpy <1.17")
     not_array = _NotAnArray(np.ones(10))
     msg = "Don't want to call array_function sum!"
diff --git a/sklearn/utils/tests/test_multiclass.py b/sklearn/utils/tests/test_multiclass.py
index 4cb7368d34c03..a627f8f55cc0b 100644
--- a/sklearn/utils/tests/test_multiclass.py
+++ b/sklearn/utils/tests/test_multiclass.py
@@ -3,7 +3,6 @@
 import scipy.sparse as sp
 from itertools import product
 import pytest
-from distutils.version import LooseVersion
 
 from scipy.sparse import issparse
 from scipy.sparse import csc_matrix
@@ -16,6 +15,7 @@
 from sklearn.utils._testing import assert_array_almost_equal
 from sklearn.utils._testing import assert_allclose
 from sklearn.utils.estimator_checks import _NotAnArray
+from sklearn.utils.fixes import parse_version
 
 from sklearn.utils.multiclass import unique_labels
 from sklearn.utils.multiclass import is_multilabel
@@ -307,7 +307,7 @@ def test_type_of_target():
 def test_type_of_target_pandas_sparse():
     pd = pytest.importorskip("pandas")
 
-    if LooseVersion(pd.__version__) >= '0.25':
+    if parse_version(pd.__version__) >= parse_version('0.25'):
         pd_sparse_array = pd.arrays.SparseArray
     else:
         pd_sparse_array = pd.SparseArray
diff --git a/sklearn/utils/validation.py b/sklearn/utils/validation.py
index 7a6ef1e05fdde..ad19429c1eb33 100644
--- a/sklearn/utils/validation.py
+++ b/sklearn/utils/validation.py
@@ -15,7 +15,6 @@
 
 import numpy as np
 import scipy.sparse as sp
-from distutils.version import LooseVersion
 from inspect import signature, isclass, Parameter
 
 from numpy.core.numeric import ComplexWarning
@@ -23,7 +22,7 @@
 
 from contextlib import suppress
 
-from .fixes import _object_dtype_isnan
+from .fixes import _object_dtype_isnan, parse_version
 from .. import get_config as _get_config
 from ..exceptions import NonBLASDotWarning, PositiveSpectrumWarning
 from ..exceptions import NotFittedError
@@ -228,7 +227,7 @@ def check_memory(memory):
     """
 
     if memory is None or isinstance(memory, str):
-        if LooseVersion(joblib.__version__) < '0.12':
+        if parse_version(joblib.__version__) < parse_version('0.12'):
             memory = joblib.Memory(cachedir=memory, verbose=0)
         else:
             memory = joblib.Memory(location=memory, verbose=0)