docs: def reconstruct_dwmds(edm, X0, W=None, n=None, r=None, X_bar=None, print_out=False, tol=1e-10, sweeps=100):
      """ Reconstruct point set using d(istributed)w(eighted) MDS.
  
      Refer to paper "Distributed Weighted-Multidimensional Scaling for Node Localization in Sensor Networks" for 
      implementation details (doi.org/10.1145/1138127.1138129)
  
      :param X0: Nxd matrix of starting points.
      :param n: Number of points of unknown position. The first n points in X0 and edm are considered unknown. 
      :param tol: Stopping criterion: when the cost is below this level, we stop. 
      :param sweeps: Maximum number of sweeps. 
      """
      from .basics import get_edm
      from .distributed_mds import get_b, get_Si
  
      N, d = X0.shape
      if r is None and n is None:
          raise ValueError('either r or n have to be given.')
      elif n is None:
          n = r.shape[0]
  
      if W is None:
          W = np.ones((N, N)) - np.eye(N)
  
      X_k = X0.copy()
  
      costs = []
      # don't have to ignore i=j, because W[i,i] is zero.
      a = np.sum(W[:n, :n], axis=1).flatten() + 2 * \
          np.sum(W[:n, n:], axis=1).flatten()
      if r is not None:
          a += r.flatten()
      for k in range(sweeps):
          S = 0
          for i in range(n):
              edm_estimated = get_edm(X_k)
              bi = get_b(i, edm_estimated, W, edm, n)
              if r is not None and X_bar is not None:
                  X_k[i] = 1 / a[i] * (r[i] * X_bar[i, :] +
                                       X_k.T.dot(bi).flatten())
                  Si = get_Si(i, edm_estimated, edm, W, n, r, X_bar[i], X_k[i])
              else:
                  X_k[i] = 1 / a[i] * X_k.T.dot(bi).flatten()
                  Si = get_Si(i, edm_estimated, edm, W, n)
              S += Si
          costs.append(S)
          if k > 1 and abs(costs[-1] - costs[-2]) < tol:
              if (print_out):
                  print('dwMDS: converged after', k)
              break
      return X_k, costs