Add implementation of fast approximation of LB problems

[luisfpereira]

Inspired by #9, this PR implements the general strategy for fast approximation of Laplace-Beltrami eigenproblems proposed in (Nasikun, 2018) (TODO: check if general methodology was proposed before and this paper only brings new ways of computing distances; if yes update references/namings)

Computation of distance is pluggable (use any from #11).
I haven't implemented yet the variants of Dijsktra they suggest. Dijsktra with corrected distances is trivial and we should implement it fast. I've implemented a weird version of Dijsktra that always selects the same number of neighbors for the support of the local functions (probably useless and behaving poorly, but fun to implement nonetheless).

Haven't explored much the local transform, probably need to correct that part for considering the support.

Notice (Magnet, 2023) uses this, but with a different distance.

We can probably bring the notion of a HierarchicalMesh based on the sampling of a mesh here (SampledBasedHierarchicalMesh?).

[luisfpereira]

btw, with the proper merging of branches (#11 ad #12), you should be able to run the following code and get some initial interesting results to verify the code may be surprisingly working:

dataset = NotebooksDataset()
mesh = TriangleMesh.from_file(dataset.get_filename("cat-00"))

# a weird metric just for fun
mesh.metric = FixedNeighborsSingleSourceDijkstra(mesh, n_neighbors=400)

local_func_constr = NasikunLocalFunctionsConstructor(min_n_samples=150)

local_func_mat = local_func_constr(mesh)

# checking all the vertices belong to the support of a local function (code is still naive to do updates if not)
# chose parameters above accordingly
print(
    np.where(np.sum(local_func_mat.toarray(), axis=-1) < 0.99)[0].shape,
    np.where(np.sum(local_func_mat = local_func_constr(mesh).toarray(), axis=-1) > 1.01)[0].shape,
)

# this is done by `NasikunLaplacianSpectrumFinder`, but probably nice to see it in a scripting way
restricted_mass_matrix = local_func_mat.T @ mesh.laplacian.mass_matrix @ local_func_mat
restricted_stiffness_matrix = local_func_mat.T @ mesh.laplacian.stiffness_matrix @ local_func_mat

eig_solver = ScipyEigsh(spectrum_size=5, sigma=-0.01)
eigenvals, restricted_eigenvecs = eig_solver(
    restricted_stiffness_matrix, M=restricted_mass_matrix
)

eigenvecs = local_func_mat @ restricted_eigenvecs

And now some visualization:

pv_mesh = pv.PolyData.from_regular_faces(points=mesh.vertices, faces=mesh.faces)


pv_mesh["scalar"] = eigenvecs[:, 1]
# uncomment to see ground truth
# pv_mesh["scalar"] = mesh.basis.vecs[:, 1]


pl = pv.Plotter()

pl.add_mesh(pv_mesh)

pl.camera.roll = 0
pl.camera.azimuth = -10
pl.camera.elevation = -20
pl.show_axes()

pl.show()