docs: def MDS(D, dim, method='simple', theta=False):
      """ recover points from euclidean distance matrix using classic MDS algorithm. 
      """
      N = D.shape[0]
      if method == 'simple':
          d1 = D[0, :]
          # buf_ = d1 * np.ones([1, N]).T + (np.ones([N, 1]) * d1).T
          buf_ = np.broadcast_to(d1, D.shape) + np.broadcast_to(d1[:, np.newaxis], D.shape)
          np.subtract(D, buf_, out=buf_)
          G = buf_ # G = (D - d1 * np.ones([1, N]).T - (np.ones([N, 1]) * d1).T)
      elif method == 'advanced':
          # s1T = np.vstack([np.ones([1, N]), np.zeros([N - 1, N])])
          s1T = np.zeros_like(D)
          s1T[0, :] = 1
          np.subtract(np.identity(N), s1T, out = s1T)
          G = np.dot(np.dot(s1T.T, D), s1T)
      elif method == 'geometric':
          J = np.identity(N) + np.full((N, N), -1.0 / float(N))
          G = np.dot(np.dot(J, D), J)
      else:
          print('Unknown method {} in MDS'.format(method))
      G *= -0.5
      factor, u = eigendecomp(G, dim)
      if (theta):
          return theta_from_eigendecomp(factor, u)
      else:
          return x_from_eigendecomp(factor, u, dim)