docs: def isomap(geom, n_components=8, eigen_solver='auto',
             random_state=None, path_method='auto',
             distance_matrix=None, graph_distance_matrix = None,
             centered_matrix=None, solver_kwds=None):
      """
      Parameters
      ----------
      geom : a Geometry object from megaman.geometry.geometry
      n_components : integer, optional
          The dimension of the projection subspace.
      eigen_solver : {'auto', 'dense', 'arpack', 'lobpcg', or 'amg'}
          'auto' :
              algorithm will attempt to choose the best method for input data
          'dense' :
              use standard dense matrix operations for the eigenvalue decomposition.
              For this method, M must be an array or matrix type.  This method should be avoided for large problems.
          'arpack' :
              use arnoldi iteration in shift-invert mode. For this method,
              M may be a dense matrix, sparse matrix, or general linear operator.
              Warning: ARPACK can be unstable for some problems.  It is best to
              try several random seeds in order to check results.
          'lobpcg' :
              Locally Optimal Block Preconditioned Conjugate Gradient Method.
              A preconditioned eigensolver for large symmetric positive definite
              (SPD) generalized eigenproblems.
          'amg' :
              AMG requires pyamg to be installed. It can be faster on very large,
              sparse problems, but may also lead to instabilities.
      random_state : int seed, RandomState instance, or None (default)
          A pseudo random number generator used for the initialization of the
          lobpcg eigen vectors decomposition when eigen_solver == 'amg'.
          By default, arpack is used.
      path_method : string, method for computing graph shortest path. One of :
          'auto', 'D', 'FW', 'BF', 'J'. See scipy.sparse.csgraph.shortest_path
          for more information.
      distance_matrix : sparse Ndarray (n_obs, n_obs), optional. Pairwise distance matrix
          sparse zeros considered 'infinite'.
      graph_distance_matrix : Ndarray (n_obs, n_obs), optional. Pairwise graph distance
          matrix. Output of graph_shortest_path.
      centered_matrix : Ndarray (n_obs, n_obs), optional. Centered version of
          graph_distance_matrix
      solver_kwds : any additional keyword arguments to pass to the selected eigen_solver
  
      Returns
      -------
      embedding : array, shape=(n_samples, n_components)
          The reduced samples.
  
      Notes
      -----
      """
      # Step 1: use geometry to calculate the distance matrix
      if ((distance_matrix is None) and (centered_matrix is None)):
          if geom.adjacency_matrix is None:
              distance_matrix = geom.compute_adjacency_matrix()
          else:
              distance_matrix = geom.adjacency_matrix
  
      # Step 2: use graph_shortest_path to construct D_G
      ## WARNING: D_G is an (NxN) DENSE matrix!!
      if ((graph_distance_matrix is None) and (centered_matrix is None)):
          graph_distance_matrix = graph_shortest_path(distance_matrix,
                                                      method=path_method,
                                                      directed=False)
  
      # Step 3: center graph distance matrix
      if centered_matrix is None:
          centered_matrix = center_matrix(graph_distance_matrix)
  
  
      # Step 4: compute d largest eigenvectors/values of centered_matrix
      lambdas, diffusion_map = eigen_decomposition(centered_matrix, n_components,
                                                   largest=True,
                                                   eigen_solver=eigen_solver,
                                                   random_state=random_state,
                                                   solver_kwds=solver_kwds)
      # Step 5:
      # return Y = [sqrt(lambda_1)*V_1, ..., sqrt(lambda_d)*V_d]
      ind = np.argsort(lambdas); ind = ind[::-1] # sort largest
      lambdas = lambdas[ind];
      diffusion_map = diffusion_map[:, ind]
      embedding = diffusion_map[:, 0:n_components] * np.sqrt(lambdas[0:n_components])
      return embedding