diff --git a/doc/history.rst b/doc/history.rst index df18ad9b..a49f98a5 100644 --- a/doc/history.rst +++ b/doc/history.rst @@ -5,6 +5,20 @@ History 0.6.0 (xxxx-xx-xx) ------------------ +* Filters approximated by Chebyshev polynomials are now implemented as separate + filters, i.e. you need to do Chebyshev(Heat(tau=1), order=10).filter(signal) to + filter a signal with a heat kernel approximated by a Chebyshev polynomial of + order 10. The reasons are multiple: + - They are arguably different filters. + Try Heat().plot() and Chebyshev(Heat(), order=1).plot(). + - It allows to visualize the approximated filters in the spectral domain, and + compare them with their continuous counterpart or other approximations. + - One can now instantiates filters which are solely defined by their Chebyshev + coefficients. That is the case when learning them with back-propagation in a + neural network setting, and it's useful to (i) visualize the learned filters + in the spectral and vertex domains, and (ii) to compute and visualize the + feature maps. + See the new tutorial on filter approximations for usage. * print(graph) and print(filters) now show valuable information. * Building a graph object is much faster. * New rectangular filter (low-pass and band-pass). diff --git a/pygsp/filters/__init__.py b/pygsp/filters/__init__.py index 6bcc474a..9bb17bf8 100644 --- a/pygsp/filters/__init__.py +++ b/pygsp/filters/__init__.py @@ -20,6 +20,7 @@ Filter.synthesize Filter.compute_frame Filter.estimate_frame_bounds + Filter.approximate Filter.plot Filter.localize @@ -85,6 +86,10 @@ **Chebyshev polynomials** +.. autosummary:: + + Chebyshev + .. autosummary:: compute_cheby_coeff @@ -130,7 +135,9 @@ 'lanczos_op' ] -__all__ = _FILTERS + _APPROXIMATIONS +__all__ = _FILTERS + _APPROXIMATIONS + ['Chebyshev'] _utils.import_classes(_FILTERS, 'filters', 'filters') -_utils.import_functions(_APPROXIMATIONS, 'filters.approximations', 'filters') +_utils.import_functions(_APPROXIMATIONS, 'filters.approximations_old', 'filters') + +from .approximations import * diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index 905ed3e8..0f0ebb23 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -1,343 +1,362 @@ # -*- coding: utf-8 -*- +from __future__ import division + import numpy as np from scipy import sparse from pygsp import utils +from . import Filter # prevent circular import in Python < 3.5 _logger = utils.build_logger(__name__) -@utils.filterbank_handler -def compute_cheby_coeff(f, m=30, N=None, *args, **kwargs): - r""" - Compute Chebyshev coefficients for a Filterbank. - - Parameters - ---------- - f : Filter - Filterbank with at least 1 filter - m : int - Maximum order of Chebyshev coeff to compute - (default = 30) - N : int - Grid order used to compute quadrature - (default = m + 1) - i : int - Index of the Filterbank element to compute - (default = 0) - - Returns - ------- - c : ndarray - Matrix of Chebyshev coefficients - - """ - G = f.G - i = kwargs.pop('i', 0) - - if not N: - N = m + 1 - - a_arange = [0, G.lmax] - - a1 = (a_arange[1] - a_arange[0]) / 2 - a2 = (a_arange[1] + a_arange[0]) / 2 - c = np.zeros(m + 1) - - tmpN = np.arange(N) - num = np.cos(np.pi * (tmpN + 0.5) / N) - for o in range(m + 1): - c[o] = 2. / N * np.dot(f._kernels[i](a1 * num + a2), - np.cos(np.pi * o * (tmpN + 0.5) / N)) - - return c - - -def cheby_op(G, c, signal, **kwargs): - r""" - Chebyshev polynomial of graph Laplacian applied to vector. - - Parameters - ---------- - G : Graph - c : ndarray or list of ndarrays - Chebyshev coefficients for a Filter or a Filterbank - signal : ndarray - Signal to filter - - Returns - ------- - r : ndarray - Result of the filtering - - """ - # Handle if we do not have a list of filters but only a simple filter in cheby_coeff. - if not isinstance(c, np.ndarray): - c = np.array(c) - - c = np.atleast_2d(c) - Nscales, M = c.shape - - if M < 2: - raise TypeError("The coefficients have an invalid shape") - - # thanks to that, we can also have 1d signal. - try: - Nv = np.shape(signal)[1] - r = np.zeros((G.N * Nscales, Nv)) - except IndexError: - r = np.zeros((G.N * Nscales)) - - a_arange = [0, G.lmax] - - a1 = float(a_arange[1] - a_arange[0]) / 2. - a2 = float(a_arange[1] + a_arange[0]) / 2. - - twf_old = signal - twf_cur = (G.L.dot(signal) - a2 * signal) / a1 - - tmpN = np.arange(G.N, dtype=int) - for i in range(Nscales): - r[tmpN + G.N*i] = 0.5 * c[i, 0] * twf_old + c[i, 1] * twf_cur - - factor = 2/a1 * (G.L - a2 * sparse.eye(G.N)) - for k in range(2, M): - twf_new = factor.dot(twf_cur) - twf_old - for i in range(Nscales): - r[tmpN + G.N*i] += c[i, k] * twf_new - - twf_old = twf_cur - twf_cur = twf_new - - return r - - -def cheby_rect(G, bounds, signal, **kwargs): - r""" - Fast filtering using Chebyshev polynomial for a perfect rectangle filter. - - Parameters - ---------- - G : Graph - bounds : array-like - The bounds of the pass-band filter - signal : array-like - Signal to filter - order : int (optional) - Order of the Chebyshev polynomial (default: 30) - - Returns - ------- - r : array-like - Result of the filtering - - """ - if not (isinstance(bounds, (list, np.ndarray)) and len(bounds) == 2): - raise ValueError('Bounds of wrong shape.') - - bounds = np.array(bounds) - - m = int(kwargs.pop('order', 30) + 1) - - try: - Nv = np.shape(signal)[1] - r = np.zeros((G.N, Nv)) - except IndexError: - r = np.zeros((G.N)) - - b1, b2 = np.arccos(2. * bounds / G.lmax - 1.) - factor = 4./G.lmax * G.L - 2.*sparse.eye(G.N) +class Chebyshev(Filter): + r"""Approximate continuous filters with Chebyshev polynomials. - T_old = signal - T_cur = factor.dot(signal) / 2. - r = (b1 - b2)/np.pi * signal + 2./np.pi * (np.sin(b1) - np.sin(b2)) * T_cur + Math which explains the polynomial filters sum_k theta_k lambda^k + Weighted sum of diffused versions of the signal + Note recursive computation. O(N) computational cost and 4N space. - for k in range(2, m): - T_new = factor.dot(T_cur) - T_old - r += 2./(k*np.pi) * (np.sin(k*b1) - np.sin(k*b2)) * T_new - T_old = T_cur - T_cur = T_new + Math to show how the coefficients are computed - return r + Evaluation methods (which can be passed when calling :meth:`Filter.evaluate` or :meth:`Filter.filter` are: - -def compute_jackson_cheby_coeff(filter_bounds, delta_lambda, m): - r""" - To compute the m+1 coefficients of the polynomial approximation of an ideal band-pass between a and b, between a range of values defined by lambda_min and lambda_max. - - Parameters - ---------- - filter_bounds : list - [a, b] - delta_lambda : list - [lambda_min, lambda_max] - m : int - - Returns - ------- - ch : ndarray - jch : ndarray - - References - ---------- - :cite:`tremblay2016compressive` - - """ - # Parameters check - if delta_lambda[0] > filter_bounds[0] or delta_lambda[1] < filter_bounds[1]: - _logger.error("Bounds of the filter are out of the lambda values") - raise() - elif delta_lambda[0] > delta_lambda[1]: - _logger.error("lambda_min is greater than lambda_max") - raise() - - # Scaling and translating to standard cheby interval - a1 = (delta_lambda[1]-delta_lambda[0])/2 - a2 = (delta_lambda[1]+delta_lambda[0])/2 - - # Scaling bounds of the band pass according to lrange - filter_bounds[0] = (filter_bounds[0]-a2)/a1 - filter_bounds[1] = (filter_bounds[1]-a2)/a1 - - # First compute cheby coeffs - ch = np.arange(float(m+1)) - ch[0] = (2/(np.pi))*(np.arccos(filter_bounds[0])-np.arccos(filter_bounds[1])) - for i in ch[1:]: - ch[i] = (2/(np.pi * i)) * \ - (np.sin(i * np.arccos(filter_bounds[0])) - np.sin(i * np.arccos(filter_bounds[1]))) - - # Then compute jackson coeffs - jch = np.arange(float(m+1)) - alpha = (np.pi/(m+2)) - for i in jch: - jch[i] = (1/np.sin(alpha)) * \ - ((1 - i/(m+2)) * np.sin(alpha) * np.cos(i * alpha) + - (1/(m+2)) * np.cos(alpha) * np.sin(i * alpha)) - - # Combine jackson and cheby coeffs - jch = ch * jch - - return ch, jch - - -def lanczos_op(f, s, order=30): - r""" - Perform the lanczos approximation of the signal s. + * recursive, defined + * direct, which returns :math:`\sum_k c_k T_k(x) s = \sum_k c_k \cos(k \arccos x) s`. Parameters ---------- - f: Filter - s : ndarray - Signal to approximate. + G : graph + filters : Filter or array-like + Either a :class:`Filter` object or a set of Chebyshev coefficients + represented as an array of size K x F, where K is the polynomial + order and F the number of filters. + K x Fout x Fin + For convenience, Fin and Fout can be omitted. order : int - Degree of the lanczos approximation. (default = 30) - - Returns - ------- - L : ndarray - lanczos approximation of s + Polynomial order. + + Examples + -------- + + Plot the basis formed by the first K Chebyshev polynomials: + + >>> import matplotlib.pyplot as plt + >>> fig, axes = plt.subplots(1, 2) + >>> + >>> G = graphs.Ring(N=20) + >>> G.compute_fourier_basis() # To be exactly orthogonal in vertex domain. + >>> G.set_coordinates('line1D') + >>> + >>> K = 5 # Polynomials of order up to K. + >>> + >>> coefficients = np.identity(K) + >>> f = filters.Chebyshev(G, coefficients) + >>> s = f.localize(G.N // 2) + >>> _ = f.plot(sum=False, eigenvalues=False, ax=axes[0]) + >>> _ = G.plot_signal(s.T, ax=axes[1]) + >>> + >>> _ = axes[0].set_title('Chebysev polynomials in the spectral domain') + >>> _ = axes[1].set_title('Chebysev polynomials in the ring graph domain') + >>> _ = axes[0].legend(['order {}'.format(order) for order in range(K)]) + >>> _ = axes[1].legend(['order {}'.format(order) for order in range(K)]) + + They are orthogonal in the vertex domain: + + >>> s = s.T.reshape((G.N, -1)) + >>> print(s.T.dot(s)) + [[20. 0. 0. 0. 0.] + [ 0. 10. 0. 0. 0.] + [ 0. 0. 10. 0. 0.] + [ 0. 0. 0. 10. 0.] + [ 0. 0. 0. 0. 10.]] """ - G = f.G - Nf = len(f.g) - - # To have the right shape for the output array depending on the signal dim - try: - Nv = np.shape(s)[1] - is2d = True - c = np.zeros((G.N*Nf, Nv)) - except IndexError: - Nv = 1 - is2d = False - c = np.zeros((G.N*Nf)) - - tmpN = np.arange(G.N, dtype=int) - for j in range(Nv): - if is2d: - V, H, _ = lanczos(G.L.toarray(), order, s[:, j]) - else: - V, H, _ = lanczos(G.L.toarray(), order, s) - - Eh, Uh = np.linalg.eig(H) - - Eh[Eh < 0] = 0 - fe = f.evaluate(Eh) - V = np.dot(V, Uh) - - for i in range(Nf): - if is2d: - c[tmpN + i*G.N, j] = np.dot(V, fe[:][i] * np.dot(V.T, s[:, j])) - else: - c[tmpN + i*G.N] = np.dot(V, fe[:][i] * np.dot(V.T, s)) - - return c - - -def lanczos(A, order, x): - r""" - TODO short description - Parameters - ---------- - A: ndarray - - Returns - ------- - """ - try: - N, M = np.shape(x) - except ValueError: - N = np.shape(x)[0] - M = 1 - x = x[:, np.newaxis] - - # normalization - q = np.divide(x, np.kron(np.ones((N, 1)), np.linalg.norm(x, axis=0))) - - # initialization - hiv = np.arange(0, order*M, order) - V = np.zeros((N, M*order)) - V[:, hiv] = q + def __init__(self, G, coefficients): + + self.G = G + + coefficients = np.asarray(coefficients) + while coefficients.ndim < 3: + coefficients = np.expand_dims(coefficients, -1) + + self.n_features_out, self.n_features_in = coefficients.shape[1:] + self.shape = (self.n_features_in, self.n_features_out) # TODO: useful? + self.n_filters = self.n_features_in * self.n_features_out + self.Nf = self.n_filters # TODO: kept for backward compatibility only. + self._coefficients = coefficients + + # That is a factory method. + @classmethod + def from_filter(cls, filters, order=30, n=None): + r"""Compute the Chebyshev coefficients which approximate the filters. + + The :math:`K+1` coefficients, where :math:`K` is the polynomial order, + to approximate the function :math:`f` are computed by the discrete + orthogonality condition as + + .. math:: a_k \approx \frac{2-\delta_{0k}}{N} + \sum_{n=0}^{N-1} T_k(x_n) f(x_n), + + where :math:`\delta_{ij}` is the Kronecker delta function and the + :math:`x_n` are the N shifted Gauss–Chebyshev zeros of :math:`T_N(x)`, + given by + + .. math:: x_n = \frac{\lambda_\text{max}}{2} + \cos\left( \frac{\pi (2k+1)}{2N} \right) + + \frac{\lambda_\text{max}}{2}. + + For any N, these approximate coefficients provide an exact + approximation to the function at :math:`x_k` with a controlled error + between those points. The exact coefficients are obtained with + :math:`N=\infty`, thus representing the function exactly at all points + in :math:`[0, \lambda_\text{max}]`. The rate of convergence depends on + the function and its smoothness. + + Parameters + ---------- + filters : filters.Filter + A filterbank (:class:`Filter`) to be approximated by a set of + Chebyshev polynomials. + order : int + The order of the Chebyshev polynomials. + n : int + The number of Gauss–Chebyshev zeros used to approximate the + coefficients. Defaults to the polynomial order plus one. + + Examples + -------- + + Chebyshev coefficients which approximate a linear function: + + >>> G = graphs.Ring() + >>> G.estimate_lmax() + >>> g = filters.Filter(G, lambda x: 2*x) + >>> h = filters.Chebyshev.from_filter(g, order=4) + >>> print(', '.join([str(int(c)) for c in h._coefficients])) + 4, 4, 0, 0, 0 + + Coefficients of smooth filters decrease rapidly: + + >>> import matplotlib.pyplot as plt + >>> taus = [5, 10, 20, 50] + >>> g = filters.Heat(G, tau=taus) + >>> h = filters.Chebyshev.from_filter(g, order=10) + >>> fig, axes = plt.subplots(1, 2) + >>> _ = g.plot(sum=False, ax=axes[0]) + >>> _ = axes[1].plot(h._coefficients.squeeze()) + >>> _ = axes[0].legend(['tau = {}'.format(tau) for tau in taus]) + >>> _ = axes[1].legend(['tau = {}'.format(tau) for tau in taus]) + + """ + lmax = filters.G.lmax + + if n is None: + n = order + 1 + + points = np.pi * (np.arange(n) + 0.5) / n + + # The Gauss–Chebyshev zeros of Tk(x), scaled to [0, lmax]. + zeros = lmax/2 * np.cos(points) + lmax/2 + + # TODO: compute with scipy.fftpack.dct(). + c = np.empty((order + 1, filters.Nf)) + for i, kernel in enumerate(filters._kernels): + for k in range(order + 1): + T_k = np.cos(k * points) # Chebyshev polynomials of order k. + c[k, i] = 2 / n * kernel(zeros).dot(T_k) + c[0, :] /= 2 + + return cls(filters.G, c) + + @staticmethod + def scale_data(x, lmax): + r"""Given numbers in [0, lmax], scale them to [-1, 1].""" + + if x.min() < 0 or x.max() > lmax: + _logger.warning('You are trying to evaluate Chebyshev ' + 'polynomials outside of their orthonormal ' + 'domain [0, lmax={:.2f}].'.format(lmax)) + + return 2 * x / lmax - 1 + + @staticmethod + def scale_operator(L, lmax): + r"""Scale an operator's eigenvalues from [0, lmax] to [-1, 1].""" + if not sparse.issparse(L): + I = np.identity(L.shape[0], dtype=L.dtype) + else: + I = sparse.identity(L.shape[0], format=L.format, dtype=L.dtype) - H = np.zeros((order + 1, M*order)) - r = np.dot(A, q) - H[0, hiv] = np.sum(q*r, axis=0) - r -= np.kron(np.ones((N, 1)), H[0, hiv])*q - H[1, hiv] = np.linalg.norm(r, axis=0) + return 2 * L / lmax - I - orth = np.zeros((order)) - orth[0] = np.linalg.norm(np.dot(V.T, V) - M) + def _evaluate(self, x, method): - for k in range(1, order): - if np.sum(np.abs(H[k, hiv + k - 1])) <= np.spacing(1): - H = H[:k - 1, _sum_ind(np.arange(k), hiv)] - V = V[:, _sum_ind(np.arange(k), hiv)] - orth = orth[:k] + x = self.scale_data(x, self.G.lmax) - return V, H, orth + c = self._coefficients + K, Fout, Fin = c.shape # #order x #features_out x #features_in + c = c.reshape((K, Fout * Fin) + (1,) * x.ndim) # For broadcasting. - H[k - 1, hiv + k] = H[k, hiv + k - 1] - v = q - q = r/np.tile(H[k - 1, k + hiv], (N, 1)) - V[:, k + hiv] = q + # Recursive faster than direct faster than clenshaw. + method = 'recursive' if method is None else method - r = np.dot(A, q) - r -= np.tile(H[k - 1, k + hiv], (N, 1))*v - H[k, k + hiv] = np.sum(np.multiply(q, r), axis=0) - r -= np.tile(H[k, k + hiv], (N, 1))*q + try: + y = getattr(self, '_evaluate_' + method)(c, x) + except AttributeError: + raise ValueError('Unknown method {}.'.format(method)) - # The next line has to be checked - r -= np.dot(V, np.dot(V.T, r)) # full reorthogonalization - H[k + 1, k + hiv] = np.linalg.norm(r, axis=0) - orth[k] = np.linalg.norm(np.dot(V.T, V) - M) + return y.reshape((Fout, Fin) + x.shape).squeeze() - H = H[np.ix_(np.arange(order), np.arange(order))] + def _filter(self, s, method, _): - return V, H, orth + # TODO: signal normalization will move to Filter.filter() + # Dimension 3: number of nodes. + if s.shape[-1] != self.G.N: + raise ValueError('The last dimension should be {}, ' + 'the number of nodes. ' + 'Got instead a signal of shape ' + '{}.'.format(self.G.N, s.shape)) -def _sum_ind(ind1, ind2): - ind = np.tile(np.ravel(ind1), (np.size(ind2), 1)).T + np.ravel(ind2) - return np.ravel(ind) + # Dimension 2: number of input features. + if s.ndim == 1: + s = np.expand_dims(s, 0) + if s.shape[-2] != self.n_features_in: + if self.n_features_in == 1 and s.ndim == 2: + # Dimension can be omitted if there's 1 input feature. + s = np.expand_dims(s, -2) + else: + raise ValueError('The second to last dimension should be {}, ' + 'the number of input features. ' + 'Got instead a signal of shape ' + '{}.'.format(self.n_features_in, s.shape)) + + # Dimension 1: number of independent signals. + if s.ndim < 3: + s = np.expand_dims(s, 0) + + if s.ndim > 3: + raise ValueError('Signals should have at most 3 dimensions: ' + '#signals x #features x #nodes.') + + assert s.ndim == 3 + assert s.shape[2] == self.G.N # Number of nodes. + assert s.shape[1] == self.n_features_in # Number of input features. + # n_signals = s.shape[0] + + L = self.scale_operator(self.G.L, self.G.lmax) + + # Recursive and clenshaw are similarly fast. + method = 'recursive' if method is None else method + + try: + return getattr(self, '_filter_' + method)(L, s) + except AttributeError: + raise ValueError('Unknown method {}.'.format(method)) + + def _evaluate_direct(self, c, x): + r"""Evaluate Fout*Fin polynomials at each value in x.""" + K, F = c.shape[:2] + result = np.zeros((F,) + x.shape) + x = np.arccos(x) + for k in range(K): + result += c[k] * np.cos(k * x) + return result + + def _evaluate_recursive(self, c, x): + """Evaluate a Chebyshev series for y. Optionally, times s. + + .. math:: p(x) = \sum_{k=0}^{K} a_k * T_k(x) * s + + Parameters + ---------- + c: array-like + set of Chebyshev coefficients. (size K x F where K is the polynomial order, F is the number of filters) + x: array-like + vector to be evaluated. (size N x 1) + vector or matrix + signal: array-like + signal (vector) to be multiplied to the result. It allows to avoid the computation of powers of matrices when what we care about is L^k s not L^k. + vector or matrix (ndarray) + + Returns + ------- + corresponding Chebyshev Series. (size F x N) + + """ + K = c.shape[0] + x0 = np.ones_like(x) + result = c[0] * x0 + if K > 1: + x1 = x + result += c[1] * x1 + for k in range(2, K): + x2 = 2 * x * x1 - x0 + result += c[k] * x2 + x0, x1 = x1, x2 + return result + + def _filter_recursive(self, L, s): + r"""Filter a signal with the 3-way recursive relation. + Time: O(M N Fin Fout K) + Space: O(4 M Fin N) + """ + c = self._coefficients + K, Fout, Fin = c.shape # #order x #features_out x #features_in + M, Fin, N = s.shape # #signals x #features x #nodes + + def mult(c, x): + """Multiply the diffused signals by the Chebyshev coefficients.""" + x.shape = (M, Fin, N) + return np.einsum('oi,min->mon', c, x) + def dot(L, x): + """One diffusion step by multiplication with the Laplacian.""" + x.shape = (M * Fin, N) + return L.__rmatmul__(x) # x @ L + + x0 = s.view() + result = mult(c[0], x0) + if K > 1: + x1 = dot(L, x0) + result += mult(c[1], x1) + for k in range(2, K): + x2 = 2 * dot(L, x1) - x0 + result += mult(c[k], x2) + x0, x1 = x1, x2 + return result + + def _filter_clenshaw(self, L, s): + c = self._coefficients + K, Fout, Fin = c.shape # #order x #features_out x #features_in + M, Fin, N = s.shape # #signals x #features x #nodes + + def mult(c, s): + """Multiply the signals by the Chebyshev coefficients.""" + return np.einsum('oi,min->mon', c, s) + def dot(L, x): + """One diffusion step by multiplication with the Laplacian.""" + x.shape = (M * Fout, N) + y = L.__rmatmul__(x) # x @ L + x.shape = (M, Fout, N) + y.shape = (M, Fout, N) + return y + + b2 = 0 + b1 = mult(c[K-1], s) if K >= 2 else np.zeros((M, Fout, N)) + for k in range(K-2, 0, -1): + b = mult(c[k], s) + 2 * dot(L, b1) - b2 + b2, b1 = b1, b + return mult(c[0], s) + dot(L, b1) - b2 + + def _evaluate_clenshaw(self, c, x): + K = c.shape[0] + b2 = 0 + b1 = c[K-1] * np.ones_like(x) if K >= 2 else 0 + for k in range(K-2, 0, -1): + b = c[k] + 2 * x * b1 - b2 + b2, b1 = b1, b + return c[0] + x * b1 - b2 diff --git a/pygsp/filters/approximations_old.py b/pygsp/filters/approximations_old.py new file mode 100644 index 00000000..905ed3e8 --- /dev/null +++ b/pygsp/filters/approximations_old.py @@ -0,0 +1,343 @@ +# -*- coding: utf-8 -*- + +import numpy as np +from scipy import sparse + +from pygsp import utils + + +_logger = utils.build_logger(__name__) + + +@utils.filterbank_handler +def compute_cheby_coeff(f, m=30, N=None, *args, **kwargs): + r""" + Compute Chebyshev coefficients for a Filterbank. + + Parameters + ---------- + f : Filter + Filterbank with at least 1 filter + m : int + Maximum order of Chebyshev coeff to compute + (default = 30) + N : int + Grid order used to compute quadrature + (default = m + 1) + i : int + Index of the Filterbank element to compute + (default = 0) + + Returns + ------- + c : ndarray + Matrix of Chebyshev coefficients + + """ + G = f.G + i = kwargs.pop('i', 0) + + if not N: + N = m + 1 + + a_arange = [0, G.lmax] + + a1 = (a_arange[1] - a_arange[0]) / 2 + a2 = (a_arange[1] + a_arange[0]) / 2 + c = np.zeros(m + 1) + + tmpN = np.arange(N) + num = np.cos(np.pi * (tmpN + 0.5) / N) + for o in range(m + 1): + c[o] = 2. / N * np.dot(f._kernels[i](a1 * num + a2), + np.cos(np.pi * o * (tmpN + 0.5) / N)) + + return c + + +def cheby_op(G, c, signal, **kwargs): + r""" + Chebyshev polynomial of graph Laplacian applied to vector. + + Parameters + ---------- + G : Graph + c : ndarray or list of ndarrays + Chebyshev coefficients for a Filter or a Filterbank + signal : ndarray + Signal to filter + + Returns + ------- + r : ndarray + Result of the filtering + + """ + # Handle if we do not have a list of filters but only a simple filter in cheby_coeff. + if not isinstance(c, np.ndarray): + c = np.array(c) + + c = np.atleast_2d(c) + Nscales, M = c.shape + + if M < 2: + raise TypeError("The coefficients have an invalid shape") + + # thanks to that, we can also have 1d signal. + try: + Nv = np.shape(signal)[1] + r = np.zeros((G.N * Nscales, Nv)) + except IndexError: + r = np.zeros((G.N * Nscales)) + + a_arange = [0, G.lmax] + + a1 = float(a_arange[1] - a_arange[0]) / 2. + a2 = float(a_arange[1] + a_arange[0]) / 2. + + twf_old = signal + twf_cur = (G.L.dot(signal) - a2 * signal) / a1 + + tmpN = np.arange(G.N, dtype=int) + for i in range(Nscales): + r[tmpN + G.N*i] = 0.5 * c[i, 0] * twf_old + c[i, 1] * twf_cur + + factor = 2/a1 * (G.L - a2 * sparse.eye(G.N)) + for k in range(2, M): + twf_new = factor.dot(twf_cur) - twf_old + for i in range(Nscales): + r[tmpN + G.N*i] += c[i, k] * twf_new + + twf_old = twf_cur + twf_cur = twf_new + + return r + + +def cheby_rect(G, bounds, signal, **kwargs): + r""" + Fast filtering using Chebyshev polynomial for a perfect rectangle filter. + + Parameters + ---------- + G : Graph + bounds : array-like + The bounds of the pass-band filter + signal : array-like + Signal to filter + order : int (optional) + Order of the Chebyshev polynomial (default: 30) + + Returns + ------- + r : array-like + Result of the filtering + + """ + if not (isinstance(bounds, (list, np.ndarray)) and len(bounds) == 2): + raise ValueError('Bounds of wrong shape.') + + bounds = np.array(bounds) + + m = int(kwargs.pop('order', 30) + 1) + + try: + Nv = np.shape(signal)[1] + r = np.zeros((G.N, Nv)) + except IndexError: + r = np.zeros((G.N)) + + b1, b2 = np.arccos(2. * bounds / G.lmax - 1.) + factor = 4./G.lmax * G.L - 2.*sparse.eye(G.N) + + T_old = signal + T_cur = factor.dot(signal) / 2. + r = (b1 - b2)/np.pi * signal + 2./np.pi * (np.sin(b1) - np.sin(b2)) * T_cur + + for k in range(2, m): + T_new = factor.dot(T_cur) - T_old + r += 2./(k*np.pi) * (np.sin(k*b1) - np.sin(k*b2)) * T_new + T_old = T_cur + T_cur = T_new + + return r + + +def compute_jackson_cheby_coeff(filter_bounds, delta_lambda, m): + r""" + To compute the m+1 coefficients of the polynomial approximation of an ideal band-pass between a and b, between a range of values defined by lambda_min and lambda_max. + + Parameters + ---------- + filter_bounds : list + [a, b] + delta_lambda : list + [lambda_min, lambda_max] + m : int + + Returns + ------- + ch : ndarray + jch : ndarray + + References + ---------- + :cite:`tremblay2016compressive` + + """ + # Parameters check + if delta_lambda[0] > filter_bounds[0] or delta_lambda[1] < filter_bounds[1]: + _logger.error("Bounds of the filter are out of the lambda values") + raise() + elif delta_lambda[0] > delta_lambda[1]: + _logger.error("lambda_min is greater than lambda_max") + raise() + + # Scaling and translating to standard cheby interval + a1 = (delta_lambda[1]-delta_lambda[0])/2 + a2 = (delta_lambda[1]+delta_lambda[0])/2 + + # Scaling bounds of the band pass according to lrange + filter_bounds[0] = (filter_bounds[0]-a2)/a1 + filter_bounds[1] = (filter_bounds[1]-a2)/a1 + + # First compute cheby coeffs + ch = np.arange(float(m+1)) + ch[0] = (2/(np.pi))*(np.arccos(filter_bounds[0])-np.arccos(filter_bounds[1])) + for i in ch[1:]: + ch[i] = (2/(np.pi * i)) * \ + (np.sin(i * np.arccos(filter_bounds[0])) - np.sin(i * np.arccos(filter_bounds[1]))) + + # Then compute jackson coeffs + jch = np.arange(float(m+1)) + alpha = (np.pi/(m+2)) + for i in jch: + jch[i] = (1/np.sin(alpha)) * \ + ((1 - i/(m+2)) * np.sin(alpha) * np.cos(i * alpha) + + (1/(m+2)) * np.cos(alpha) * np.sin(i * alpha)) + + # Combine jackson and cheby coeffs + jch = ch * jch + + return ch, jch + + +def lanczos_op(f, s, order=30): + r""" + Perform the lanczos approximation of the signal s. + + Parameters + ---------- + f: Filter + s : ndarray + Signal to approximate. + order : int + Degree of the lanczos approximation. (default = 30) + + Returns + ------- + L : ndarray + lanczos approximation of s + + """ + G = f.G + Nf = len(f.g) + + # To have the right shape for the output array depending on the signal dim + try: + Nv = np.shape(s)[1] + is2d = True + c = np.zeros((G.N*Nf, Nv)) + except IndexError: + Nv = 1 + is2d = False + c = np.zeros((G.N*Nf)) + + tmpN = np.arange(G.N, dtype=int) + for j in range(Nv): + if is2d: + V, H, _ = lanczos(G.L.toarray(), order, s[:, j]) + else: + V, H, _ = lanczos(G.L.toarray(), order, s) + + Eh, Uh = np.linalg.eig(H) + + Eh[Eh < 0] = 0 + fe = f.evaluate(Eh) + V = np.dot(V, Uh) + + for i in range(Nf): + if is2d: + c[tmpN + i*G.N, j] = np.dot(V, fe[:][i] * np.dot(V.T, s[:, j])) + else: + c[tmpN + i*G.N] = np.dot(V, fe[:][i] * np.dot(V.T, s)) + + return c + + +def lanczos(A, order, x): + r""" + TODO short description + + Parameters + ---------- + A: ndarray + + Returns + ------- + """ + try: + N, M = np.shape(x) + except ValueError: + N = np.shape(x)[0] + M = 1 + x = x[:, np.newaxis] + + # normalization + q = np.divide(x, np.kron(np.ones((N, 1)), np.linalg.norm(x, axis=0))) + + # initialization + hiv = np.arange(0, order*M, order) + V = np.zeros((N, M*order)) + V[:, hiv] = q + + H = np.zeros((order + 1, M*order)) + r = np.dot(A, q) + H[0, hiv] = np.sum(q*r, axis=0) + r -= np.kron(np.ones((N, 1)), H[0, hiv])*q + H[1, hiv] = np.linalg.norm(r, axis=0) + + orth = np.zeros((order)) + orth[0] = np.linalg.norm(np.dot(V.T, V) - M) + + for k in range(1, order): + if np.sum(np.abs(H[k, hiv + k - 1])) <= np.spacing(1): + H = H[:k - 1, _sum_ind(np.arange(k), hiv)] + V = V[:, _sum_ind(np.arange(k), hiv)] + orth = orth[:k] + + return V, H, orth + + H[k - 1, hiv + k] = H[k, hiv + k - 1] + v = q + q = r/np.tile(H[k - 1, k + hiv], (N, 1)) + V[:, k + hiv] = q + + r = np.dot(A, q) + r -= np.tile(H[k - 1, k + hiv], (N, 1))*v + H[k, k + hiv] = np.sum(np.multiply(q, r), axis=0) + r -= np.tile(H[k, k + hiv], (N, 1))*q + + # The next line has to be checked + r -= np.dot(V, np.dot(V.T, r)) # full reorthogonalization + H[k + 1, k + hiv] = np.linalg.norm(r, axis=0) + orth[k] = np.linalg.norm(np.dot(V.T, V) - M) + + H = H[np.ix_(np.arange(order), np.arange(order))] + + return V, H, orth + + +def _sum_ind(ind1, ind2): + ind = np.tile(np.ravel(ind1), (np.size(ind2), 1)).T + np.ravel(ind2) + return np.ravel(ind) diff --git a/pygsp/filters/filter.py b/pygsp/filters/filter.py index 544fd331..39aa2df6 100644 --- a/pygsp/filters/filter.py +++ b/pygsp/filters/filter.py @@ -1,10 +1,11 @@ # -*- coding: utf-8 -*- import numpy as np +from scipy import sparse from pygsp import utils # prevent circular import in Python < 3.5 -from . import approximations +from . import approximations_old as approximations _logger = utils.build_logger(__name__) @@ -61,8 +62,9 @@ def __init__(self, G, kernels): kernels = [kernels] self._kernels = kernels - # Only used by subclasses to instantiate a single filterbank. + # This constructor is only used by subclasses to instantiate a single filterbank. self.n_features_in, self.n_features_out = (1, len(kernels)) + self.shape = (self.n_features_in, self.n_features_out) # TODO: useful? self.n_filters = self.n_features_in * self.n_features_out self.Nf = self.n_filters # TODO: kept for backward compatibility only. @@ -78,7 +80,7 @@ def __repr__(self): s += '{}={}, '.format(key, value) return '{}({})'.format(self.__class__.__name__, s[:-2]) - def evaluate(self, x): + def evaluate(self, x, method=None): r"""Evaluate the kernels at given frequencies. Parameters @@ -105,13 +107,65 @@ def evaluate(self, x): [] """ + x = np.asarray(x) # For iterables. Sparse makes no sense here. + return self._evaluate(x, method) + + def _evaluate(self, x, _): + r"""Default implementation for filters defined as kernel functions.""" # Avoid to copy data as with np.array([g(x) for g in self._kernels]). y = np.empty([self.Nf] + list(x.shape)) for i, kernel in enumerate(self._kernels): y[i] = kernel(x) return y - def filter(self, s, method='chebyshev', order=30): + def approximate(self, method, **kwargs): + r"""Return a filter which approximates this filter. + + While approximations might loose accuracy, they allow for much faster + computations. + + Parameters + ---------- + method : str + Approximation method. Only 'Chebyshev' is supported for now. + kwargs : dict + Parameters for the approximation method. + + Examples + -------- + + Approximate a filter with Chebyshev polynomials of various orders: + + >>> import matplotlib.pyplot as plt + >>> fig, ax = plt.subplots(1, 1) + >>> + >>> G = graphs.Ring() + >>> G.compute_fourier_basis() + >>> f1 = filters.Heat(G) + >>> _ = f1.plot(eigenvalues=True, linewidth=3, label='continuous', ax=ax) + >>> + >>> for order in range(1, 5): + ... f2 = f1.approximate('Chebyshev', order=order) + ... l = 'Chebyshev order {}'.format(order) + ... _ = f2.plot(eigenvalues=False, label=l, linestyle='dashed', ax=ax) + >>> + >>> _ = ax.set_title('Approximation for various polynomial orders') + >>> _ = ax.legend() + + Approximate a filterbank with Chebyshev polynomials: + + >>> G = graphs.Ring() + >>> G.compute_fourier_basis() + >>> f1 = filters.Itersine(G) + >>> f2 = f1.approximate('Chebyshev', order=20) + >>> _ = f1.plot(title='Continuous filterbank') + >>> _ = f2.plot(title='Approximated filterbank') + + """ + from . import approximations + return getattr(approximations, method).from_filter(self, **kwargs) + + def filter(self, s, method=None, order=30): r"""Filter signals (analysis or synthesis). A signal is defined as a rank-3 tensor of shape ``(N_NODES, N_SIGNALS, @@ -164,7 +218,7 @@ def filter(self, s, method='chebyshev', order=30): s : ndarray Graph signals, a tensor of shape ``(N_NODES, N_SIGNALS, N_FEATURES)``, where ``N_NODES`` and ``N_SIGNALS`` are the number - of nodes and signals of the signal tensor that pas passed in, and + of nodes and signals of the signal tensor that was passed in, and ``N_FEATURES`` is either 1 (synthesis) or the number of filters in the filter bank (analysis). @@ -232,6 +286,19 @@ def filter(self, s, method='chebyshev', order=30): True """ + if not sparse.issparse(s): + s = np.asarray(s) # For iterables. + + s = self._filter(s, method, order) + + # Return a 1D signal if e.g. a 1D signal was filtered by one filter. + return s.squeeze() + + def _filter(self, s, method='chebyshev', order=30): + r"""Default implementation for filters defined as kernel functions.""" + + method = 'chebyshev' if method is None else method + if s.shape[0] != self.G.N: raise ValueError('First dimension should be the number of nodes ' 'G.N = {}, got {}.'.format(self.G.N, s.shape)) @@ -296,8 +363,7 @@ def filter(self, s, method='chebyshev', order=30): else: raise ValueError('Unknown method {}.'.format(method)) - # Return a 1D signal if e.g. a 1D signal was filtered by one filter. - return s.squeeze() + return s def analyze(self, s, method='chebyshev', order=30): r"""Convenience alias to :meth:`filter`.""" @@ -354,8 +420,12 @@ def localize(self, i, **kwargs): >>> _ = G.plot(s, highlight=DELTA) """ - s = np.zeros(self.G.N) - s[i] = 1 + # Localize each filter: g[in, out].localize(i) for all (in, out). + s = np.zeros((self.n_features_in, self.n_features_in, self.G.N)) + fin = range(self.n_features_in) + s[fin, fin, i] = 1 + # TODO: remove once all signals have the new shape. + s = s.squeeze() return np.sqrt(self.G.N) * self.filter(s, **kwargs) def estimate_frame_bounds(self, min=0, max=None, N=1000, diff --git a/pygsp/graphs/graph.py b/pygsp/graphs/graph.py index 0f10a69f..39e4362d 100644 --- a/pygsp/graphs/graph.py +++ b/pygsp/graphs/graph.py @@ -716,7 +716,7 @@ def lmax(self): self.estimate_lmax() return self._lmax - def estimate_lmax(self, recompute=False): + def estimate_lmax(self, tol=5e-3, recompute=False): r"""Estimate the Laplacian's largest eigenvalue (cached). The result is cached and accessible by the :attr:`lmax` property. @@ -727,6 +727,10 @@ def estimate_lmax(self, recompute=False): Parameters ---------- + tol : float + Relative accuracy when computing eigenvalues (stopping criterion). + The larger the tolerance the faster the computation but the less + accurate the estimation. A value of 0 implies machine precision. recompute : boolean Force to recompute the largest eigenvalue. Default is false. @@ -747,21 +751,23 @@ def estimate_lmax(self, recompute=False): >>> G.compute_fourier_basis() >>> print('{:.2f}'.format(G.lmax)) 13.78 - >>> G = graphs.Logo() >>> G.estimate_lmax(recompute=True) >>> print('{:.2f}'.format(G.lmax)) 13.92 + >>> G.estimate_lmax(tol=0, recompute=True) + >>> print('{:.2f}'.format(G.lmax)) + 13.78 """ if hasattr(self, '_lmax') and not recompute: return try: - lmax = sparse.linalg.eigsh(self.L, k=1, tol=5e-3, + lmax = sparse.linalg.eigsh(self.L, k=1, tol=tol, ncv=min(self.N, 10), return_eigenvectors=False) lmax = lmax[0] - lmax *= 1.01 # Increase by 1 percent to be robust to errors. + lmax *= 1 + 2*tol # Be robust to errors. except sparse.linalg.ArpackNoConvergence: self.logger.warning('Lanczos method did not converge. ' diff --git a/pygsp/tests/test_all.py b/pygsp/tests/test_all.py index 94ad4c94..bad490b3 100755 --- a/pygsp/tests/test_all.py +++ b/pygsp/tests/test_all.py @@ -14,10 +14,13 @@ from pygsp.tests import test_docstrings from pygsp.tests import test_plotting +loader = unittest.TestLoader() suites = [] + suites.append(test_graphs.suite) suites.append(test_filters.suite) +suites.append(test_filters.suite_approximations) suites.append(test_utils.suite) suites.append(test_docstrings.suite) suites.append(test_plotting.suite) # TODO: can SIGSEGV if not last diff --git a/pygsp/tests/test_filters.py b/pygsp/tests/test_filters.py index 224f3905..da7b3bdd 100644 --- a/pygsp/tests/test_filters.py +++ b/pygsp/tests/test_filters.py @@ -273,4 +273,160 @@ def test_approximations(self): self.assertRaises(ValueError, f.filter, self._signal, method='lanczos') +class TestApproximations(unittest.TestCase): + + @classmethod + def setUpClass(cls): + cls._G = graphs.Logo() + cls._G.compute_fourier_basis() + cls._rs = np.random.RandomState(42) + cls._signal = cls._rs.uniform(size=cls._G.N) + + def test_scaling(self, N=40): + x = np.linspace(0, self._G.lmax, N) + x = filters.Chebyshev.scale_data(x, self._G.lmax) + self.assertEqual(x.min(), -1) + self.assertEqual(x.max(), 1) + L = filters.Chebyshev.scale_operator(self._G.L, self._G.lmax) + + def test_chebyshev_basis(self, K=5, c=2, N=100): + r""" + Test that the evaluation of the Chebyshev series yields the expected + basis. We only test the first two elements here. The higher-order + polynomials are compared with the trigonometric definition. + """ + f = filters.Chebyshev(self._G, c * np.identity(K)) + x = np.linspace(0, self._G.lmax, N) + y = f.evaluate(x) + np.testing.assert_equal(y[0], c) + np.testing.assert_allclose(y[1], np.linspace(-c, c, N)) + + def test_evaluation_methods(self, K=30, F=5, N=100): + r"""Test that all evaluation methods return the same results.""" + coefficients = self._rs.uniform(size=(K, F, F)) + f = filters.Chebyshev(self._G, coefficients) + x = np.linspace(0, self._G.lmax, N) + y1 = f.evaluate(x, method='recursive') + y2 = f.evaluate(x, method='direct') + y3 = f.evaluate(x, method='clenshaw') + np.testing.assert_allclose(y1, y2) + np.testing.assert_allclose(y1, y3) + # Evaluate on n-dimensional arrays. + x = self._rs.uniform(0, self._G.lmax, size=(3, 1, 19)) + y1 = f.evaluate(x, method='recursive') + y2 = f.evaluate(x, method='direct') + y3 = f.evaluate(x, method='clenshaw') + np.testing.assert_allclose(y1, y2) + np.testing.assert_allclose(y1, y3) + # Unknown method. + self.assertRaises(ValueError, f.evaluate, x, method='unk') + + def test_filter_identity(self, M=10, c=2.3): + r"""Test that filtering with c0 only scales the signal.""" + x = self._rs.uniform(size=(M, self._G.N)) + f = filters.Chebyshev(self._G, c) + y = f.filter(x, method='recursive') + np.testing.assert_equal(y, c * x) + # Test with dense Laplacian. + L = self._G.L + self._G.L = L.toarray() + y = f.filter(x, method='recursive') + self._G.L = L + np.testing.assert_equal(y, c * x) + + def test_filter_methods(self, K=30, Fin=5, Fout=6, M=100): + r"""Test that all filter methods return the same results.""" + coefficients = self._rs.uniform(size=(K, Fout, Fin)) + x = self._rs.uniform(size=(M, Fin, self._G.N)) + f = filters.Chebyshev(self._G, coefficients) + y1 = f.filter(x, method='recursive') + y2 = f.filter(x, method='clenshaw') + self.assertTupleEqual(y1.shape, (M, Fout, self._G.N)) + np.testing.assert_allclose(y1, y2, rtol=1e-5) + # Unknown method. + self.assertRaises(ValueError, f.filter, x, method='unk') + + def test_coefficients(self, K=10, slope=3.14): + r"""Test that the computed Chebyshev coefficients are correct.""" + # Identity filter. + f = filters.Heat(self._G, tau=0) + f = filters.Chebyshev.from_filter(f, order=K) + c = f._coefficients.squeeze() + np.testing.assert_allclose(c, [1] + K*[0], atol=1e-12) + # Linear filter. + f = filters.Filter(self._G, lambda x: slope*x) + f = filters.Chebyshev.from_filter(f, order=K) + c1 = f._coefficients.squeeze() + c2 = slope / 2 * self._G.lmax + np.testing.assert_allclose(c1, 2*[c2] + (K-1)*[0], atol=1e-12) + + def test_approximations(self, N=100, K=20): + r"""Test that the approximations are not far from the exact filters.""" + # Evaluation. + x = self._rs.uniform(0, self._G.lmax, N) + f1 = filters.Heat(self._G) + y1 = f1.evaluate(x) + f2 = f1.approximate('Chebyshev', order=K) + y2 = f2.evaluate(x) + np.testing.assert_allclose(y2, y1.squeeze()) + # Filtering. + x = self._rs.uniform(size=(1, 1, self._G.N)) + y1 = f1.filter(x.T).T + y2 = f2.filter(x) + np.testing.assert_allclose(y2.squeeze(), y1) + + def test_shape_normalization(self): + """Test that signal's shapes are properly normalized.""" + # TODO: should also test filters which are not approximations. + + def test_normalization(M, Fin, Fout, K=7): + + def test_shape(y, M, Fout, N=self._G.N): + """Test that filtered signals are squeezed.""" + if Fout == 1 and M == 1: + self.assertEqual(y.shape, (N,)) + elif Fout == 1: + self.assertEqual(y.shape, (M, N)) + elif M == 1: + self.assertEqual(y.shape, (Fout, N)) + else: + self.assertEqual(y.shape, (M, Fout, N)) + + coefficients = self._rs.uniform(size=(K, Fout, Fin)) + f = filters.Chebyshev(self._G, coefficients) + assert f.shape == (Fin, Fout) + assert (f.n_features_in, f.n_features_out) == (Fin, Fout) + + x = self._rs.uniform(size=(M, Fin, self._G.N)) + y = f.filter(x) + test_shape(y, M, Fout) + + if Fin == 1 or M == 1: + # It only makes sense to squeeze if one dimension is unitary. + x = x.squeeze() + y = f.filter(x) + test_shape(y, M, Fout) + + # Test all possible correct combinations of input and output signals. + for M in [1, 9]: + for Fin in [1, 3]: + for Fout in [1, 5]: + test_normalization(M, Fin, Fout) + + # Test failure cases. + M, Fin, Fout, K = 9, 3, 5, 7 + coefficients = self._rs.uniform(size=(K, Fout, Fin)) + f = filters.Chebyshev(self._G, coefficients) + x = self._rs.uniform(size=(M, Fin, 2)) + self.assertRaises(ValueError, f.filter, x) + x = self._rs.uniform(size=(M, 2, self._G.N)) + self.assertRaises(ValueError, f.filter, x) + x = self._rs.uniform(size=(2, self._G.N)) + self.assertRaises(ValueError, f.filter, x) + x = self._rs.uniform(size=(self._G.N)) + self.assertRaises(ValueError, f.filter, x) + x = self._rs.uniform(size=(2, M, Fin, self._G.N)) + self.assertRaises(ValueError, f.filter, x) + suite = unittest.TestLoader().loadTestsFromTestCase(TestCase) +suite_approximations = unittest.TestLoader().loadTestsFromTestCase(TestApproximations) diff --git a/pygsp/tests/test_graphs.py b/pygsp/tests/test_graphs.py index d2cdb87c..cd183971 100644 --- a/pygsp/tests/test_graphs.py +++ b/pygsp/tests/test_graphs.py @@ -360,5 +360,4 @@ def test_imgpatches(self): def test_grid2dimgpatches(self): graphs.Grid2dImgPatches(img=self._img, patch_shape=(3, 3)) - -suite = unittest.TestLoader().loadTestsFromTestCase(TestCase) +suite = unittest.TestLoader().loadTestsFromTestCase(TestCase) \ No newline at end of file diff --git a/pygsp/tests/test_plotting.py b/pygsp/tests/test_plotting.py index 936415d9..a823b0d0 100644 --- a/pygsp/tests/test_plotting.py +++ b/pygsp/tests/test_plotting.py @@ -98,6 +98,7 @@ def test(G): G = graphs.Torus(Nv=5) test(G) + def test_indices(self): def test(G): diff --git a/pygsp/tests/test_utils.py b/pygsp/tests/test_utils.py index cbd824b1..670510c2 100644 --- a/pygsp/tests/test_utils.py +++ b/pygsp/tests/test_utils.py @@ -130,5 +130,4 @@ def test_distanz(x, y): # test_distanz(x, y) - -suite = unittest.TestLoader().loadTestsFromTestCase(TestCase) +suite = unittest.TestLoader().loadTestsFromTestCase(TestCase) \ No newline at end of file