From 69aff82a512726e80337f536736d530d2cb26848 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Sun, 11 Mar 2018 13:35:42 +0100 Subject: [PATCH 01/25] motivation to represent approximated filters as proper filters --- doc/history.rst | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/doc/history.rst b/doc/history.rst index e71e00bf..b075a568 100644 --- a/doc/history.rst +++ b/doc/history.rst @@ -2,6 +2,26 @@ History ======= +0.6.0 (2018-03-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. + 0.5.1 (2017-12-15) ------------------ From 22e335ed90d0fde8d8f55ddcdc76b34a48848f4a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Sun, 11 Mar 2018 14:28:49 +0100 Subject: [PATCH 02/25] move old implementation of approximate filtering --- pygsp/filters/__init__.py | 2 +- pygsp/filters/{approximations.py => approximations_old.py} | 0 pygsp/filters/filter.py | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) rename pygsp/filters/{approximations.py => approximations_old.py} (100%) diff --git a/pygsp/filters/__init__.py b/pygsp/filters/__init__.py index aaa330e6..de7209bb 100644 --- a/pygsp/filters/__init__.py +++ b/pygsp/filters/__init__.py @@ -113,4 +113,4 @@ __all__ = _FILTERS + _APPROXIMATIONS _utils.import_classes(_FILTERS, 'filters', 'filters') -_utils.import_functions(_APPROXIMATIONS, 'filters.approximations', 'filters') +_utils.import_functions(_APPROXIMATIONS, 'filters.approximations_old', 'filters') diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations_old.py similarity index 100% rename from pygsp/filters/approximations.py rename to pygsp/filters/approximations_old.py diff --git a/pygsp/filters/filter.py b/pygsp/filters/filter.py index b4d9a6e4..b42942f5 100644 --- a/pygsp/filters/filter.py +++ b/pygsp/filters/filter.py @@ -4,7 +4,7 @@ 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__) From 0ca71ba386ccb363cdff8b2fc293e590778b04c4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Sun, 11 Mar 2018 14:54:03 +0100 Subject: [PATCH 03/25] Chebyshev approximation: sketch implementation --- pygsp/filters/__init__.py | 8 +++- pygsp/filters/approximations.py | 69 +++++++++++++++++++++++++++++++++ pygsp/filters/filter.py | 9 +++++ 3 files changed, 85 insertions(+), 1 deletion(-) create mode 100644 pygsp/filters/approximations.py diff --git a/pygsp/filters/__init__.py b/pygsp/filters/__init__.py index de7209bb..39726683 100644 --- a/pygsp/filters/__init__.py +++ b/pygsp/filters/__init__.py @@ -67,6 +67,10 @@ **Chebyshev polynomials** +.. autosummary:: + + Chebyshev + .. autosummary:: compute_cheby_coeff @@ -110,7 +114,9 @@ 'lanczos_op' ] -__all__ = _FILTERS + _APPROXIMATIONS +__all__ = _FILTERS + _APPROXIMATIONS + ['Chebyshev'] _utils.import_classes(_FILTERS, 'filters', 'filters') _utils.import_functions(_APPROXIMATIONS, 'filters.approximations_old', 'filters') + +from .approximations import * diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py new file mode 100644 index 00000000..e4a869d9 --- /dev/null +++ b/pygsp/filters/approximations.py @@ -0,0 +1,69 @@ +# -*- 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__) + + +class Chebyshev(Filter): + r"""Approximate continuous filters with Chebyshev polynomials. + + 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. + + Math to show how the coefficients are computed + + Parameters + ---------- + 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. + order : int + Polynomial order. + + """ + + def __init__(self, G, filters, order=30): + + self.G = G + self.order = order + + try: + self._compute_coefficients(filters) + self.Nf = filters.Nf + except: + self._coefficients = np.asarray(filters) + self.Nf = self._coefficients[1] + + def _evaluate(self, x, method='recursive'): + + if x.min() < 0 or x.max() > self.G.lmax: + _logger.warning('You are trying to evaluate Chebyshev ' + 'polynomials outside of their orthonormal ' + 'domain [0, {:.2f}]'.format(self.G.lmax)) + + x = np.asarray(x) + x = 2 * x / self.G.lmax - 1 # [0, lmax] => [-1, 1] + + return self._evaluate_polynomials(x) + + def _filter(self, s, method='recursive', _=None): + # method = 'clenshaw' in constructor or filter? + L = rescale_L(G.L) + return self._apply_polynomials(L, s) + + def _compute_coefficients(self, filters): + pass + + def _evaluate_polynomials(self, y, s=1): + pass diff --git a/pygsp/filters/filter.py b/pygsp/filters/filter.py index b42942f5..0913a590 100644 --- a/pygsp/filters/filter.py +++ b/pygsp/filters/filter.py @@ -87,6 +87,10 @@ def evaluate(self, x): [] """ + return self._evaluate(x) + + 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, len(x))) for i, kernel in enumerate(self._kernels): @@ -214,6 +218,11 @@ def filter(self, s, method='chebyshev', order=30): True """ + return self._filter(s, method, order) + + def _filter(self, s, method='chebyshev', order=30): + r"""Default implementation for filters defined as kernel functions.""" + 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)) From 8b3e6a7062844526f5619a1ae80a5a5da525fcd0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Sun, 11 Mar 2018 16:15:25 +0100 Subject: [PATCH 04/25] tests: only build the suites when running all tests The goal is to have multiple TestCase in one module without the need to concatenate them in one suite per module. --- pygsp/tests/test_all.py | 10 ++++++---- pygsp/tests/test_filters.py | 3 --- pygsp/tests/test_graphs.py | 3 --- pygsp/tests/test_plotting.py | 3 --- pygsp/tests/test_utils.py | 3 --- 5 files changed, 6 insertions(+), 16 deletions(-) diff --git a/pygsp/tests/test_all.py b/pygsp/tests/test_all.py index e6b675b5..869a8ee3 100755 --- a/pygsp/tests/test_all.py +++ b/pygsp/tests/test_all.py @@ -36,14 +36,16 @@ def setup(doctest): 'np': numpy, } +loader = unittest.TestLoader() suites = [] -suites.append(test_graphs.suite) -suites.append(test_filters.suite) -suites.append(test_utils.suite) +suites.append(loader.loadTestsFromModule(test_graphs)) +suites.append(loader.loadTestsFromModule(test_filters)) +suites.append(loader.loadTestsFromModule(test_utils)) suites.append(test_docstrings('pygsp', '.py', setup)) suites.append(test_docstrings('.', '.rst')) # No setup to not forget imports. -suites.append(test_plotting.suite) # TODO: can SIGSEGV if not last +# TODO: can SIGSEGV if not last. +suites.append(loader.loadTestsFromModule(test_plotting)) suite = unittest.TestSuite(suites) diff --git a/pygsp/tests/test_filters.py b/pygsp/tests/test_filters.py index 9807bd5a..18bd67c8 100644 --- a/pygsp/tests/test_filters.py +++ b/pygsp/tests/test_filters.py @@ -228,6 +228,3 @@ def test_approximations(self): np.testing.assert_allclose(c_exact, c_cheby) self.assertRaises(ValueError, f.filter, self._signal, method='lanczos') - - -suite = unittest.TestLoader().loadTestsFromTestCase(TestCase) diff --git a/pygsp/tests/test_graphs.py b/pygsp/tests/test_graphs.py index dc51fec5..66336cee 100644 --- a/pygsp/tests/test_graphs.py +++ b/pygsp/tests/test_graphs.py @@ -301,6 +301,3 @@ def test_imgpatches(self): def test_grid2dimgpatches(self): graphs.Grid2dImgPatches(img=self._img, patch_shape=(3, 3)) - - -suite = unittest.TestLoader().loadTestsFromTestCase(TestCase) diff --git a/pygsp/tests/test_plotting.py b/pygsp/tests/test_plotting.py index 6706b75d..f38ec1e3 100644 --- a/pygsp/tests/test_plotting.py +++ b/pygsp/tests/test_plotting.py @@ -105,6 +105,3 @@ def test(G): test(G) G = graphs.Torus(Nv=5) test(G) - - -suite = unittest.TestLoader().loadTestsFromTestCase(TestCase) diff --git a/pygsp/tests/test_utils.py b/pygsp/tests/test_utils.py index cbd824b1..dc585241 100644 --- a/pygsp/tests/test_utils.py +++ b/pygsp/tests/test_utils.py @@ -129,6 +129,3 @@ def test_distanz(x, y): # test_tree_depths(A, root) # test_distanz(x, y) - - -suite = unittest.TestLoader().loadTestsFromTestCase(TestCase) From 3e70aeeed4fbc284ca5ca1082ed43a7c893985d3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Sun, 11 Mar 2018 17:23:07 +0100 Subject: [PATCH 05/25] Chebyshev recursive evaluation --- pygsp/filters/approximations.py | 52 +++++++++++++++++++++++++++++---- pygsp/filters/filter.py | 2 ++ pygsp/tests/test_filters.py | 25 ++++++++++++++++ 3 files changed, 74 insertions(+), 5 deletions(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index e4a869d9..5051f48c 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -43,27 +43,69 @@ def __init__(self, G, filters, order=30): self.Nf = filters.Nf except: self._coefficients = np.asarray(filters) - self.Nf = self._coefficients[1] + self.Nf = self._coefficients.shape[1] def _evaluate(self, x, method='recursive'): if x.min() < 0 or x.max() > self.G.lmax: _logger.warning('You are trying to evaluate Chebyshev ' 'polynomials outside of their orthonormal ' - 'domain [0, {:.2f}]'.format(self.G.lmax)) + 'domain [0, {:.2f}].'.format(self.G.lmax)) - x = np.asarray(x) x = 2 * x / self.G.lmax - 1 # [0, lmax] => [-1, 1] return self._evaluate_polynomials(x) def _filter(self, s, method='recursive', _=None): # method = 'clenshaw' in constructor or filter? - L = rescale_L(G.L) + + M, M = L.shape + I = sparse.identity(M, format='csr', dtype=L.dtype) + L = 2 * L - self.G.lmax / 2 - I + return self._apply_polynomials(L, s) def _compute_coefficients(self, filters): + r"""Compute the coefficients of the Chebyshev series approximating the filters. + + Some implementations define c_0 / 2. + """ pass def _evaluate_polynomials(self, y, s=1): - pass + """Evaluate a Chebyshev series for y. Optionally, times s. + + .. math: p(y) = \sum_{k=0}^{K} a_k * T_k(y) * 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) + y: 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 = self._coefficients.shape[1] + c = self._coefficients + # Reshaping the coefficients to use broadcasting. + c = c.reshape(c.shape + (1,) * y.ndim) + + y0 = np.ones_like(y) + result = c[0] * y0.dot(s) + if K > 1: + y1 = y + result += c[1] * y1 + for k in range(2, K): + y2 = 2 * y.dot(s) * y1 - y0 + result += c[k] * y2 + y0, y1 = y1, y2 + return result diff --git a/pygsp/filters/filter.py b/pygsp/filters/filter.py index 0913a590..c619e22c 100644 --- a/pygsp/filters/filter.py +++ b/pygsp/filters/filter.py @@ -87,6 +87,7 @@ def evaluate(self, x): [] """ + x = np.asarray(x) return self._evaluate(x) def _evaluate(self, x): @@ -218,6 +219,7 @@ def filter(self, s, method='chebyshev', order=30): True """ + s = np.asarray(s) return self._filter(s, method, order) def _filter(self, s, method='chebyshev', order=30): diff --git a/pygsp/tests/test_filters.py b/pygsp/tests/test_filters.py index 18bd67c8..3af1244c 100644 --- a/pygsp/tests/test_filters.py +++ b/pygsp/tests/test_filters.py @@ -228,3 +228,28 @@ def test_approximations(self): np.testing.assert_allclose(c_exact, c_cheby) 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_chebyshev_basis(self): + r""" + Test that the evaluation of the Chebyshev series yields the expected + basis. + """ + K = 5 + c = 2 + f = filters.Chebyshev(self._G, c * np.identity(K)) + N = 100 + 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, 100)) + # TODO: do it for higher orders with the trigonometric definition From 2c057344cfcee05671af354bb8f308fd4122f5b0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Sun, 11 Mar 2018 17:45:52 +0100 Subject: [PATCH 06/25] Chebyshev direct evaluation --- pygsp/filters/approximations.py | 45 +++++++++++++++++++++++---------- pygsp/filters/filter.py | 10 +++++--- pygsp/tests/test_filters.py | 18 ++++++++----- 3 files changed, 49 insertions(+), 24 deletions(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index 5051f48c..7e701cae 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -21,6 +21,11 @@ class Chebyshev(Filter): Math to show how the coefficients are computed + Evaluation methods (which can be passed when calling :meth:`Filter.evaluate` or :meth:`Filter.filter` are: + + * recursive, defined + * direct, which returns :math:`\sum_k c_k T_k(x) s = \sum_k c_k \cos(k \arccos x) s`. + Parameters ---------- G : graph @@ -45,7 +50,7 @@ def __init__(self, G, filters, order=30): self._coefficients = np.asarray(filters) self.Nf = self._coefficients.shape[1] - def _evaluate(self, x, method='recursive'): + def _evaluate(self, x, method): if x.min() < 0 or x.max() > self.G.lmax: _logger.warning('You are trying to evaluate Chebyshev ' @@ -54,16 +59,18 @@ def _evaluate(self, x, method='recursive'): x = 2 * x / self.G.lmax - 1 # [0, lmax] => [-1, 1] - return self._evaluate_polynomials(x) + method = 'recursive' if method is None else method + + return getattr(self, '_evaluate_' + method)(x) - def _filter(self, s, method='recursive', _=None): + def _filter(self, s, method, _): # method = 'clenshaw' in constructor or filter? M, M = L.shape I = sparse.identity(M, format='csr', dtype=L.dtype) L = 2 * L - self.G.lmax / 2 - I - return self._apply_polynomials(L, s) + return self._evaluate_direct(L, s) def _compute_coefficients(self, filters): r"""Compute the coefficients of the Chebyshev series approximating the filters. @@ -72,7 +79,17 @@ def _compute_coefficients(self, filters): """ pass - def _evaluate_polynomials(self, y, s=1): + def _evaluate_direct(self, x, s=1): + K, F = self._coefficients.shape + c = self._coefficients + c = c.reshape(c.shape + (1,) * x.ndim) + result = np.zeros((F, x.shape[0])) + x_arccos = np.arccos(x) + for k in range(K): + result += c[k] * np.cos(k * x_arccos).dot(s) + return result + + def _evaluate_recursive(self, x, s=1): """Evaluate a Chebyshev series for y. Optionally, times s. .. math: p(y) = \sum_{k=0}^{K} a_k * T_k(y) * s @@ -94,18 +111,18 @@ def _evaluate_polynomials(self, y, s=1): """ - K = self._coefficients.shape[1] + K = self._coefficients.shape[0] c = self._coefficients # Reshaping the coefficients to use broadcasting. - c = c.reshape(c.shape + (1,) * y.ndim) + c = c.reshape(c.shape + (1,) * x.ndim) - y0 = np.ones_like(y) - result = c[0] * y0.dot(s) + x0 = np.ones_like(x) + result = c[0] * x0.dot(s) if K > 1: - y1 = y - result += c[1] * y1 + x1 = x + result += c[1] * x1 for k in range(2, K): - y2 = 2 * y.dot(s) * y1 - y0 - result += c[k] * y2 - y0, y1 = y1, y2 + x2 = 2 * x.dot(s) * x1 - x0 + result += c[k] * x2 + x0, x1 = x1, x2 return result diff --git a/pygsp/filters/filter.py b/pygsp/filters/filter.py index c619e22c..f6df0e75 100644 --- a/pygsp/filters/filter.py +++ b/pygsp/filters/filter.py @@ -60,7 +60,7 @@ def __init__(self, G, kernels): self.Nf = len(kernels) - def evaluate(self, x): + def evaluate(self, x, method=None): r"""Evaluate the kernels at given frequencies. Parameters @@ -88,9 +88,9 @@ def evaluate(self, x): """ x = np.asarray(x) - return self._evaluate(x) + return self._evaluate(x, method) - def _evaluate(self, x): + 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, len(x))) @@ -98,7 +98,7 @@ def _evaluate(self, x): y[i] = kernel(x) return y - def filter(self, s, method='chebyshev', order=30): + 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, @@ -225,6 +225,8 @@ def filter(self, s, method='chebyshev', order=30): 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)) diff --git a/pygsp/tests/test_filters.py b/pygsp/tests/test_filters.py index 3af1244c..72b62a60 100644 --- a/pygsp/tests/test_filters.py +++ b/pygsp/tests/test_filters.py @@ -239,17 +239,23 @@ def setUpClass(cls): cls._rs = np.random.RandomState(42) cls._signal = cls._rs.uniform(size=cls._G.N) - def test_chebyshev_basis(self): + def test_chebyshev_basis(self, K=5, c=2, N=100): r""" Test that the evaluation of the Chebyshev series yields the expected - basis. + basis. We only test the first two elements here. The higher-order + polynomials are compared with the trigonometric definition. """ - K = 5 - c = 2 f = filters.Chebyshev(self._G, c * np.identity(K)) - N = 100 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, 100)) - # TODO: do it for higher orders with the trigonometric definition + + def test_evaluation_methods(self, K=30, F=5, N=100): + r"""Test that all evaluation methods return the same results.""" + coefficients = np.random.RandomState(42).uniform(size=(K, 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') + np.testing.assert_allclose(y1, y2) From 5e563a2337e7c8cbf0d8d806245d39f6307da1f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Sun, 11 Mar 2018 18:33:10 +0100 Subject: [PATCH 07/25] evaluate on n-dimensional arrays --- pygsp/filters/approximations.py | 3 ++- pygsp/tests/test_filters.py | 7 ++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index 7e701cae..77762f59 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -59,6 +59,7 @@ def _evaluate(self, x, method): x = 2 * x / self.G.lmax - 1 # [0, lmax] => [-1, 1] + # The recursive method is the fastest. method = 'recursive' if method is None else method return getattr(self, '_evaluate_' + method)(x) @@ -83,7 +84,7 @@ def _evaluate_direct(self, x, s=1): K, F = self._coefficients.shape c = self._coefficients c = c.reshape(c.shape + (1,) * x.ndim) - result = np.zeros((F, x.shape[0])) + result = np.zeros((F,) + x.shape) x_arccos = np.arccos(x) for k in range(K): result += c[k] * np.cos(k * x_arccos).dot(s) diff --git a/pygsp/tests/test_filters.py b/pygsp/tests/test_filters.py index 72b62a60..09638fcc 100644 --- a/pygsp/tests/test_filters.py +++ b/pygsp/tests/test_filters.py @@ -249,7 +249,7 @@ def test_chebyshev_basis(self, K=5, c=2, N=100): 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, 100)) + 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.""" @@ -259,3 +259,8 @@ def test_evaluation_methods(self, K=30, F=5, N=100): y1 = f.evaluate(x, method='recursive') y2 = f.evaluate(x, method='direct') np.testing.assert_allclose(y1, y2) + # Evaluate on n-dimensional arrays. + x = np.random.RandomState(42).uniform(0, self._G.lmax, size=(3, 1, 19)) + y1 = f.evaluate(x, method='recursive') + y2 = f.evaluate(x, method='direct') + np.testing.assert_allclose(y1, y2) From 9f576ddd8a46015b8b6dae03fb59f67b19a14914 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Wed, 14 Mar 2018 08:57:16 +0100 Subject: [PATCH 08/25] first implementation of Chebyshev recursive filter --- pygsp/filters/approximations.py | 50 ++++++++++++++++++++++++--------- pygsp/filters/filter.py | 6 ++-- pygsp/tests/test_filters.py | 13 +++++++++ 3 files changed, 53 insertions(+), 16 deletions(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index 77762f59..448e436f 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -55,7 +55,7 @@ def _evaluate(self, x, method): if x.min() < 0 or x.max() > self.G.lmax: _logger.warning('You are trying to evaluate Chebyshev ' 'polynomials outside of their orthonormal ' - 'domain [0, {:.2f}].'.format(self.G.lmax)) + 'domain [0, G.lmax={:.2f}].'.format(self.G.lmax)) x = 2 * x / self.G.lmax - 1 # [0, lmax] => [-1, 1] @@ -67,11 +67,18 @@ def _evaluate(self, x, method): def _filter(self, s, method, _): # method = 'clenshaw' in constructor or filter? - M, M = L.shape - I = sparse.identity(M, format='csr', dtype=L.dtype) - L = 2 * L - self.G.lmax / 2 - I + L = self.G.L + if not sparse.issparse(L): + I = np.identity(self.G.N, dtype=L.dtype) + else: + I = sparse.identity(self.G.N, format=L.format, dtype=L.dtype) - return self._evaluate_direct(L, s) + L = 2 * L / self.G.lmax - I # [0, lmax] => [-1, 1] + + # The recursive method is the fastest. + method = 'recursive' if method is None else method + + return getattr(self, '_filter_' + method)(L, s) def _compute_coefficients(self, filters): r"""Compute the coefficients of the Chebyshev series approximating the filters. @@ -80,17 +87,17 @@ def _compute_coefficients(self, filters): """ pass - def _evaluate_direct(self, x, s=1): + def _evaluate_direct(self, x): K, F = self._coefficients.shape c = self._coefficients c = c.reshape(c.shape + (1,) * x.ndim) result = np.zeros((F,) + x.shape) - x_arccos = np.arccos(x) + x = np.arccos(x) for k in range(K): - result += c[k] * np.cos(k * x_arccos).dot(s) + result += c[k] * np.cos(k * x) return result - def _evaluate_recursive(self, x, s=1): + def _evaluate_recursive(self, x): """Evaluate a Chebyshev series for y. Optionally, times s. .. math: p(y) = \sum_{k=0}^{K} a_k * T_k(y) * s @@ -112,18 +119,33 @@ def _evaluate_recursive(self, x, s=1): """ - K = self._coefficients.shape[0] c = self._coefficients - # Reshaping the coefficients to use broadcasting. - c = c.reshape(c.shape + (1,) * x.ndim) + K = c.shape[0] + c = c.reshape(c.shape + (1,) * x.ndim) # For broadcasting. x0 = np.ones_like(x) - result = c[0] * x0.dot(s) + result = c[0] * x0 if K > 1: x1 = x result += c[1] * x1 for k in range(2, K): - x2 = 2 * x.dot(s) * x1 - x0 + x2 = 2 * x * x1 - x0 + result += c[k] * x2 + x0, x1 = x1, x2 + return result + + def _filter_recursive(self, L, s): + c = self._coefficients + K = c.shape[0] + c = c.reshape(c.shape + (1,) * s.ndim) # For broadcasting. + + x0 = s + result = c[0] * x0 + if K > 1: + x1 = L.dot(s) + result += c[1] * x1 + for k in range(2, K): + x2 = 2 * L.dot(x1) - x0 result += c[k] * x2 x0, x1 = x1, x2 return result diff --git a/pygsp/filters/filter.py b/pygsp/filters/filter.py index f6df0e75..69898203 100644 --- a/pygsp/filters/filter.py +++ b/pygsp/filters/filter.py @@ -1,6 +1,7 @@ # -*- coding: utf-8 -*- import numpy as np +from scipy import sparse from pygsp import utils # prevent circular import in Python < 3.5 @@ -87,7 +88,7 @@ def evaluate(self, x, method=None): [] """ - x = np.asarray(x) + x = np.asarray(x) # For iterables. Sparse makes no sense here. return self._evaluate(x, method) def _evaluate(self, x, _): @@ -219,7 +220,8 @@ def filter(self, s, method=None, order=30): True """ - s = np.asarray(s) + if not sparse.issparse(s): + s = np.asarray(s) # For iterables. return self._filter(s, method, order) def _filter(self, s, method='chebyshev', order=30): diff --git a/pygsp/tests/test_filters.py b/pygsp/tests/test_filters.py index 09638fcc..2972738c 100644 --- a/pygsp/tests/test_filters.py +++ b/pygsp/tests/test_filters.py @@ -264,3 +264,16 @@ def test_evaluation_methods(self, K=30, F=5, N=100): y1 = f.evaluate(x, method='recursive') y2 = f.evaluate(x, method='direct') np.testing.assert_allclose(y1, y2) + + def test_filter_identity(self, c=2.3): + r"""Test that filtering with c0 only scales the signal.""" + x = np.random.uniform(size=(self._G.N, 2, 3)) + f = filters.Chebyshev(self._G, [[c]]) + y = f.filter(x, method='recursive').squeeze() + 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').squeeze() + self._G.L = L + np.testing.assert_equal(y, c * x) From 8e36852c27900286c5389c2a933eee2b8c62be94 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Thu, 15 Mar 2018 18:55:29 +0100 Subject: [PATCH 09/25] Chebyshev: 3D signals and 2D coefficients --- pygsp/filters/approximations.py | 54 ++++++++++++++++++++++----------- pygsp/tests/test_filters.py | 22 +++++++++----- 2 files changed, 52 insertions(+), 24 deletions(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index 448e436f..9290d1d6 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -33,6 +33,8 @@ class Chebyshev(Filter): 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 Polynomial order. @@ -48,6 +50,8 @@ def __init__(self, G, filters, order=30): self.Nf = filters.Nf except: self._coefficients = np.asarray(filters) + while self._coefficients.ndim < 3: + self._coefficients = np.expand_dims(self._coefficients, -1) self.Nf = self._coefficients.shape[1] def _evaluate(self, x, method): @@ -66,6 +70,8 @@ def _evaluate(self, x, method): def _filter(self, s, method, _): # method = 'clenshaw' in constructor or filter? + # Might be faster with signals in fortran-contiguous format. + # s: N_SIGNALS x N_FEATURES x N_NODES L = self.G.L if not sparse.issparse(L): @@ -88,14 +94,15 @@ def _compute_coefficients(self, filters): pass def _evaluate_direct(self, x): - K, F = self._coefficients.shape + r"""Evaluate Fout*Fin polynomials at each value in x.""" c = self._coefficients - c = c.reshape(c.shape + (1,) * x.ndim) - result = np.zeros((F,) + x.shape) + K, Fout, Fin = c.shape # #order x #features_out x #features_in + c = c.reshape((K, Fout * Fin) + (1,) * x.ndim) # For broadcasting. + result = np.zeros((Fout * Fin,) + x.shape) x = np.arccos(x) for k in range(K): result += c[k] * np.cos(k * x) - return result + return result.reshape((Fout, Fin) + x.shape).squeeze() def _evaluate_recursive(self, x): """Evaluate a Chebyshev series for y. Optionally, times s. @@ -120,8 +127,8 @@ def _evaluate_recursive(self, x): """ c = self._coefficients - K = c.shape[0] - c = c.reshape(c.shape + (1,) * x.ndim) # For broadcasting. + K, Fout, Fin = c.shape # #order x #features_out x #features_in + c = c.reshape((K, Fout * Fin) + (1,) * x.ndim) # For broadcasting. x0 = np.ones_like(x) result = c[0] * x0 @@ -132,20 +139,33 @@ def _evaluate_recursive(self, x): x2 = 2 * x * x1 - x0 result += c[k] * x2 x0, x1 = x1, x2 - return result + return result.reshape((Fout, Fin) + x.shape).squeeze() 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 = c.shape[0] - c = c.reshape(c.shape + (1,) * s.ndim) # For broadcasting. - - x0 = s - result = c[0] * x0 + 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 = (N, Fin, M) + return np.einsum('oi,nim->nom', c, x) + def dot(L, x): + """One diffusion step by multiplication with the Laplacian.""" + x.shape = (N, Fin * M) + return L.dot(x) + + x0 = np.asarray(s.T, order='C') + result = mult(c[0], x0) if K > 1: - x1 = L.dot(s) - result += c[1] * x1 + x1 = dot(L, x0) + result += mult(c[1], x1) for k in range(2, K): - x2 = 2 * L.dot(x1) - x0 - result += c[k] * x2 + x2 = 2 * dot(L, x1) - x0 + result += mult(c[k], x2) x0, x1 = x1, x2 - return result + return result.T diff --git a/pygsp/tests/test_filters.py b/pygsp/tests/test_filters.py index 2972738c..091e9e93 100644 --- a/pygsp/tests/test_filters.py +++ b/pygsp/tests/test_filters.py @@ -253,27 +253,35 @@ def test_chebyshev_basis(self, K=5, c=2, N=100): def test_evaluation_methods(self, K=30, F=5, N=100): r"""Test that all evaluation methods return the same results.""" - coefficients = np.random.RandomState(42).uniform(size=(K, F)) + 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') np.testing.assert_allclose(y1, y2) # Evaluate on n-dimensional arrays. - x = np.random.RandomState(42).uniform(0, self._G.lmax, size=(3, 1, 19)) + x = self._rs.uniform(0, self._G.lmax, size=(3, 1, 19)) y1 = f.evaluate(x, method='recursive') y2 = f.evaluate(x, method='direct') np.testing.assert_allclose(y1, y2) - def test_filter_identity(self, c=2.3): + def test_filter_identity(self, M=10, c=2.3): r"""Test that filtering with c0 only scales the signal.""" - x = np.random.uniform(size=(self._G.N, 2, 3)) - f = filters.Chebyshev(self._G, [[c]]) - y = f.filter(x, method='recursive').squeeze() + x = self._rs.uniform(size=(M, 1, 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').squeeze() + 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) + y = f.filter(x) + self.assertTupleEqual(y.shape, (M, Fout, self._G.N)) From ccd590b48022709d550c067d445b6ed127202d58 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Thu, 15 Mar 2018 21:30:40 +0100 Subject: [PATCH 10/25] Chebyshev: Clenshaw evaluation --- pygsp/filters/approximations.py | 79 ++++++++++++++++++++++----------- 1 file changed, 54 insertions(+), 25 deletions(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index 9290d1d6..9821b92e 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -63,15 +63,18 @@ def _evaluate(self, x, method): x = 2 * x / self.G.lmax - 1 # [0, lmax] => [-1, 1] - # The recursive method is the fastest. + 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. + + # Recursive faster than direct faster than clenshaw. method = 'recursive' if method is None else method - return getattr(self, '_evaluate_' + method)(x) + y = getattr(self, '_evaluate_' + method)(c, x) + + return y.reshape((Fout, Fin) + x.shape).squeeze() def _filter(self, s, method, _): - # method = 'clenshaw' in constructor or filter? - # Might be faster with signals in fortran-contiguous format. - # s: N_SIGNALS x N_FEATURES x N_NODES L = self.G.L if not sparse.issparse(L): @@ -81,7 +84,7 @@ def _filter(self, s, method, _): L = 2 * L / self.G.lmax - I # [0, lmax] => [-1, 1] - # The recursive method is the fastest. + # Recursive and clenshaw are similarly fast. method = 'recursive' if method is None else method return getattr(self, '_filter_' + method)(L, s) @@ -93,18 +96,16 @@ def _compute_coefficients(self, filters): """ pass - def _evaluate_direct(self, x): + def _evaluate_direct(self, c, x): r"""Evaluate Fout*Fin polynomials at each value in x.""" - 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. - result = np.zeros((Fout * Fin,) + x.shape) + 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.reshape((Fout, Fin) + x.shape).squeeze() + return result - def _evaluate_recursive(self, x): + def _evaluate_recursive(self, c, x): """Evaluate a Chebyshev series for y. Optionally, times s. .. math: p(y) = \sum_{k=0}^{K} a_k * T_k(y) * s @@ -125,11 +126,7 @@ def _evaluate_recursive(self, x): corresponding Chebyshev Series. (size F x N) """ - - 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. - + K = c.shape[0] x0 = np.ones_like(x) result = c[0] * x0 if K > 1: @@ -139,7 +136,7 @@ def _evaluate_recursive(self, x): x2 = 2 * x * x1 - x0 result += c[k] * x2 x0, x1 = x1, x2 - return result.reshape((Fout, Fin) + x.shape).squeeze() + return result def _filter_recursive(self, L, s): r"""Filter a signal with the 3-way recursive relation. @@ -152,14 +149,14 @@ def _filter_recursive(self, L, s): def mult(c, x): """Multiply the diffused signals by the Chebyshev coefficients.""" - x.shape = (N, Fin, M) - return np.einsum('oi,nim->nom', c, x) + 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 = (N, Fin * M) - return L.dot(x) + x.shape = (M * Fin, N) + return L.__rmul__(x) # x @ L - x0 = np.asarray(s.T, order='C') + x0 = s.view() result = mult(c[0], x0) if K > 1: x1 = dot(L, x0) @@ -168,4 +165,36 @@ def dot(L, x): x2 = 2 * dot(L, x1) - x0 result += mult(c[k], x2) x0, x1 = x1, x2 - return result.T + 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.__rmul__(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 From 65346bfd1fa8c5c5e64eacb376b481f2403a3400 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Fri, 16 Mar 2018 10:58:29 +0100 Subject: [PATCH 11/25] catch unknown filter or evaluate method --- pygsp/filters/approximations.py | 10 ++++++++-- pygsp/tests/test_filters.py | 14 ++++++++++++-- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index 9821b92e..4774c2a4 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -70,7 +70,10 @@ def _evaluate(self, x, method): # Recursive faster than direct faster than clenshaw. method = 'recursive' if method is None else method - y = getattr(self, '_evaluate_' + method)(c, x) + try: + y = getattr(self, '_evaluate_' + method)(c, x) + except AttributeError: + raise ValueError('Unknown method {}.'.format(method)) return y.reshape((Fout, Fin) + x.shape).squeeze() @@ -87,7 +90,10 @@ def _filter(self, s, method, _): # Recursive and clenshaw are similarly fast. method = 'recursive' if method is None else method - return getattr(self, '_filter_' + method)(L, s) + try: + return getattr(self, '_filter_' + method)(L, s) + except AttributeError: + raise ValueError('Unknown method {}.'.format(method)) def _compute_coefficients(self, filters): r"""Compute the coefficients of the Chebyshev series approximating the filters. diff --git a/pygsp/tests/test_filters.py b/pygsp/tests/test_filters.py index 091e9e93..c7220b4f 100644 --- a/pygsp/tests/test_filters.py +++ b/pygsp/tests/test_filters.py @@ -258,12 +258,18 @@ def test_evaluation_methods(self, K=30, F=5, N=100): 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.""" @@ -283,5 +289,9 @@ def test_filter_methods(self, K=30, Fin=5, Fout=6, M=100): coefficients = self._rs.uniform(size=(K, Fout, Fin)) x = self._rs.uniform(size=(M, Fin, self._G.N)) f = filters.Chebyshev(self._G, coefficients) - y = f.filter(x) - self.assertTupleEqual(y.shape, (M, Fout, self._G.N)) + 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) + # Unknown method. + self.assertRaises(ValueError, f.filter, x, method='unk') From 6d45b02843992369758db08749a5db77a9192631 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Fri, 16 Mar 2018 12:54:00 +0100 Subject: [PATCH 12/25] scaling functions shall be callable from outside --- pygsp/filters/approximations.py | 33 ++++++++++++++++++++++----------- pygsp/tests/test_filters.py | 7 +++++++ 2 files changed, 29 insertions(+), 11 deletions(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index 4774c2a4..c9788d73 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -54,14 +54,31 @@ def __init__(self, G, filters, order=30): self._coefficients = np.expand_dims(self._coefficients, -1) self.Nf = self._coefficients.shape[1] - def _evaluate(self, x, method): - if x.min() < 0 or x.max() > self.G.lmax: + @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, G.lmax={:.2f}].'.format(self.G.lmax)) + '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) - x = 2 * x / self.G.lmax - 1 # [0, lmax] => [-1, 1] + return 2 * L / lmax - I + + def _evaluate(self, x, method): + + x = self.scale_data(x, self.G.lmax) c = self._coefficients K, Fout, Fin = c.shape # #order x #features_out x #features_in @@ -79,13 +96,7 @@ def _evaluate(self, x, method): def _filter(self, s, method, _): - L = self.G.L - if not sparse.issparse(L): - I = np.identity(self.G.N, dtype=L.dtype) - else: - I = sparse.identity(self.G.N, format=L.format, dtype=L.dtype) - - L = 2 * L / self.G.lmax - I # [0, lmax] => [-1, 1] + L = self.scale_operator(self.G.L, self.G.lmax) # Recursive and clenshaw are similarly fast. method = 'recursive' if method is None else method diff --git a/pygsp/tests/test_filters.py b/pygsp/tests/test_filters.py index c7220b4f..64b0c241 100644 --- a/pygsp/tests/test_filters.py +++ b/pygsp/tests/test_filters.py @@ -239,6 +239,13 @@ def setUpClass(cls): 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 From cb6cbd1f5ed2c8afc7052a193a7fc8a5d8fbca12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Sun, 18 Mar 2018 12:03:01 +0100 Subject: [PATCH 13/25] Chebyshev: coefficients computation --- pygsp/filters/approximations.py | 100 +++++++++++++++++++++++++------- pygsp/filters/filter.py | 21 +++++++ pygsp/tests/test_filters.py | 31 +++++++++- 3 files changed, 130 insertions(+), 22 deletions(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index c9788d73..0af87c13 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -40,20 +40,85 @@ class Chebyshev(Filter): """ - def __init__(self, G, filters, order=30): + def __init__(self, G, coefficients): self.G = G - self.order = order - try: - self._compute_coefficients(filters) - self.Nf = filters.Nf - except: - self._coefficients = np.asarray(filters) - while self._coefficients.ndim < 3: - self._coefficients = np.expand_dims(self._coefficients, -1) - self.Nf = self._coefficients.shape[1] + coefficients = np.asarray(coefficients) + while coefficients.ndim < 3: + coefficients = np.expand_dims(coefficients, -1) + + self.Nout, self.Nin = coefficients.shape[1:] + self.Nf = self.Nout * self.Nin + 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 + -------- + >>> G = graphs.Ring(50) + >>> 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 + + """ + 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): @@ -106,13 +171,6 @@ def _filter(self, s, method, _): except AttributeError: raise ValueError('Unknown method {}.'.format(method)) - def _compute_coefficients(self, filters): - r"""Compute the coefficients of the Chebyshev series approximating the filters. - - Some implementations define c_0 / 2. - """ - pass - def _evaluate_direct(self, c, x): r"""Evaluate Fout*Fin polynomials at each value in x.""" K, F = c.shape[:2] @@ -125,13 +183,13 @@ def _evaluate_direct(self, c, x): def _evaluate_recursive(self, c, x): """Evaluate a Chebyshev series for y. Optionally, times s. - .. math: p(y) = \sum_{k=0}^{K} a_k * T_k(y) * 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) - y: array-like + x: array-like vector to be evaluated. (size N x 1) vector or matrix signal: array-like @@ -171,7 +229,7 @@ def mult(c, x): def dot(L, x): """One diffusion step by multiplication with the Laplacian.""" x.shape = (M * Fin, N) - return L.__rmul__(x) # x @ L + return L.__rmatmul__(x) # x @ L x0 = s.view() result = mult(c[0], x0) @@ -195,7 +253,7 @@ def mult(c, s): def dot(L, x): """One diffusion step by multiplication with the Laplacian.""" x.shape = (M * Fout, N) - y = L.__rmul__(x) # x @ L + y = L.__rmatmul__(x) # x @ L x.shape = (M, Fout, N) y.shape = (M, Fout, N) return y diff --git a/pygsp/filters/filter.py b/pygsp/filters/filter.py index 69898203..c92a15e2 100644 --- a/pygsp/filters/filter.py +++ b/pygsp/filters/filter.py @@ -99,6 +99,27 @@ def _evaluate(self, x, _): y[i] = kernel(x) return y + def approximate(self, method, **kwargs): + r"""Returns 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 + -------- + TODO: approx plot from notebook (needs new plotting) + + """ + 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). diff --git a/pygsp/tests/test_filters.py b/pygsp/tests/test_filters.py index 64b0c241..9e16c799 100644 --- a/pygsp/tests/test_filters.py +++ b/pygsp/tests/test_filters.py @@ -299,6 +299,35 @@ def test_filter_methods(self, K=30, Fin=5, Fout=6, M=100): 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) + 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) From c9a277ae4ee70c964226ad57d22319cfff0314d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Mon, 2 Apr 2018 16:33:04 +0200 Subject: [PATCH 14/25] doc: examples of filter approximations --- pygsp/filters/__init__.py | 1 + pygsp/filters/approximations.py | 52 ++++++++++++++++++++++++++++++++- pygsp/filters/filter.py | 32 ++++++++++++++++++-- 3 files changed, 81 insertions(+), 4 deletions(-) diff --git a/pygsp/filters/__init__.py b/pygsp/filters/__init__.py index 39726683..35d628e4 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 diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index 0af87c13..8336fd63 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -38,6 +38,41 @@ class Chebyshev(Filter): order : int 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 @ 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.]] + """ def __init__(self, G, coefficients): @@ -92,13 +127,28 @@ def from_filter(cls, filters, order=30, n=None): Examples -------- - >>> G = graphs.Ring(50) + + 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 diff --git a/pygsp/filters/filter.py b/pygsp/filters/filter.py index 376402b9..c30e744a 100644 --- a/pygsp/filters/filter.py +++ b/pygsp/filters/filter.py @@ -100,7 +100,7 @@ def _evaluate(self, x, _): return y def approximate(self, method, **kwargs): - r"""Returns a filter which approximates this filter. + r"""Return a filter which approximates this filter. While approximations might loose accuracy, they allow for much faster computations. @@ -114,7 +114,33 @@ def approximate(self, method, **kwargs): Examples -------- - TODO: approx plot from notebook (needs new plotting) + + 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 @@ -173,7 +199,7 @@ def filter(self, s, method=None, 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). From dd44dba33fc18d3c1e752a6665bc7f419d189354 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Tue, 3 Apr 2018 10:50:18 +0200 Subject: [PATCH 15/25] allow user to choose the tolerance when estimating lmax --- pygsp/graphs/graph.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/pygsp/graphs/graph.py b/pygsp/graphs/graph.py index 7e4633b3..b764e620 100644 --- a/pygsp/graphs/graph.py +++ b/pygsp/graphs/graph.py @@ -540,7 +540,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. @@ -551,6 +551,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. @@ -571,21 +575,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. ' From cc9f508899049a74cf6e7c55944321edf1f2edec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Tue, 3 Apr 2018 16:02:18 +0200 Subject: [PATCH 16/25] localize and plot 2D (n to m features) filterbanks --- pygsp/filters/approximations.py | 6 ++++-- pygsp/filters/filter.py | 14 +++++++++++--- pygsp/plotting.py | 4 +++- 3 files changed, 18 insertions(+), 6 deletions(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index 8336fd63..2b20fe0b 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -83,8 +83,10 @@ def __init__(self, G, coefficients): while coefficients.ndim < 3: coefficients = np.expand_dims(coefficients, -1) - self.Nout, self.Nin = coefficients.shape[1:] - self.Nf = self.Nout * self.Nin + 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. diff --git a/pygsp/filters/filter.py b/pygsp/filters/filter.py index c30e744a..1a29c488 100644 --- a/pygsp/filters/filter.py +++ b/pygsp/filters/filter.py @@ -59,7 +59,11 @@ def __init__(self, G, kernels): kernels = [kernels] self._kernels = kernels - self.Nf = len(kernels) + # 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. def evaluate(self, x, method=None): r"""Evaluate the kernels at given frequencies. @@ -398,8 +402,12 @@ def localize(self, i, **kwargs): >>> G.plot_signal(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/plotting.py b/pygsp/plotting.py index ccb228bc..f0cff5c9 100644 --- a/pygsp/plotting.py +++ b/pygsp/plotting.py @@ -369,7 +369,9 @@ def _plot_filter(filters, n, eigenvalues, sum, ax, **kwargs): ax.axvline(x=e, color=[0.9]*3, linewidth=1) x = np.linspace(0, G.lmax, n) - y = filters.evaluate(x).T + # Plot all filters on one figure. A user can plot a single filterbank with + # filter[1].plot() or a single filter with filter[1, 3].plot(). + y = filters.evaluate(x).reshape((-1, n)).T ax.plot(x, y, **kwargs) # TODO: plot highlighted eigenvalues From 5ca1ae05387afe36de62973767ae89c61d0d0964 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Tue, 3 Apr 2018 16:07:07 +0200 Subject: [PATCH 17/25] normalize signal's shapes to #signal x #features x #nodes --- pygsp/filters/approximations.py | 35 +++++++++++++++++++++ pygsp/filters/filter.py | 9 ++++-- pygsp/tests/test_filters.py | 55 ++++++++++++++++++++++++++++++++- 3 files changed, 95 insertions(+), 4 deletions(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index 2b20fe0b..7f05dd78 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -213,6 +213,41 @@ def _evaluate(self, x, method): def _filter(self, s, method, _): + # 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)) + + # 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: + # 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. diff --git a/pygsp/filters/filter.py b/pygsp/filters/filter.py index 1a29c488..bddf4212 100644 --- a/pygsp/filters/filter.py +++ b/pygsp/filters/filter.py @@ -273,7 +273,11 @@ def filter(self, s, method=None, order=30): """ if not sparse.issparse(s): s = np.asarray(s) # For iterables. - return self._filter(s, method, order) + + 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.""" @@ -344,8 +348,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`.""" diff --git a/pygsp/tests/test_filters.py b/pygsp/tests/test_filters.py index 9e16c799..0496ffe2 100644 --- a/pygsp/tests/test_filters.py +++ b/pygsp/tests/test_filters.py @@ -280,7 +280,7 @@ def test_evaluation_methods(self, K=30, F=5, N=100): 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, 1, self._G.N)) + 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) @@ -331,3 +331,56 @@ def test_approximations(self, N=100, K=20): 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) From ab8ff760b8463c6808b9158e3b6609157e48a0f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Tue, 3 Apr 2018 16:46:53 +0200 Subject: [PATCH 18/25] invalid syntax for python 2.7 and 3.4 --- pygsp/filters/approximations.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index 7f05dd78..64031980 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -66,7 +66,7 @@ class Chebyshev(Filter): They are orthogonal in the vertex domain: >>> s = s.T.reshape((G.N, -1)) - >>> print(s.T @ s) + >>> print(s.T.dot(s)) [[20. 0. 0. 0. 0.] [ 0. 10. 0. 0. 0.] [ 0. 0. 10. 0. 0.] From f0b66a37197f8699b0a462243fdf467a32ae14f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Wed, 4 Apr 2018 00:57:47 +0200 Subject: [PATCH 19/25] much faster graph construction (avoid constructing matrices) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Benchmark on the whole test suite (i.e. make test): * Before: Ran 140 tests in 336.127s * Avoid building LIL and TRIL: Ran 140 tests in 310.229s * Avoid building the CSR if not necessary: Ran 140 tests in 303.354s Benchmark pygsp.graphs.Graph(W) on large graph (800k nodes, 3M edges): * Before: 14.7 s ± 693 ms * After: 207 ms ± 15.4 ms --- doc/tutorials/intro.rst | 2 +- pygsp/graphs/graph.py | 21 ++++++++++++++++----- 2 files changed, 17 insertions(+), 6 deletions(-) diff --git a/doc/tutorials/intro.rst b/doc/tutorials/intro.rst index d8b09e37..a3d9b654 100644 --- a/doc/tutorials/intro.rst +++ b/doc/tutorials/intro.rst @@ -68,7 +68,7 @@ We can retrieve our weight matrix, which is stored in a sparse format. >>> (G.W == W).all() True >>> type(G.W) - + We can access the `graph Laplacian`_ diff --git a/pygsp/graphs/graph.py b/pygsp/graphs/graph.py index b764e620..ef99e9d0 100644 --- a/pygsp/graphs/graph.py +++ b/pygsp/graphs/graph.py @@ -71,20 +71,31 @@ def __init__(self, W, gtype='unknown', lap_type='combinatorial', if len(W.shape) != 2 or W.shape[0] != W.shape[1]: raise ValueError('W has incorrect shape {}'.format(W.shape)) + # CSR sparse matrices are the most efficient for matrix multiplication. + # They are the sole sparse matrix type to support eliminate_zeros(). + if sparse.isspmatrix_csr(W): + self.W = W + else: + self.W = sparse.csr_matrix(W) + # Don't keep edges of 0 weight. Otherwise Ne will not correspond to the # real number of edges. Problematic when e.g. plotting. - W = sparse.csr_matrix(W) - W.eliminate_zeros() + self.W.eliminate_zeros() self.N = W.shape[0] - self.W = sparse.lil_matrix(W) + + # TODO: why would we ever want this? + # For large matrices it slows the graph construction by a factor 100. + # self.W = sparse.lil_matrix(self.W) # Don't count edges two times if undirected. # Be consistent with the size of the differential operator. if self.is_directed(): self.Ne = self.W.nnz else: - self.Ne = sparse.tril(W).nnz + diagonal = np.count_nonzero(self.W.diagonal()) + off_diagonal = self.W.nnz - diagonal + self.Ne = off_diagonal // 2 + diagonal self.check_weights() @@ -634,7 +645,7 @@ def get_edge_list(self): else: v_in, v_out = sparse.tril(self.W).nonzero() weights = self.W[v_in, v_out] - weights = weights.toarray().squeeze() + weights = np.asarray(weights).squeeze() # TODO G.ind_edges = sub2ind(size(G.W), G.v_in, G.v_out) From 717ceb58de936c2973f4447694cda56290552b3f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Tue, 10 Apr 2018 22:00:20 +0200 Subject: [PATCH 20/25] bug: wrong error message when n_features_in=1 with invalid 3D signal shape --- pygsp/filters/approximations.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index 64031980..c8c79398 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -226,7 +226,7 @@ def _filter(self, s, method, _): if s.ndim == 1: s = np.expand_dims(s, 0) if s.shape[-2] != self.n_features_in: - if self.n_features_in == 1: + 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: From 9b4de7c940b48966011cbf6d2679c82acf33269b Mon Sep 17 00:00:00 2001 From: Nathanael Perraudin Date: Thu, 5 Jul 2018 11:05:07 +0200 Subject: [PATCH 21/25] merge with master --- pygsp/graphs/graph.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pygsp/graphs/graph.py b/pygsp/graphs/graph.py index a007ed35..9995e69a 100644 --- a/pygsp/graphs/graph.py +++ b/pygsp/graphs/graph.py @@ -88,8 +88,7 @@ def __init__(self, W, lap_type='combinatorial', coords=None, plotting={}): else: diagonal = np.count_nonzero(self.W.diagonal()) off_diagonal = self.W.nnz - diagonal - - self.n_edges = off_diagonal // 2 + diagonal + self.n_edges = off_diagonal // 2 + diagonal self.check_weights() From ac7c5ec59d24eadf07c838b83334b816ac5bf040 Mon Sep 17 00:00:00 2001 From: Nathanael Perraudin Date: Thu, 5 Jul 2018 12:55:01 +0200 Subject: [PATCH 22/25] add the repr back for filter --- pygsp/filters/filter.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/pygsp/filters/filter.py b/pygsp/filters/filter.py index e9139103..64bf99f6 100644 --- a/pygsp/filters/filter.py +++ b/pygsp/filters/filter.py @@ -65,16 +65,16 @@ def __init__(self, G, kernels): self.n_filters = self.n_features_in * self.n_features_out self.Nf = self.n_filters # TODO: kept for backward compatibility only. - # def _get_extra_repr(self): - # return dict() - - # def __repr__(self): - # attrs = {'in': self.n_features_in, 'out': self.n_features_out} - # attrs.update(self._get_extra_repr()) - # s = '' - # for key, value in attrs.items(): - # s += '{}={}, '.format(key, value) - # return '{}({})'.format(self.__class__.__name__, s[:-2]) + def _get_extra_repr(self): + return dict() + + def __repr__(self): + attrs = {'in': self.n_features_in, 'out': self.n_features_out} + attrs.update(self._get_extra_repr()) + s = '' + for key, value in attrs.items(): + s += '{}={}, '.format(key, value) + return '{}({})'.format(self.__class__.__name__, s[:-2]) def evaluate(self, x, method=None): r"""Evaluate the kernels at given frequencies. From 579994a8be2eed33b205b53896037bf32372020f Mon Sep 17 00:00:00 2001 From: Nathanael Perraudin Date: Wed, 25 Jul 2018 18:58:50 +0200 Subject: [PATCH 23/25] fix small bug if n_eigenvector is set --- pygsp/graphs/fourier.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pygsp/graphs/fourier.py b/pygsp/graphs/fourier.py index 938c2b91..55cc9c11 100644 --- a/pygsp/graphs/fourier.py +++ b/pygsp/graphs/fourier.py @@ -141,8 +141,9 @@ def compute_fourier_basis(self, n_eigenvectors=None, recompute=False): assert self._e[-1] <= 2 assert np.max(self._e) == self._e[-1] - self._lmax = self._e[-1] - self._mu = np.max(np.abs(self._U)) + if n_eigenvectors == self.N: + self._lmax = self._e[-1] + self._mu = np.max(np.abs(self._U)) def gft(self, s): r"""Compute the graph Fourier transform. From a12f5a45b821658a624a1d427428c463e3795940 Mon Sep 17 00:00:00 2001 From: Nathanael Perraudin Date: Mon, 10 Dec 2018 11:27:02 +0100 Subject: [PATCH 24/25] fix new test system? --- pygsp/tests/test_all.py | 1 + pygsp/tests/test_filters.py | 3 +++ pygsp/tests/test_graphs.py | 2 ++ pygsp/tests/test_utils.py | 2 ++ 4 files changed, 8 insertions(+) diff --git a/pygsp/tests/test_all.py b/pygsp/tests/test_all.py index 6d42a7c1..bad490b3 100755 --- a/pygsp/tests/test_all.py +++ b/pygsp/tests/test_all.py @@ -20,6 +20,7 @@ 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 2a9e779c..da7b3bdd 100644 --- a/pygsp/tests/test_filters.py +++ b/pygsp/tests/test_filters.py @@ -427,3 +427,6 @@ def test_shape(y, M, Fout, N=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 7b2d9d5d..cd183971 100644 --- a/pygsp/tests/test_graphs.py +++ b/pygsp/tests/test_graphs.py @@ -359,3 +359,5 @@ def test_imgpatches(self): def test_grid2dimgpatches(self): graphs.Grid2dImgPatches(img=self._img, patch_shape=(3, 3)) + +suite = unittest.TestLoader().loadTestsFromTestCase(TestCase) \ No newline at end of file diff --git a/pygsp/tests/test_utils.py b/pygsp/tests/test_utils.py index dc585241..670510c2 100644 --- a/pygsp/tests/test_utils.py +++ b/pygsp/tests/test_utils.py @@ -129,3 +129,5 @@ def test_distanz(x, y): # test_tree_depths(A, root) # test_distanz(x, y) + +suite = unittest.TestLoader().loadTestsFromTestCase(TestCase) \ No newline at end of file From c0868318f26572ebf5702425c82e7b8bce3f77a8 Mon Sep 17 00:00:00 2001 From: Nathanael Perraudin Date: Mon, 10 Dec 2018 12:56:46 +0100 Subject: [PATCH 25/25] fix plots --- pygsp/filters/approximations.py | 6 +++--- pygsp/filters/filter.py | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/pygsp/filters/approximations.py b/pygsp/filters/approximations.py index c8c79398..0f0ebb23 100644 --- a/pygsp/filters/approximations.py +++ b/pygsp/filters/approximations.py @@ -55,8 +55,8 @@ class Chebyshev(Filter): >>> 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]) + >>> _ = 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') @@ -146,7 +146,7 @@ def from_filter(cls, filters, order=30, n=None): >>> 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]) + >>> _ = 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]) diff --git a/pygsp/filters/filter.py b/pygsp/filters/filter.py index 32adc22a..39aa2df6 100644 --- a/pygsp/filters/filter.py +++ b/pygsp/filters/filter.py @@ -142,12 +142,12 @@ def approximate(self, method, **kwargs): >>> G = graphs.Ring() >>> G.compute_fourier_basis() >>> f1 = filters.Heat(G) - >>> f1.plot(eigenvalues=True, linewidth=3, label='continuous', ax=ax) + >>> _ = 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) + ... _ = f2.plot(eigenvalues=False, label=l, linestyle='dashed', ax=ax) >>> >>> _ = ax.set_title('Approximation for various polynomial orders') >>> _ = ax.legend() @@ -158,8 +158,8 @@ def approximate(self, method, **kwargs): >>> G.compute_fourier_basis() >>> f1 = filters.Itersine(G) >>> f2 = f1.approximate('Chebyshev', order=20) - >>> f1.plot(title='Continuous filterbank') - >>> f2.plot(title='Approximated filterbank') + >>> _ = f1.plot(title='Continuous filterbank') + >>> _ = f2.plot(title='Approximated filterbank') """ from . import approximations