From 747aabab0ec536880a225a01be4c6dced90747f0 Mon Sep 17 00:00:00 2001 From: evenmatencio Date: Thu, 30 May 2024 09:42:57 +0200 Subject: [PATCH 01/10] add gfss l2 cost --- src/ruptures/costs/costgfssl2.py | 88 ++++++++++++++++++++++++++++++++ tests/test_costs.py | 25 +++++++++ 2 files changed, 113 insertions(+) create mode 100644 src/ruptures/costs/costgfssl2.py diff --git a/src/ruptures/costs/costgfssl2.py b/src/ruptures/costs/costgfssl2.py new file mode 100644 index 00000000..bfc49cbd --- /dev/null +++ b/src/ruptures/costs/costgfssl2.py @@ -0,0 +1,88 @@ +r"""CostL2 (least squared deviation) after GFSS rotation (low-pass graph +spectral filtering)""" + +import numpy as np + +from scipy.linalg import eigh +from ruptures.costs import NotEnoughPoints +from ruptures.base import BaseCost + + +class CostGFSSL2(BaseCost): + """Applies the GFSS rotation to the whole signal before computing the + standard L2 cost for fixed variance gaussian hypothesis.""" + + model = "gfss_l2_cost" + min_size = 1 + + def __init__(self, laplacian_mat, cut_sparsity) -> None: + """ + Args: + laplacian_mat (array): the discrete Laplacian matrix of the graph: D - W + where D is the diagonal matrix diag(d_i) of the node degrees and W the adjacency matrix + cut_sparsity (float): frequency threshold of the GFSS spectral filter + """ + self.cut_sparsity = cut_sparsity + self.graph_laplacian_mat = laplacian_mat + self.signal = None + self.gft_square_cumsum = None + self.gft_cumsum = None + super().__init__() + + def filter(self, freqs, eps=0.00001): + """Applies the GFSS filter to the input (spatial) frequences. + NOTE: the frequences must be in increasing order. + + Args: + freqs (array): ordered frequences to filter. + eps (float, optional): threshold for non zero values. Defaults to 0.00001. + + Returns: + filtered_freqs (array): the output of the filter. + """ + nb_zeros = np.sum(freqs < eps) + filtered_freqs = np.minimum(1, np.sqrt(self.cut_sparsity / freqs[nb_zeros:])) + return np.concatenate([np.zeros(nb_zeros), filtered_freqs]) + + def fit(self, signal): + """Performs pre-computations for per-segment approximation cost. + + NOTE: the number of dimensions of the signal and their ordering + must match those of the nodes of the graph. + The function eigh used below returns the eigenvector corresponding to + the ith eigenvalue in the ith column eigvect[:, i] + + Args: + signal (array): of shape [n_samples, n_dim]. + """ + self.signal = signal + # Computation of the GFSS + eigvals, eigvects = eigh(self.graph_laplacian_mat) + filter_matrix = np.diag(self.filter(eigvals), k=0) + gft = filter_matrix.dot(eigvects.T.dot(signal.T)).T + # Computation of the per-segment cost utils + self.gft_square_cumsum = np.concatenate( + [np.zeros((1, signal.shape[1])), np.cumsum(gft**2, axis=0)], axis=0 + ) + self.gft_cumsum = np.concatenate( + [np.zeros((1, signal.shape[1])), np.cumsum(gft, axis=0)], axis=0 + ) + return self + + def error(self, start, end): + """Return the L2 approximation cost on the segment [start:end] where + end is excluded, over the filtered signal. + + Args: + start (int): start of the segment + end (int): end of the segment + + Returns: + float: segment cost + """ + if end - start < self.min_size: + raise NotEnoughPoints + + sub_square_sum = self.gft_square_cumsum[end] - self.gft_square_cumsum[start] + sub_sum = self.gft_cumsum[end] - self.gft_cumsum[start] + return np.sum(sub_square_sum - (sub_sum**2) / (end - start)) diff --git a/tests/test_costs.py b/tests/test_costs.py index 6a680d66..e21ade09 100644 --- a/tests/test_costs.py +++ b/tests/test_costs.py @@ -3,6 +3,7 @@ from ruptures import Binseg from ruptures.costs import CostLinear, CostNormal, cost_factory from ruptures.costs.costml import CostMl +from ruptures.costs.costgfssl2 import CostGFSSL2 from ruptures.datasets import pw_constant from ruptures.exceptions import NotEnoughPoints @@ -237,3 +238,27 @@ def test_costl2_small_data(): "got {computed_break_dict}.", ) assert expected_break_dict == computed_break_dict, err_msg + + +def test_costgfssl2(signal_bkps_5D_noisy): + """Test the GFSS cost over a very simple star graph for a signal with one + mean change.""" + signal, bkps = signal_bkps_5D_noisy + n_samples = signal.shape[0] + star_graph_laplacian_mat = np.array( + [ + [4, -1, -1, -1, -1], + [-1, 1, 0, 0, 0], + [-1, 0, 1, 0, 0], + [-1, 0, 0, 1, 0], + [-1, 0, 0, 0, 1], + ] + ) + cut_sparsity = 1 + c = CostGFSSL2(star_graph_laplacian_mat, cut_sparsity).fit(signal) + c.error(0, n_samples // 2) + c.error(100, n_samples) + c.error(10, n_samples // 4) + c.sum_of_costs(bkps) + with pytest.raises(NotEnoughPoints): + c.error(10, 10) From d92e869df6ee08a26bbc95278e717269617a5bd4 Mon Sep 17 00:00:00 2001 From: evenmatencio Date: Thu, 30 May 2024 15:52:19 +0200 Subject: [PATCH 02/10] first version of gfss cost user guide --- docs/user-guide/costs/costgfssl2.md | 85 +++++++++++++++++++++++++++++ 1 file changed, 85 insertions(+) create mode 100644 docs/user-guide/costs/costgfssl2.md diff --git a/docs/user-guide/costs/costgfssl2.md b/docs/user-guide/costs/costgfssl2.md new file mode 100644 index 00000000..ed810b03 --- /dev/null +++ b/docs/user-guide/costs/costgfssl2.md @@ -0,0 +1,85 @@ +# Graph Fourier Scan Statistic least squared deviation (`CostGFSSL2`) + +## Description + +This cost function is specifically designed to detect localized mean-shifts in a graph signal. It relies on the least squared deviation of a graph signal filtered with a low-pass graph spectral filter $h$. + +Formally, let $G = (V, E, W)$ a graph containing $p =|V|$ nodes, with a weighted adjacency matrix $W$. We define its combinatorial Laplacian matrix $L$ as classicaly done in Graph Signal Processing [[Shuman2013](#Shuman2013)]: + +$$ +L = D - W = U \Lambda U^T, +$$ + +where + +- $D$ is the diagional degree matrix of the graph: $D = \text{diag}(d_1, \ldots, d_p)$. +- $U$ is the orthogonal matrix whose columns $\{u_i\}_{i=1}^p$ are the eigenvectors of $L$ +- $\Lambda = \text{diag}(\lambda_1, \ldots, \lambda_p)$ contains the eigenvalues of $L$. + +Let $\{y_t\}_t, ~ y_t \in \mathbb{R}^p$, a multivariate signal on an interval $I$. The Graph Fourier Scan Statistic (GFSS) [[Sharpnack2016](#Sharpnack2016)] of a graph signal $y_t$ is $\|G(y_t)\|_2^2$ with: + +$$ +G(y) = \sum_{i=2}^p h(\lambda_i) (u_i^T y) u_i \quad \text{ and } \quad h(\lambda) = \min \left\{ 1, \sqrt{\frac{\rho}{\lambda}} \right\} +$$ + +where $\rho$ is the so-called cut-sparsity. The cost function [`CostGFSSL2`][ruptures.costs.costrbf.CostGFSSL2] over interval $I$ is given by: + +$$ +c(y_{I}) = \sum_{t\in I} \|G(y_t - \bar{y})\|_2^2 +$$ + +where $\bar{y}$ is the mean of $\{y_t\}_{t\in I}$. Note that $G: \mathbb{R}^p \rightarrow \mathbb{R}^p$ is linear. + +## Usage + +Start with the usual imports and create a graph signal. Note that the signal dimension must equal the number of nodes of the underlying graph. + +*Note*: the below graph and signal are meaningless. For relevant use-cases, see [Graph signal change point detection](../../examples/merging-cost-functions.ipynb). + +```python +import numpy as np +import networkx as nx +import matplotlib.pylab as plt +import ruptures as rpt + +# creation of the graph +nb_nodes = 30 +G = nx.gnp_random_graph(n=nb_nodes, p=0.5) + +# creation of the signal +n, dim = 500, nb_nodes # number of samples, dimension +n_bkps, sigma = 3, 1 # number of change points, noise standart deviation +signal, bkps = rpt.pw_constant(n, dim, n_bkps, noise_std=sigma) +``` + +Then create a [`CostGFSSL2`][ruptures.costs.costrbf.CostGFSSL2] instance and print the cost of the sub-signal `signal[50:150]`. The initialization of the class requires the Laplacian matrix of the underlying graph and the value of the cut-sparsity $\rho$ should be set based on the eigenvalues of $L$ (see [Graph signal change point detection](../../examples/merging-cost-functions.ipynb)). + +```python +# creation of the cost function instance +rho = 1 # the cut-sparsity +c = rpt.costs.CostGFSSL2(nx.laplacian_matrix(G).toarray(), rho) +c.fit(signal) +print(c.error(50, 150)) +``` + +You can also compute the sum of costs for a given list of change points. + +```python +print(c.sum_of_costs(bkps)) +print(c.sum_of_costs([10, 100, 200, 250, n])) +``` + +In order to use this cost class in a change point detection algorithm (inheriting from [`BaseEstimator`][ruptures.base.BaseEstimator]), pass a [`CostGFSSL2`][ruptures.costs.costrbf.CostGFSSL2] instance (through the argument `custom_cost`). + +```python +c = rpt.costs.CostL2() +algo = rpt.Dynp(custom_cost=c) +``` + +## References + +[Sharpnack2016] +Sharpnack, J., Rinaldo, A., and Singh, A. (2016). Detecting Anomalous Activity on Networks With the Graph Fourier Scan Statistic. EEE Transactions on Signal Processing, 64(2):364–379. + +[Shuman2013] +Shuman, D. I., Narang, S. K., Frossard, P., Ortega, A., and Vandergheynst, P. (2013). The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains. EEE Signal Processing Magazine, 30(3):83–98. From d47d7da5e5c37481a7580b3afddfac76b01c0941 Mon Sep 17 00:00:00 2001 From: evenmatencio Date: Thu, 30 May 2024 17:39:40 +0200 Subject: [PATCH 03/10] syntax and spell checks --- docs/user-guide/costs/costgfssl2.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/user-guide/costs/costgfssl2.md b/docs/user-guide/costs/costgfssl2.md index ed810b03..a5f3b33d 100644 --- a/docs/user-guide/costs/costgfssl2.md +++ b/docs/user-guide/costs/costgfssl2.md @@ -4,7 +4,7 @@ This cost function is specifically designed to detect localized mean-shifts in a graph signal. It relies on the least squared deviation of a graph signal filtered with a low-pass graph spectral filter $h$. -Formally, let $G = (V, E, W)$ a graph containing $p =|V|$ nodes, with a weighted adjacency matrix $W$. We define its combinatorial Laplacian matrix $L$ as classicaly done in Graph Signal Processing [[Shuman2013](#Shuman2013)]: +Formally, let $G = (V, E, W)$ a graph containing $p =|V|$ nodes, with a weighted adjacency matrix $W$. We define its combinatorial Laplacian matrix $L$ as classically done in Graph Signal Processing [[Shuman2013](#Shuman2013)]: $$ L = D - W = U \Lambda U^T, @@ -19,10 +19,10 @@ where Let $\{y_t\}_t, ~ y_t \in \mathbb{R}^p$, a multivariate signal on an interval $I$. The Graph Fourier Scan Statistic (GFSS) [[Sharpnack2016](#Sharpnack2016)] of a graph signal $y_t$ is $\|G(y_t)\|_2^2$ with: $$ -G(y) = \sum_{i=2}^p h(\lambda_i) (u_i^T y) u_i \quad \text{ and } \quad h(\lambda) = \min \left\{ 1, \sqrt{\frac{\rho}{\lambda}} \right\} +G(y) = \sum_{i=2}^p h(\lambda_i) (u_i^T y) u_i \quad \text{ and } \quad h(\lambda) = \min \left( 1, \sqrt{\frac{\rho}{\lambda}} \right) $$ -where $\rho$ is the so-called cut-sparsity. The cost function [`CostGFSSL2`][ruptures.costs.costrbf.CostGFSSL2] over interval $I$ is given by: +where $\rho$ is the so-called cut-sparsity. The cost function [`CostGFSSL2`][ruptures.costs.costgfssl2.CostGFSSL2] over interval $I$ is given by: $$ c(y_{I}) = \sum_{t\in I} \|G(y_t - \bar{y})\|_2^2 @@ -48,11 +48,11 @@ G = nx.gnp_random_graph(n=nb_nodes, p=0.5) # creation of the signal n, dim = 500, nb_nodes # number of samples, dimension -n_bkps, sigma = 3, 1 # number of change points, noise standart deviation +n_bkps, sigma = 3, 1 # number of change points, noise standard deviation signal, bkps = rpt.pw_constant(n, dim, n_bkps, noise_std=sigma) ``` -Then create a [`CostGFSSL2`][ruptures.costs.costrbf.CostGFSSL2] instance and print the cost of the sub-signal `signal[50:150]`. The initialization of the class requires the Laplacian matrix of the underlying graph and the value of the cut-sparsity $\rho$ should be set based on the eigenvalues of $L$ (see [Graph signal change point detection](../../examples/merging-cost-functions.ipynb)). +Then create a [`CostGFSSL2`][ruptures.costs.costgfssl2.CostGFSSL2] instance and print the cost of the sub-signal `signal[50:150]`. The initialization of the class requires the Laplacian matrix of the underlying graph and the value of the cut-sparsity $\rho$ should be set based on the eigenvalues of $L$ (see [Graph signal change point detection](../../examples/merging-cost-functions.ipynb)). ```python # creation of the cost function instance @@ -69,7 +69,7 @@ print(c.sum_of_costs(bkps)) print(c.sum_of_costs([10, 100, 200, 250, n])) ``` -In order to use this cost class in a change point detection algorithm (inheriting from [`BaseEstimator`][ruptures.base.BaseEstimator]), pass a [`CostGFSSL2`][ruptures.costs.costrbf.CostGFSSL2] instance (through the argument `custom_cost`). +In order to use this cost class in a change point detection algorithm (inheriting from [`BaseEstimator`][ruptures.base.BaseEstimator]), pass a [`CostGFSSL2`][ruptures.costs.costgfssl2.CostGFSSL2] instance (through the argument `custom_cost`). ```python c = rpt.costs.CostL2() From 9a3ac24b01f80ae151b89ecd0c06b947da8ba16c Mon Sep 17 00:00:00 2001 From: evenmatencio Date: Thu, 30 May 2024 17:40:13 +0200 Subject: [PATCH 04/10] first version of gfss example usage --- docs/examples/graph-signal-cpd.ipynb | 184 +++++++++++++++++++++++++++ 1 file changed, 184 insertions(+) create mode 100644 docs/examples/graph-signal-cpd.ipynb diff --git a/docs/examples/graph-signal-cpd.ipynb b/docs/examples/graph-signal-cpd.ipynb new file mode 100644 index 00000000..c0ccb892 --- /dev/null +++ b/docs/examples/graph-signal-cpd.ipynb @@ -0,0 +1,184 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Graph signals change point detection with the Graph Fourier Scan Statistic: a low-pass band filter\n", + "\n", + "## Introduction\n", + "\n", + "Graph signal processing (GSP) is the study of multivariate signals $y_t \\in \\mathbb{R}^d$ lying on the nodes of a graph $\\mathcal{G} = (\\mathcal{V}, \\mathcal{E}, \\mathcal{W})$ (see for instance [[Stankovic2019](#Stankovic2019), [Shuman2013](#Shuman2013)]). As in standard signal processing, it is possible to define a Graph Fourier Transform and to generalize the notion of spectral filtering. The intuition behind the graph spectral frequencies is that a signal whose (graph) spectrum is located at low frequencies is \"smoother\" than a signal whose energy is concentrated on high frequencies. By \"smoother\", we refer to the notion of smoothness with respect to the structure of the graph [[Shuman2013](#Shuman2013)]: the smoother a signal, the closer its values on neighbor nodes. \n", + "\n", + "Thus, by applying a low-pass filter on the graph spectral domain, one is likely to remove from a graph signal the noise and/or uncorrelated information across the nodes. The authors of [[Ferrari2019](#Ferrari2019)] leverage this idea to define the Graph Fourier Scan Statistic (GFSS) algorithm (derived from the statistic introduced in [[Sharpnack2016](#Sharpnack2016)]). When a graph structure is available, this is one possible way of using notions coming from the field of GSP to enhance change point detection for multivariate signals.\n", + "\n", + "In what follows, we focus on the above approach and we show how to apply it with `ruptures`. This example relies on the class [CostGFSSL2](../user-guide/costs/costgfssl2.md), which results from the combination of the GFSS and the least squared deviation.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Illustration and objectives\n", + "\n", + "First, we briefly illustrate the behavior of the GFSS and we justify the usage of the cost function `CostGFSSL2`. For a formal definition, please see [CostGFSSL2](../user-guide/costs/costgfssl2.md). \n", + "\n", + "The application of the GFSS amounts to a low-pass graph spectral filtering parametrized by the so-called cut-sparsity $\\rho$. The corresponding filter is displayed below:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "rho = 1\n", + "filter = lambda x, rho: np.minimum(1, np.sqrt(rho / x))\n", + "x = np.linspace(0, 10, 100)\n", + "filtered = [1] + list(filter(x[1:], rho))\n", + "\n", + "fig, ax = plt.subplots(1, 1, figsize=(6, 3))\n", + "ax.plot(x, filtered, label=\"GFSS filter\")\n", + "ax.axvline(x=rho, linestyle=\"--\", c=\"k\", label=\"$\\\\rho$\")\n", + "ax.set_xlabel(\"eigenvalues $\\lambda$\")\n", + "ax.legend()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As explained in the [introduction](#introduction), applying a low-pass filter to a graph signal shrinks the high-frequency components, thus attenuating the signal components that are not smooth with respect to the graph structure. Based on this statement, we deduce two potential (related) benefits of applying the cost `CostGFSSL2`:\n", + "\n", + "1. as in [[Ferrari2019](#Ferrari2019)], detecting mean changes that are localized on clusters of the graph only. Formally if we denote $m_t(i)$ the mean of the process at node $i$ and $\\mathcal{C}$ a well connected subset of $\\mathcal{V}$ (a cluster), one may try to detect $t_r$ such that:\n", + "\n", + "$$\n", + "y_t = m_t + e_t \\quad \\text{ with } ~ m_t = \n", + "\\begin{cases}\n", + " m & \\forall t < t_r \\\\ \n", + " m + \\delta & \\forall t \\geq t_r \n", + "\\end{cases}\n", + "\\quad \\text{ and } ~ \\delta_i = \n", + "\\begin{cases}\n", + " c > 0 & \\forall i \\in \\mathcal{C} \\\\ \n", + " 0 & \\text{ otherwise } \n", + "\\end{cases}\n", + "$$\n", + "\n", + "2. attenuating changes induced by spatially white noise (with high variance) or that may be due to individual dysfunctions of the observed system, for instance a geographical censors network, a social network..." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup\n", + "\n", + "We generate a synthetic graph matching the above description and we define a signal over it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import ruptures as rpt # our package\n", + "import networkx as nx # for graph utils" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Graph generation\n", + "\n", + "nb_nodes = 120\n", + "cluster_nb = 6\n", + "mean_cluster_size = 20\n", + "inter_density = 0.02 # density of inter-clusters edges\n", + "intra_density = 0.9 # density of intra-clusters edges\n", + "graph_seed = 9 # for reproducibility\n", + "G = nx.gaussian_random_partition_graph(\n", + " n=nb_nodes,\n", + " s=mean_cluster_size,\n", + " v=2 * mean_cluster_size,\n", + " p_in=intra_density,\n", + " p_out=inter_density,\n", + " seed=graph_seed,\n", + ")\n", + "coord = nx.spring_layout(G) # for plotting" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Vizualization of the graph clusters\n", + "\n", + "clusters_seed = 20\n", + "clusters = nx.algorithms.community.louvain.louvain_communities(G, seed=clusters_seed)\n", + "colors_dct = {0: \"r\", 1: \"b\", 2: \"g\", 3: \"orange\", 4: \"purple\", 5: \"brown\"}\n", + "cluster_idx_arr = np.zeros((nb_nodes))\n", + "\n", + "for cl_ind in range(len(clusters)):\n", + " for node_ind in list(clusters[cl_ind]):\n", + " cluster_idx_arr[node_ind] = cl_ind\n", + "\n", + "colors_l = [colors_dct[cluster_idx_arr[node_ind]] for node_ind in range(nb_nodes)]\n", + "\n", + "fig, ax = plt.subplots(1, 1, figsize=(8, 5))\n", + "ax.set_title(\"Clusters vizualization\")\n", + "nx.draw_networkx(G, pos=coord, with_labels=True, node_color=colors_l, ax=ax)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## References\n", + "\n", + "[Ferrari2019]\n", + "Ferrari, A., Richard, C., and Verduci, L. (2019). Distributed Change Detection in Streaming Graph Signals. IEEE 8th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), pages 166–170.\n", + "\n", + "[Sharpnack2016]\n", + "Sharpnack, J., Rinaldo, A., and Singh, A. (2016). Detecting Anomalous Activity on Networks With the Graph Fourier Scan Statistic. EEE Transactions on Signal Processing, 64(2):364–379.\n", + "\n", + "[Shuman2013]\n", + "Shuman, D. I., Narang, S. K., Frossard, P., Ortega, A., and Vandergheynst, P. (2013). The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains. EEE Signal Processing Magazine, 30(3):83–98.\n", + "\n", + "[Stankovic2019]\n", + "Ljubisa Stankovic, Danilo P. Mandic, Milos Dakovic, Ilia Kisil, Ervin Sejdic, and Anthony G. Constantinides (2019). Understanding the Basis of Graph Signal Processing via an Intuitive Example-Driven Approach [Lecture Notes]. IEEE Signal Processing Magazine, 36(6):133–145.\n", + "\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From f11f68eb0c2a990d6469627832ede7e8ce7029f1 Mon Sep 17 00:00:00 2001 From: evenmatencio Date: Fri, 31 May 2024 14:40:42 +0200 Subject: [PATCH 05/10] spell checks and naming changes --- docs/user-guide/costs/costgfssl2.md | 2 +- src/ruptures/costs/costgfssl2.py | 24 ++++++++++++------------ 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/docs/user-guide/costs/costgfssl2.md b/docs/user-guide/costs/costgfssl2.md index a5f3b33d..5d8f2201 100644 --- a/docs/user-guide/costs/costgfssl2.md +++ b/docs/user-guide/costs/costgfssl2.md @@ -12,7 +12,7 @@ $$ where -- $D$ is the diagional degree matrix of the graph: $D = \text{diag}(d_1, \ldots, d_p)$. +- $D$ is the diagonal degree matrix of the graph: $D = \text{diag}(d_1, \ldots, d_p)$. - $U$ is the orthogonal matrix whose columns $\{u_i\}_{i=1}^p$ are the eigenvectors of $L$ - $\Lambda = \text{diag}(\lambda_1, \ldots, \lambda_p)$ contains the eigenvalues of $L$. diff --git a/src/ruptures/costs/costgfssl2.py b/src/ruptures/costs/costgfssl2.py index bfc49cbd..031ef4fc 100644 --- a/src/ruptures/costs/costgfssl2.py +++ b/src/ruptures/costs/costgfssl2.py @@ -25,16 +25,16 @@ def __init__(self, laplacian_mat, cut_sparsity) -> None: self.cut_sparsity = cut_sparsity self.graph_laplacian_mat = laplacian_mat self.signal = None - self.gft_square_cumsum = None - self.gft_cumsum = None + self.gfss_square_cumsum = None + self.gfss_cumsum = None super().__init__() def filter(self, freqs, eps=0.00001): - """Applies the GFSS filter to the input (spatial) frequences. - NOTE: the frequences must be in increasing order. + """Applies the GFSS filter to the input (spatial) frequencies. + NOTE: the frequencies must be in increasing order. Args: - freqs (array): ordered frequences to filter. + freqs (array): ordered frequencies to filter. eps (float, optional): threshold for non zero values. Defaults to 0.00001. Returns: @@ -59,13 +59,13 @@ def fit(self, signal): # Computation of the GFSS eigvals, eigvects = eigh(self.graph_laplacian_mat) filter_matrix = np.diag(self.filter(eigvals), k=0) - gft = filter_matrix.dot(eigvects.T.dot(signal.T)).T + gfss = filter_matrix.dot(eigvects.T.dot(signal.T)).T # Computation of the per-segment cost utils - self.gft_square_cumsum = np.concatenate( - [np.zeros((1, signal.shape[1])), np.cumsum(gft**2, axis=0)], axis=0 + self.gfss_square_cumsum = np.concatenate( + [np.zeros((1, signal.shape[1])), np.cumsum(gfss**2, axis=0)], axis=0 ) - self.gft_cumsum = np.concatenate( - [np.zeros((1, signal.shape[1])), np.cumsum(gft, axis=0)], axis=0 + self.gfss_cumsum = np.concatenate( + [np.zeros((1, signal.shape[1])), np.cumsum(gfss, axis=0)], axis=0 ) return self @@ -83,6 +83,6 @@ def error(self, start, end): if end - start < self.min_size: raise NotEnoughPoints - sub_square_sum = self.gft_square_cumsum[end] - self.gft_square_cumsum[start] - sub_sum = self.gft_cumsum[end] - self.gft_cumsum[start] + sub_square_sum = self.gfss_square_cumsum[end] - self.gfss_square_cumsum[start] + sub_sum = self.gfss_cumsum[end] - self.gfss_cumsum[start] return np.sum(sub_square_sum - (sub_sum**2) / (end - start)) From eec9efaec094235988fd0c1d09241209b1fe44d8 Mon Sep 17 00:00:00 2001 From: evenmatencio Date: Fri, 31 May 2024 14:45:10 +0200 Subject: [PATCH 06/10] code cells added, to be documented --- docs/examples/graph-signal-cpd.ipynb | 215 ++++++++++++++++++++++++++- 1 file changed, 207 insertions(+), 8 deletions(-) diff --git a/docs/examples/graph-signal-cpd.ipynb b/docs/examples/graph-signal-cpd.ipynb index c0ccb892..ab050704 100644 --- a/docs/examples/graph-signal-cpd.ipynb +++ b/docs/examples/graph-signal-cpd.ipynb @@ -19,7 +19,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Illustration and objectives\n", + "### Illustration\n", "\n", "First, we briefly illustrate the behavior of the GFSS and we justify the usage of the cost function `CostGFSSL2`. For a formal definition, please see [CostGFSSL2](../user-guide/costs/costgfssl2.md). \n", "\n", @@ -51,9 +51,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "As explained in the [introduction](#introduction), applying a low-pass filter to a graph signal shrinks the high-frequency components, thus attenuating the signal components that are not smooth with respect to the graph structure. Based on this statement, we deduce two potential (related) benefits of applying the cost `CostGFSSL2`:\n", + "### Objectives\n", "\n", - "1. as in [[Ferrari2019](#Ferrari2019)], detecting mean changes that are localized on clusters of the graph only. Formally if we denote $m_t(i)$ the mean of the process at node $i$ and $\\mathcal{C}$ a well connected subset of $\\mathcal{V}$ (a cluster), one may try to detect $t_r$ such that:\n", + "As explained above, applying a low-pass filter to a graph signal shrinks the high-frequency components, thus attenuating the signal components that are not smooth with respect to the graph structure. Based on this statement, we deduce two potential (related) benefits of applying the cost `CostGFSSL2`:\n", + "\n", + "1. *Scenario 1*: as in [[Ferrari2019](#Ferrari2019)], detecting mean changes that are localized on clusters of the graph only. Formally if we denote $m_t(i)$ the mean of the process at node $i$ and $\\mathcal{C}$ a well connected subset of $\\mathcal{V}$ (a cluster), one may try to detect $t_r$ such that:\n", "\n", "$$\n", "y_t = m_t + e_t \\quad \\text{ with } ~ m_t = \n", @@ -68,7 +70,7 @@ "\\end{cases}\n", "$$\n", "\n", - "2. attenuating changes induced by spatially white noise (with high variance) or that may be due to individual dysfunctions of the observed system, for instance a geographical censors network, a social network..." + "2. *Scenario 2*: attenuating changes induced by spatially white noise (with high variance) or that may be due to individual dysfunctions of the observed system, for instance a geographical censors network, a social network... in graphs showing clusters." ] }, { @@ -77,7 +79,7 @@ "source": [ "## Setup\n", "\n", - "We generate a synthetic graph matching the above description and we define a signal over it." + "We generate a synthetic graph matching the above description, namely with a highly clusterized structure." ] }, { @@ -87,7 +89,9 @@ "outputs": [], "source": [ "import ruptures as rpt # our package\n", - "import networkx as nx # for graph utils" + "import networkx as nx # for graph utils\n", + "\n", + "from ruptures.costs.costgfssl2 import CostGFSSL2 # the gfss based cost function" ] }, { @@ -102,7 +106,7 @@ "cluster_nb = 6\n", "mean_cluster_size = 20\n", "inter_density = 0.02 # density of inter-clusters edges\n", - "intra_density = 0.9 # density of intra-clusters edges\n", + "intra_density = 0.95 # density of intra-clusters edges\n", "graph_seed = 9 # for reproducibility\n", "G = nx.gaussian_random_partition_graph(\n", " n=nb_nodes,\n", @@ -112,7 +116,7 @@ " p_out=inter_density,\n", " seed=graph_seed,\n", ")\n", - "coord = nx.spring_layout(G) # for plotting" + "coord = nx.spring_layout(G, seed=graph_seed) # for plotting consistency" ] }, { @@ -139,6 +143,196 @@ "nx.draw_networkx(G, pos=coord, with_labels=True, node_color=colors_l, ax=ax)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Robustness with respect to noise in Scenario 1\n", + "\n", + "We aim at detecting mean changes localized in a cluster only, as described in the scenario 1 presented in [the objectives](#objectives)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "here explain PELT & penalty coefficient & rho" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "- to explain that we use PELT, assuming we do not know the number of ruptures\n", + "- say that we do not bother to correctly compute the penalty coefficient\n", + "- say that we look for a good cut-sparsityz" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.linalg import eigh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "eigvals, eigvects = eigh(nx.laplacian_matrix(G).toarray())\n", + "print(eigvals)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "## SIGNAL GENERATION\n", + "\n", + "n_dims = nb_nodes\n", + "dims_with_change = np.array([dim for dim in clusters[-1]])\n", + "n_dims_with_change = len(dims_with_change)\n", + "n_samples = 100\n", + "signal_seed = 10\n", + "\n", + "print(f\"Number of nodes: {nb_nodes}, \")\n", + "print(f\"The dimensions with changes are: {[dim for dim in dims_with_change]} \")\n", + "\n", + "bic_L2_pen = n_dims / 2 * np.log(n_samples)\n", + "pen = 2 * bic_L2_pen\n", + "\n", + "noise_std_values = [0.1, 1, 2, 3, 4]\n", + "\n", + "rho = 1\n", + "print(f\"We set the cut-sparsity rho = {rho}. \\n\")\n", + "\n", + "for noise_std in noise_std_values:\n", + "\n", + " signal_with_change, gt_bkps = rpt.pw_constant(\n", + " n_samples, n_dims_with_change, 1, noise_std=noise_std, seed=signal_seed\n", + " )\n", + " signal, _ = rpt.pw_constant(\n", + " n_samples, n_dims, 0, noise_std=noise_std, seed=signal_seed\n", + " )\n", + " signal[:, dims_with_change] = signal_with_change\n", + "\n", + " ## CHANGE POINT DETECTION\n", + "\n", + " # with graph and GFSS\n", + " laplacian_matrix = nx.laplacian_matrix(\n", + " G\n", + " ).toarray() # extract the laplacian matrix as an numpy array\n", + " cost_rpt_pelt = CostGFSSL2(laplacian_matrix, rho)\n", + " graph_algo = rpt.Pelt(custom_cost=cost_rpt_pelt, jump=1, min_size=1).fit(signal)\n", + " my_graph_bkps = graph_algo.predict(pen=pen)\n", + "\n", + " # without graph\n", + " algo = rpt.Pelt(model=\"l2\", jump=1, min_size=1).fit(signal)\n", + " my_bkps = algo.predict(pen=pen)\n", + "\n", + " print(\n", + " \"NOISE_STD=\",\n", + " noise_std,\n", + " \"\\tGROUNDTRUTH\",\n", + " gt_bkps,\n", + " \"\\tWITH GRAPH: \",\n", + " my_graph_bkps,\n", + " \"\\tWITHOUT GRAPH: \",\n", + " my_bkps,\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Influence of $\\rho$ in Scenario 2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "## SIGNAL GENERATION AND HYPER-PARAMETERS TUNING\n", + "\n", + "rng = np.random.default_rng(seed=10)\n", + "\n", + "# signal hyper-parameters\n", + "n_dims = nb_nodes # number of dimensions equals the number of nodes\n", + "n_dims_with_change = n_dims // 20 # number of nodes undergoing a mean change\n", + "dims_with_change = np.arange(n_dims)\n", + "rng.shuffle(dims_with_change) # randomizing the nodes with change\n", + "dims_with_change = dims_with_change[:n_dims_with_change]\n", + "\n", + "# building the signal itself\n", + "n_samples = 100\n", + "signal_seed = 10\n", + "noise_std = 1\n", + "\n", + "signal_with_change, gt_bkps = rpt.pw_constant(\n", + " n_samples, n_dims_with_change, 1, noise_std=noise_std, seed=signal_seed\n", + ")\n", + "signal, _ = rpt.pw_constant(n_samples, n_dims, 0, noise_std=noise_std, seed=signal_seed)\n", + "signal[:, dims_with_change] = signal_with_change\n", + "\n", + "print(f\"The dimensions with changes are: {[dim for dim in dims_with_change]}\")\n", + "\n", + "# algorithm hyper-parameters\n", + "bic_L2_pen = n_dims / 2 * np.log(n_samples)\n", + "pen = 2 * bic_L2_pen\n", + "\n", + "rho_values = [\n", + " eigvals[1] / 10,\n", + " eigvals[1] / 2,\n", + " eigvals[1],\n", + " 2 * eigvals[1],\n", + " eigvals[60],\n", + " eigvals[-1],\n", + "]\n", + "\n", + "for rho in rho_values:\n", + "\n", + " ## CHANGE POINT DETECTION\n", + "\n", + " # with graph and GFSS\n", + " laplacian_matrix = nx.laplacian_matrix(\n", + " G\n", + " ).toarray() # extract the laplacian matrix as an numpy array\n", + " cost_rpt_pelt = CostGFSSL2(laplacian_matrix, rho)\n", + " graph_algo = rpt.Pelt(custom_cost=cost_rpt_pelt, jump=1, min_size=1).fit(signal)\n", + " my_graph_bkps = graph_algo.predict(pen=pen)\n", + "\n", + " # without graph\n", + " algo = rpt.Pelt(model=\"l2\", jump=1, min_size=1).fit(signal)\n", + " my_bkps = algo.predict(pen=pen)\n", + "\n", + " print(\n", + " \"RHO=\",\n", + " round(rho, ndigits=3),\n", + " \"\\tGROUNDTRUTH\",\n", + " gt_bkps,\n", + " \"\\tWITH GRAPH: \",\n", + " my_graph_bkps,\n", + " \"\\tWITHOUT GRAPH: \",\n", + " my_bkps,\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Conclusions" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -158,6 +352,11 @@ "Ljubisa Stankovic, Danilo P. Mandic, Milos Dakovic, Ilia Kisil, Ervin Sejdic, and Anthony G. Constantinides (2019). Understanding the Basis of Graph Signal Processing via an Intuitive Example-Driven Approach [Lecture Notes]. IEEE Signal Processing Magazine, 36(6):133–145.\n", "\n" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] } ], "metadata": { From 7246a4c5f975e58ba94bd3e0c124efd5e746cdb2 Mon Sep 17 00:00:00 2001 From: Nicolas Cecchi Date: Fri, 31 May 2024 16:20:16 +0200 Subject: [PATCH 07/10] Change filter formula --- src/ruptures/costs/costgfssl2.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/ruptures/costs/costgfssl2.py b/src/ruptures/costs/costgfssl2.py index 031ef4fc..676db61a 100644 --- a/src/ruptures/costs/costgfssl2.py +++ b/src/ruptures/costs/costgfssl2.py @@ -42,7 +42,7 @@ def filter(self, freqs, eps=0.00001): """ nb_zeros = np.sum(freqs < eps) filtered_freqs = np.minimum(1, np.sqrt(self.cut_sparsity / freqs[nb_zeros:])) - return np.concatenate([np.zeros(nb_zeros), filtered_freqs]) + return np.concatenate([np.ones(nb_zeros), filtered_freqs]) def fit(self, signal): """Performs pre-computations for per-segment approximation cost. @@ -86,3 +86,9 @@ def error(self, start, end): sub_square_sum = self.gfss_square_cumsum[end] - self.gfss_square_cumsum[start] sub_sum = self.gfss_cumsum[end] - self.gfss_cumsum[start] return np.sum(sub_square_sum - (sub_sum**2) / (end - start)) + + + +# %% + +# %% From 78163cd16d478b14c538c46d9b14aff496cb1dfc Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 31 May 2024 14:39:48 +0000 Subject: [PATCH 08/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/ruptures/costs/costgfssl2.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ruptures/costs/costgfssl2.py b/src/ruptures/costs/costgfssl2.py index 676db61a..c9b6adc9 100644 --- a/src/ruptures/costs/costgfssl2.py +++ b/src/ruptures/costs/costgfssl2.py @@ -88,7 +88,6 @@ def error(self, start, end): return np.sum(sub_square_sum - (sub_sum**2) / (end - start)) - # %% # %% From 29becacd11d1e4bf44edb73e3cba50039b398fec Mon Sep 17 00:00:00 2001 From: evenmatencio Date: Mon, 3 Jun 2024 11:42:11 +0200 Subject: [PATCH 09/10] complete explanations, to be adapted with real data --- docs/examples/graph-signal-cpd.ipynb | 113 +++++++++++++-------------- 1 file changed, 53 insertions(+), 60 deletions(-) diff --git a/docs/examples/graph-signal-cpd.ipynb b/docs/examples/graph-signal-cpd.ipynb index ab050704..4151706f 100644 --- a/docs/examples/graph-signal-cpd.ipynb +++ b/docs/examples/graph-signal-cpd.ipynb @@ -10,7 +10,7 @@ "\n", "Graph signal processing (GSP) is the study of multivariate signals $y_t \\in \\mathbb{R}^d$ lying on the nodes of a graph $\\mathcal{G} = (\\mathcal{V}, \\mathcal{E}, \\mathcal{W})$ (see for instance [[Stankovic2019](#Stankovic2019), [Shuman2013](#Shuman2013)]). As in standard signal processing, it is possible to define a Graph Fourier Transform and to generalize the notion of spectral filtering. The intuition behind the graph spectral frequencies is that a signal whose (graph) spectrum is located at low frequencies is \"smoother\" than a signal whose energy is concentrated on high frequencies. By \"smoother\", we refer to the notion of smoothness with respect to the structure of the graph [[Shuman2013](#Shuman2013)]: the smoother a signal, the closer its values on neighbor nodes. \n", "\n", - "Thus, by applying a low-pass filter on the graph spectral domain, one is likely to remove from a graph signal the noise and/or uncorrelated information across the nodes. The authors of [[Ferrari2019](#Ferrari2019)] leverage this idea to define the Graph Fourier Scan Statistic (GFSS) algorithm (derived from the statistic introduced in [[Sharpnack2016](#Sharpnack2016)]). When a graph structure is available, this is one possible way of using notions coming from the field of GSP to enhance change point detection for multivariate signals.\n", + "Thus, by applying a low-pass filter on the graph spectral domain, one is likely to remove from a graph signal the noise and/or uncorrelated information across the nodes. The authors of [[Ferrari2019](#Ferrari2019)] leverage this idea to define the Graph Fourier Scan Statistic (GFSS) algorithm (derived from the statistic introduced in [[Sharpnack2016](#Sharpnack2016)]). When a graph structure is available, this is one possible way of using notions coming from the field of GSP to enhance change point detection (CPD) for multivariate signals.\n", "\n", "In what follows, we focus on the above approach and we show how to apply it with `ruptures`. This example relies on the class [CostGFSSL2](../user-guide/costs/costgfssl2.md), which results from the combination of the GFSS and the least squared deviation.\n" ] @@ -55,7 +55,7 @@ "\n", "As explained above, applying a low-pass filter to a graph signal shrinks the high-frequency components, thus attenuating the signal components that are not smooth with respect to the graph structure. Based on this statement, we deduce two potential (related) benefits of applying the cost `CostGFSSL2`:\n", "\n", - "1. *Scenario 1*: as in [[Ferrari2019](#Ferrari2019)], detecting mean changes that are localized on clusters of the graph only. Formally if we denote $m_t(i)$ the mean of the process at node $i$ and $\\mathcal{C}$ a well connected subset of $\\mathcal{V}$ (a cluster), one may try to detect $t_r$ such that:\n", + "1. *Scenario 1*: as in [[Ferrari2019](#Ferrari2019)], detecting mean changes that are localized on clusters of the graph only. Formally, let $e_t \\sim \\mathcal{N}(0, \\sigma^2 I)$ a white noise process. If we denote $m_t(i)$ the mean of the process at node $i$ and $\\mathcal{C}$ a well connected subset of $\\mathcal{V}$ (a cluster), one may try to detect $t_r$ such that:\n", "\n", "$$\n", "y_t = m_t + e_t \\quad \\text{ with } ~ m_t = \n", @@ -70,7 +70,7 @@ "\\end{cases}\n", "$$\n", "\n", - "2. *Scenario 2*: attenuating changes induced by spatially white noise (with high variance) or that may be due to individual dysfunctions of the observed system, for instance a geographical censors network, a social network... in graphs showing clusters." + "2. *Scenario 2*: attenuating changes induced by spatially white noise (with high variance) or changes that may be due to individual dysfunctions of the observed system, for instance a geographical censors network, a social network... in graphs showing clusters." ] }, { @@ -79,7 +79,7 @@ "source": [ "## Setup\n", "\n", - "We generate a synthetic graph matching the above description, namely with a highly clusterized structure." + "We generate a synthetic graph matching the above description, namely with a highly clusterized structure. The following graph contains $120$ nodes split between $6$ clusters." ] }, { @@ -149,32 +149,11 @@ "source": [ "#### Robustness with respect to noise in Scenario 1\n", "\n", - "We aim at detecting mean changes localized in a cluster only, as described in the scenario 1 presented in [the objectives](#objectives)." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "here explain PELT & penalty coefficient & rho" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "- to explain that we use PELT, assuming we do not know the number of ruptures\n", - "- say that we do not bother to correctly compute the penalty coefficient\n", - "- say that we look for a good cut-sparsityz" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from scipy.linalg import eigh" + "We aim at detecting mean changes localized over one cluster only, as described in the scenario 1 presented in [the objectives](#objectives). Therefore, the generated signal undergoes a single mean-change over only one of the graph clusters, while it remains unchanged (in mean) over the other nodes. In this experiment, we progressively increase the noise level of the generated signal and we compare the change points that are detected by the standard [`CostL2`](../../user-guide/costs/costl2) and those obtained with the [`CostGFSSL`](../user-guide/costs/costgfssl2.md). We expect the GFSS filtering to increase the robustness of the change-point detection.\n", + "\n", + "We assume to be in the most general case, namely as if we did not know the number of changes to be detected. In this framework, we use the [Pelt](../user-guide/detection/pelt.md) algorithm and we set the penalty coefficient $\\beta$ to an arbitrary constant such that the experiments we show are relevant. The rigorous calibration of this parameter would require more work but is not the topic of this example. \n", + "\n", + "The last parameter to be chosen is the cut-sparsity $\\rho$ used in the definition of the GFSS [[Sharpnack2016](#Sharpnack2016), [Ferrari2019](#Ferrari2019)]. Empirically, as this parameter acts as a frequency threshold, one should look at the eigenvalues (spatial frequencies) to get an idea for the magnitude of $\\rho$. For the relevancy of our two experiments we set: $\\beta = 500$ and $\\rho = \\lambda_1 / 2$ where $\\lambda_1$ is the first eigenvalue." ] }, { @@ -183,8 +162,11 @@ "metadata": {}, "outputs": [], "source": [ + "# Some eignevalues\n", + "from scipy.linalg import eigh\n", + "\n", "eigvals, eigvects = eigh(nx.laplacian_matrix(G).toarray())\n", - "print(eigvals)" + "print([eigvals[i] for i in [0, 1, 5, 10, 20, 50, 100, -1]])" ] }, { @@ -193,25 +175,24 @@ "metadata": {}, "outputs": [], "source": [ - "## SIGNAL GENERATION\n", + "## SIGNAL GENERATION AND HYPER-PARAMETERS SETUP\n", "\n", "n_dims = nb_nodes\n", "dims_with_change = np.array([dim for dim in clusters[-1]])\n", "n_dims_with_change = len(dims_with_change)\n", "n_samples = 100\n", "signal_seed = 10\n", + "noise_std_values = [0.1, 1, 2, 3, 4]\n", + "\n", + "# Algorithm hyper-parameters\n", + "pen = 500\n", + "rho = eigvals[1] / 2\n", "\n", "print(f\"Number of nodes: {nb_nodes}, \")\n", "print(f\"The dimensions with changes are: {[dim for dim in dims_with_change]} \")\n", - "\n", - "bic_L2_pen = n_dims / 2 * np.log(n_samples)\n", - "pen = 2 * bic_L2_pen\n", - "\n", - "noise_std_values = [0.1, 1, 2, 3, 4]\n", - "\n", - "rho = 1\n", "print(f\"We set the cut-sparsity rho = {rho}. \\n\")\n", "\n", + "print(\"RESULTS: \\n-------- \\n\")\n", "for noise_std in noise_std_values:\n", "\n", " signal_with_change, gt_bkps = rpt.pw_constant(\n", @@ -225,10 +206,8 @@ " ## CHANGE POINT DETECTION\n", "\n", " # with graph and GFSS\n", - " laplacian_matrix = nx.laplacian_matrix(\n", - " G\n", - " ).toarray() # extract the laplacian matrix as an numpy array\n", - " cost_rpt_pelt = CostGFSSL2(laplacian_matrix, rho)\n", + " laplacian_mat = nx.laplacian_matrix(G).toarray()\n", + " cost_rpt_pelt = CostGFSSL2(laplacian_mat, rho)\n", " graph_algo = rpt.Pelt(custom_cost=cost_rpt_pelt, jump=1, min_size=1).fit(signal)\n", " my_graph_bkps = graph_algo.predict(pen=pen)\n", "\n", @@ -241,9 +220,9 @@ " noise_std,\n", " \"\\tGROUNDTRUTH\",\n", " gt_bkps,\n", - " \"\\tWITH GRAPH: \",\n", + " \"\\tWITH GFSS: \",\n", " my_graph_bkps,\n", - " \"\\tWITHOUT GRAPH: \",\n", + " \"\\tSTANDARD L2: \",\n", " my_bkps,\n", " )" ] @@ -252,7 +231,18 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### Influence of $\\rho$ in Scenario 2" + "With the current value of the penalty coefficient and the chosen value of $\\rho$, we observe that the [CostGFSSL](../user-guide/costs/costgfssl2.md) cost function only detects the targeted change points while the standard [CostL2](../../user-guide/costs/costl2) cost is very sensitive to high noise levels. This result is necessarily dependent on the values of the hyper-parameters (other parametrization would lead to good results with the standard L2 cost function). In the [conclusion](#conclusions), we eventually justify why the above result should be seen as satisfying." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Influence of $\\rho$ in Scenario 2\n", + "\n", + "In the scenario 2, we assume that a very small number of nodes undergo a mean-change during the observation of the signal. The location of such nodes is random among the different clusters of the graph and may account for the occasional breakdown of some censors in a geographical network or an industrial system. In this framework, one may want to detect a precise physical phenomena and to avoid detecting random events like censor breakdown. \n", + "\n", + "The purpose of the following experiment is to compare the outputs of the CPD algorithm when applied with [`CostL2`](../../user-guide/costs/costl2) and [`CostGFSSL`](../user-guide/costs/costgfssl2.md) with respect to the value of $\\rho$. We use the same value of $\\beta$ than above and we explore different values of $\\rho$." ] }, { @@ -261,13 +251,13 @@ "metadata": {}, "outputs": [], "source": [ - "## SIGNAL GENERATION AND HYPER-PARAMETERS TUNING\n", + "## SIGNAL GENERATION AND HYPER-PARAMETERS SETUP\n", "\n", "rng = np.random.default_rng(seed=10)\n", "\n", "# signal hyper-parameters\n", "n_dims = nb_nodes # number of dimensions equals the number of nodes\n", - "n_dims_with_change = n_dims // 20 # number of nodes undergoing a mean change\n", + "n_dims_with_change = n_dims // 25 # number of nodes undergoing a mean change\n", "dims_with_change = np.arange(n_dims)\n", "rng.shuffle(dims_with_change) # randomizing the nodes with change\n", "dims_with_change = dims_with_change[:n_dims_with_change]\n", @@ -286,8 +276,7 @@ "print(f\"The dimensions with changes are: {[dim for dim in dims_with_change]}\")\n", "\n", "# algorithm hyper-parameters\n", - "bic_L2_pen = n_dims / 2 * np.log(n_samples)\n", - "pen = 2 * bic_L2_pen\n", + "pen = 500\n", "\n", "rho_values = [\n", " eigvals[1] / 10,\n", @@ -303,10 +292,10 @@ " ## CHANGE POINT DETECTION\n", "\n", " # with graph and GFSS\n", - " laplacian_matrix = nx.laplacian_matrix(\n", + " laplacian_mat = nx.laplacian_matrix(\n", " G\n", " ).toarray() # extract the laplacian matrix as an numpy array\n", - " cost_rpt_pelt = CostGFSSL2(laplacian_matrix, rho)\n", + " cost_rpt_pelt = CostGFSSL2(laplacian_mat, rho)\n", " graph_algo = rpt.Pelt(custom_cost=cost_rpt_pelt, jump=1, min_size=1).fit(signal)\n", " my_graph_bkps = graph_algo.predict(pen=pen)\n", "\n", @@ -319,9 +308,9 @@ " round(rho, ndigits=3),\n", " \"\\tGROUNDTRUTH\",\n", " gt_bkps,\n", - " \"\\tWITH GRAPH: \",\n", + " \"\\tWITH GFSS: \",\n", " my_graph_bkps,\n", - " \"\\tWITHOUT GRAPH: \",\n", + " \"\\tSTANDARD L2: \",\n", " my_bkps,\n", " )" ] @@ -330,7 +319,16 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Conclusions" + "When increasing the value of $\\rho$, the [CostGFSSL](../user-guide/costs/costgfssl2.md) cost function detects the \"unwanted\" mean-changes, but for low $\\rho$ the results are as expected. Similarly to above, this result is dependent on the parametrization, the proportion of nodes undergoing a mean change and the connectivity of the graph. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Conclusions\n", + "\n", + "- to say that we found a common set of hyper-parameters values for which both experiments give the expected results, even though this set of parameters depends on the graph itsefl. " ] }, { @@ -352,11 +350,6 @@ "Ljubisa Stankovic, Danilo P. Mandic, Milos Dakovic, Ilia Kisil, Ervin Sejdic, and Anthony G. Constantinides (2019). Understanding the Basis of Graph Signal Processing via an Intuitive Example-Driven Approach [Lecture Notes]. IEEE Signal Processing Magazine, 36(6):133–145.\n", "\n" ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [] } ], "metadata": { From f7236d6b03492a8de028b45512ab2d9206f23a8c Mon Sep 17 00:00:00 2001 From: evenmatencio Date: Mon, 3 Jun 2024 13:24:41 +0200 Subject: [PATCH 10/10] set value to rho if not given + sum from 1 --- src/ruptures/costs/costgfssl2.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/ruptures/costs/costgfssl2.py b/src/ruptures/costs/costgfssl2.py index c9b6adc9..3fcafa30 100644 --- a/src/ruptures/costs/costgfssl2.py +++ b/src/ruptures/costs/costgfssl2.py @@ -15,12 +15,13 @@ class CostGFSSL2(BaseCost): model = "gfss_l2_cost" min_size = 1 - def __init__(self, laplacian_mat, cut_sparsity) -> None: - """ + def __init__(self, laplacian_mat, cut_sparsity=None) -> None: + """_ Args: laplacian_mat (array): the discrete Laplacian matrix of the graph: D - W where D is the diagonal matrix diag(d_i) of the node degrees and W the adjacency matrix - cut_sparsity (float): frequency threshold of the GFSS spectral filter + cut_sparsity (float): frequency threshold of the GFSS spectral filter. Defaults to None. + If not specified, will be set to the first non null eigenvalue. """ self.cut_sparsity = cut_sparsity self.graph_laplacian_mat = laplacian_mat @@ -29,7 +30,7 @@ def __init__(self, laplacian_mat, cut_sparsity) -> None: self.gfss_cumsum = None super().__init__() - def filter(self, freqs, eps=0.00001): + def filter(self, freqs, eps=0.000001): """Applies the GFSS filter to the input (spatial) frequencies. NOTE: the frequencies must be in increasing order. @@ -41,8 +42,12 @@ def filter(self, freqs, eps=0.00001): filtered_freqs (array): the output of the filter. """ nb_zeros = np.sum(freqs < eps) + # Set the cut-sparsity if None + if self.cut_sparsity is None: + self.cut_sparsity = freqs[nb_zeros] + # Apply filtering filtered_freqs = np.minimum(1, np.sqrt(self.cut_sparsity / freqs[nb_zeros:])) - return np.concatenate([np.ones(nb_zeros), filtered_freqs]) + return np.concatenate([np.zeros(nb_zeros), filtered_freqs]) def fit(self, signal): """Performs pre-computations for per-segment approximation cost. @@ -86,8 +91,3 @@ def error(self, start, end): sub_square_sum = self.gfss_square_cumsum[end] - self.gfss_square_cumsum[start] sub_sum = self.gfss_cumsum[end] - self.gfss_cumsum[start] return np.sum(sub_square_sum - (sub_sum**2) / (end - start)) - - -# %% - -# %%