code: def _doubly_stochastic(P, tol=1e-3):
      # cleaner implementation of btaba/sinkhorn_knopp
  
      max_iter = 1000
      c = 1 / P.sum(axis=0)
      r = 1 / (P @ c)
      P_eps = P
  
      for it in range(max_iter):
          if (np.abs(P_eps.sum(axis=1) - 1) < tol).all() and (
              np.abs(P_eps.sum(axis=0) - 1) < tol
          ).all():
              # All column/row sums ~= 1 within threshold
              break
  
          c = 1 / (r @ P)
          r = 1 / (P @ c)
          P_eps = r[:, None] * P * c
  
      return P_eps
  
imports: import pytest
  import numpy
  from scipy.optimize import quadratic_assignment, OptimizeWarning
  from scipy.optimize._qap import _calc_score
  from numpy.testing import assert_equal, assert_, assert_warns
label_desc: Use image node_id 3 for calling a global function with example usage: _doubly_stochastic(P, tol) and returns: P_eps