def METHOD_NAME(G, fill=10, ndim=3):
    # https://en.wikipedia.org/wiki/Multidimensional_scaling
    n = len(G)
    distmat = nx.adjacency_matrix(G, weight='distance').A
    distmat[distmat == 0] = fill
    D2 = distmat ** 2
    J = np.eye(n) - 1/n * np.ones(n) * np.ones(n).T
    # Must be matmul
    B = - 0.5 * J @ D2 @ J
    S = np.linalg.eig(B)
    order = np.argsort(S[0])[::-1]
    vals = S[0][order]
    vecs = S[1][:, order]
    M = np.real(vecs[:, :ndim])
    L = np.diag(np.sqrt(vals[:ndim]))
    return M @ L