Source code for pycgsp.graph

# #############################################################################
# __init__.py
# ===========
# Author : Matthieu Simeoni [matthieu.simeoni@gmail.com]
# #############################################################################

r"""
Routines for building graphs from tessellations/point sets in :math:`\mathbb{R}^3`.
"""

import numpy as np
import scipy.linalg as linalg
import scipy.sparse as sp
import scipy.sparse.linalg as splinalg
import scipy.spatial as spatial
from pygsp.graphs import Graph
import healpy as hp
from typing import Tuple


[docs]def cvxhull_graph(R: np.ndarray, cheb_normalized: bool = True, compute_differential_operator: bool = True) \ -> Tuple[Graph, float]: r""" Build the convex hull graph of a point set in :math:`\mathbb{R}^3`. The graph edges have exponential-decay weighting. Definitions of the graph Laplacians: .. math:: L = I - D^{-1/2} W D^{-1/2},\qquad L_{n} = (2 / \mu_{\max}) L - I Parameters ---------- R : :py:class:`~numpy.ndarray` (N,3) Cartesian coordinates of point set with size N. All points must be **distinct**. cheb_normalized : bool Rescale Laplacian spectrum to [-1, 1]. compute_differential_operator : bool Computes the graph gradient. Returns ------- G : :py:class:`~pygsp.graphs.Graph` If ``cheb_normalized = True``, ``G.Ln`` is created (Chebyshev Laplacian :math:`L_{n}` above) If ``compute_differential_operator = True``, ``G.D`` is created and contains the gradient. rho : float Scale parameter :math:`\rho` corresponding to the average distance of a point on the graph to its nearest neighbors. Examples -------- .. plot:: import numpy as np from pycgsp.graph import cvxhull_graph from pygsp.plotting import plot_graph theta, phi = np.linspace(0,np.pi,6, endpoint=False)[1:], np.linspace(0,2*np.pi,9, endpoint=False) theta, phi = np.meshgrid(theta, phi) x,y,z = np.cos(phi)*np.sin(theta), np.sin(phi)*np.sin(theta), np.cos(theta) R = np.stack((x.flatten(), y.flatten(), z.flatten()), axis=-1) G, _ = cvxhull_graph(R) plot_graph(G) Warnings -------- In the newest version of PyGSP (> 0.5.1) the convention is changed: ``Graph.D`` is the divergence operator and ``Graph.D.transpose()`` the gradient (see routine `Graph.compute_differential_operator <https://pygsp.readthedocs.io/en/latest/reference/graphs.html#pygsp.graphs.Graph.compute_differential_operator>`_). The code should be adapted when this new version is released. """ # Form convex hull to extract nearest neighbors. Each row in # cvx_hull.simplices is a triangle of connected points. cvx_hull = spatial.ConvexHull(R) cols = np.roll(cvx_hull.simplices, shift=1, axis=-1).reshape(-1) rows = cvx_hull.simplices.reshape(-1) # Form sparse affinity matrix from extracted pairs W = sp.coo_matrix((cols * 0 + 1, (rows, cols)), shape=(cvx_hull.vertices.size, cvx_hull.vertices.size)) # Symmetrize the matrix to obtain an undirected graph. extended_row = np.concatenate([W.row, W.col]) extended_col = np.concatenate([W.col, W.row]) W.row, W.col = extended_row, extended_col W.data = np.concatenate([W.data, W.data]) W = W.tocsr().tocoo() # Delete potential duplicate pairs # Weight matrix elements according to the exponential kernel distance = linalg.norm(cvx_hull.points[W.row, :] - cvx_hull.points[W.col, :], axis=-1) rho = np.mean(distance) W.data = np.exp(- (distance / rho) ** 2) W = W.tocsc() G = _graph_laplacian(W, R, compute_differential_operator=compute_differential_operator, cheb_normalized=cheb_normalized) return G, rho
def _graph_laplacian(W, R, compute_differential_operator=False, cheb_normalized=False): # Form Graph Laplacian G = Graph(W, gtype='undirected', lap_type='normalized', coords=R) G.compute_laplacian(lap_type='normalized') # Stored in G.L, sparse matrix, csc ordering if compute_differential_operator is True: G.compute_differential_operator() # stored in G.D, also accessible via G.grad() or G.div() (for the adjoint). else: pass if cheb_normalized: D_max = splinalg.eigsh(G.L, k=1, return_eigenvectors=False) Ln = (2 / D_max[0]) * G.L - sp.identity(W.shape[0], dtype=np.float, format='csc') G.Ln = Ln else: pass return G
[docs]def healpix_nngraph(nside: int, cheb_normalized: bool = True, compute_differential_operator: bool = True) \ -> Tuple[Graph, float]: r""" Build the nearest neighbour graph of a HEALPix spherical point set. The graph edges have exponential-decay weighting. Definitions of the graph Laplacians: .. math:: L = I - D^{-1/2} W D^{-1/2},\qquad L_{n} = (2 / \mu_{\max}) L - I Parameters ---------- nside: int Parameter NSIDE of the `HEALPix discretisation scheme <https://healpix.jpl.nasa.gov/>`_. cheb_normalized : bool Rescale Laplacian spectrum to [-1, 1]. compute_differential_operator : bool Computes the graph gradient. Returns ------- G : :py:class:`~pygsp.graphs.Graph` If ``cheb_normalized = True``, ``G.Ln`` is created (Chebyshev Laplacian :math:`L_{n}` above) If ``compute_differential_operator = True``, ``G.D`` is created and contains the gradient. rho : float Scale parameter :math:`\rho` corresponding to the average distance of a point on the graph to its nearest neighbors. Examples -------- .. plot:: from pycgsp.graph import healpix_nngraph from pygsp.plotting import plot_graph G, _ = healpix_nngraph(nside=2) plot_graph(G) Warnings -------- In the newest version of PyGSP (> 0.5.1) the convention is changed: ``Graph.D`` is the divergence operator and ``Graph.D.transpose()`` the gradient (see routine `Graph.compute_differential_operator <https://pygsp.readthedocs.io/en/latest/reference/graphs.html#pygsp.graphs.Graph.compute_differential_operator>`_). The code should be adapted when this new version is released. """ npix = hp.nside2npix(nside) x, y, z = hp.pix2vec(nside, np.arange(npix)) R = np.stack((x, y, z), axis=-1) cols = hp.get_all_neighbours(nside, np.arange(npix)).transpose().reshape(-1) cols[cols == -1] = npix - 1 rows = np.repeat(np.arange(npix), 8, axis=-1).transpose().reshape(-1) # Form sparse affinity matrix from extracted pairs W = sp.coo_matrix((cols * 0 + 1, (rows, cols)), shape=(npix, npix)) # Symmetrize the matrix to obtain an undirected graph. extended_row = np.concatenate([W.row, W.col]) extended_col = np.concatenate([W.col, W.row]) W.row, W.col = extended_row, extended_col W.data = np.concatenate([W.data, W.data]) W = W.tocsr().tocoo() # Delete potential duplicate pairs # Weight matrix elements according to the exponential kernel distance = linalg.norm(R[W.row, :] - R[W.col, :], axis=-1) rho = np.mean(distance) W.data = np.exp(- (distance / rho) ** 2) W = W.tocsc() G = _graph_laplacian(W, R, compute_differential_operator=compute_differential_operator, cheb_normalized=cheb_normalized) return G, rho