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