diff --git a/README.md b/README.md index 9e86c30..bd569f1 100644 --- a/README.md +++ b/README.md @@ -43,8 +43,7 @@ Current algorithms come from the field of: The package is aimed at execution performance (using JIT compilation and standing on the shoulders of giants like numpy -and scipy) while also keeping code readable and maintainable. This includes comments as well as architectural choices. -This might not be perfect, but we are trying! +and scipy) while also keeping code readable and maintainable. This includes comments as well as architectural choices. @@ -156,7 +155,9 @@ are currently under development. An overview with sources can be seen here: | [Fast SST](https://github.com/Lucew/changepoynt/blob/master/docs/fast-sst.md) | [Weber et al.](https://doi.org/10.1109/ACCESS.2025.3640386) | Stable :heavy_check_mark: | | IKA-SST | [Idé](https://epubs.siam.org/doi/abs/10.1137/1.9781611972771.54) | Stable :heavy_check_mark: | | [Fast IKA-SST](https://github.com/Lucew/changepoynt/blob/master/docs/fast-sst.md) | [Weber et al.](https://doi.org/10.1109/ACCESS.2025.3640386) | Stable :heavy_check_mark: | +| Multivariate (IKA-)SST (MSST) | | Stable :heavy_check_mark: | | ESST | [Boelter & Weber et al.](https://ntrs.nasa.gov/citations/20250002705) | Stable :heavy_check_mark: | +| Multivariate ESST (MESST) | | Stable :heavy_check_mark: | | RuLSIF | [Liu et al.](https://www.sciencedirect.com/science/article/pii/S0893608013000270) | Stable :heavy_check_mark: | | uLSIF | [Liu et al.](https://www.sciencedirect.com/science/article/pii/S0893608013000270) | Stable :heavy_check_mark: | | KLIEP | [Liu et al.](https://www.sciencedirect.com/science/article/pii/S0893608013000270) | Planned | diff --git a/changepoynt/algorithms/messt.py b/changepoynt/algorithms/messt.py new file mode 100644 index 0000000..9571e05 --- /dev/null +++ b/changepoynt/algorithms/messt.py @@ -0,0 +1,204 @@ +import os +from typing import Callable +from functools import partial + +import numba as nb +import numpy as np + +from changepoynt.algorithms.base_algorithm import Algorithm +from changepoynt.utils import normalization +import changepoynt.utils.block_linalg as blg +import changepoynt.algorithms.esst as esst + + +class MESST(Algorithm): + """ + This class implements the Multivariate Enhanced Singular Spectrum Transformation (MESST) + + The univariate algorithm ESST has been published together with Sarah Boelter (UMN) + Boelter, Sarah*, Lucas Weber*, et al. "Fault Prediction in Planetary Drilling Using Subspace Analysis Techniques." + International Conference on Intelligent Autonomous Systems (IAS-19). 2025. + (* is equal contribution) + https://ntrs.nasa.gov/citations/20250002705 + + It is a research algorithm, please use with caution. + """ + + def __init__(self, window_length: int, n_windows: int = None, lag: int = None, rank: int = 5, + scale: bool = True, method: str = 'rsvd', random_rank: int = None, scoring_step: int = 1, + use_fast_hankel: bool = False, threads: int = None) -> None: + """ + Experimental change point detection method evaluation the prevalence of change points within a signal + by comparing the difference in eigenvectors between to points in time. + + :param window_length: This specifies the length of the time series (in samples), which will be used to extract + the representative "eigensignals" to be compared before and after the lag. The window length should be big + enough to cover any wanted patterns (e.g., bigger than the periodicity of periodic signals). + + :param n_windows: This specifies the number of consecutive time windows used to extract the "eigensignals" from + the given time series. It should be big enough to cover different parts of the target behavior. If one does not + specify this value, we use the rule of thumb and take as many time windows as you specified the length of the + window (rule of thumb). + + :param lag: This parameter specifies the distance of the comparison for behaviors. In easy terms it tells the + algorithms how far it should look into the future to find change in behavior to the current signal. If you do + not specify this parameter, we use a rule of thumb and look ahead half of the window length you specified to + cover the behavior. + + :param rank: This parameter specifies the amount of "eigensignals" which will be used to measure the + dissimilarity of the signal in the future behavior. As a rule of thumb, we take the five most dominant + "eigensignals" if you do not specify otherwise. + + :param scale: Due to numeric stability, we REALLY RECOMMEND scaling the signal into a restricted value range. Per + default, we use a min max scaling to ensure a restricted value range. In the presence of extreme outliers, this + could cause problems, as the signal will be squished. + + :param random_rank: To use the randomized singular value decomposition, one needs to provide a + randomized rank (size of second matrix dimension for the randomized matrix) as specified in [3]. The lower + this value, the faster the computation but the higher the error (as the approximation gets worse). + + :param scoring_step: The distance between scoring steps in samples (e.g., 2 would half the computation). + + :param use_fast_hankel: Whether to deploy the fast hankel matrix product. + + :param threads: The number of threads the fast hankel matrix product is allowed to use. Default is the half of + the number of cpu cores your system has available. + """ + + # save the specified parameters into instance variables + self.window_length = window_length + self.n_windows = n_windows + self.rank = rank + self.scale = scale + self.random_rank = random_rank + self.lag = lag + self.scoring_step = scoring_step + self.use_fast_hankel = use_fast_hankel + self.method = method + self.threads = threads + + # set some default values when they have not been specified + if self.n_windows is None: + self.n_windows = self.window_length//2 + if self.lag is None: + # rule of thumb + self.lag = self.n_windows + if self.random_rank is None: + # compute the rank as specified in [3] and + # https://scikit-learn.org/stable/modules/generated/sklearn.utils.extmath.randomized_svd.html + self.random_rank = min(self.rank + 10, self.window_length, self.n_windows) + if self.threads is None: + self.threads = os.cpu_count()//2 + + # specify the methods and their corresponding functions as lambda functions expecting only the hankel matrix + self.methods = {'rsvd': partial(esst.left_entropy, rank=self.rank, random_rank=self.random_rank, + threads=self.threads, method=self.method)} + if self.method not in self.methods: + raise ValueError(f'Method {self.method} not defined. Possible methods: {list(self.methods.keys())}.') + + # set up the methods we use for the construction of the hankel matrix (either it is the fft representation + # of the other one) + if use_fast_hankel and self.method != 'rsvd': + raise ValueError(f'method {self.method} is not defined with use_fast_hankel=True') + + def transform(self, time_series: np.ndarray) -> np.ndarray: + """ + This function calculates the anomaly score for each sample within the time series. + + It also does some assertion regarding time series length. + + :param time_series: 1D array containing the time series to be scored + :return: anomaly score + """ + + # check the dimensions of the input array + assert time_series.ndim > 1, "Time series needs to be an N-D array. Currently it is 1-D." + + # compute the starting point of the scoring (past and future hankel need to fit) + starting_point = self.window_length + self.n_windows + self.lag + assert starting_point < time_series.shape[0], "The time series is too short to score any points." + + # scale the time series (or just copy it if already scaled) + time_series = time_series.copy() + if self.scale: + for idx in range(time_series.shape[1]): + time_series[:, idx] = normalization.min_max_scaling(time_series[:, idx], 1, 2, inplace=True) + + # get the different methods + scoring_function = self.methods[self.method] + return _transform(time_series, starting_point, self.window_length, self.n_windows, self.lag, + self.scoring_step, scoring_function, self.use_fast_hankel) + + +def _transform(time_series: np.ndarray, start_idx: int, window_length: int, n_windows: int, lag: int, scoring_step: int, + scoring_function: Callable, use_fast_hankel: bool) -> np.ndarray: + """ + Compute heavy and hopefully jit compilable score computation for the SST method. It does not do any parameter + checking and can throw cryptic errors. It's only used for internal use as a private function. + + :param time_series: 1D array containing the time series to be scored + :param start_idx: integer defining the start sample index for the score computation + :param window_length: the size of the time series window for each column of the hankel matrix + :param n_windows: column number of the hankel matrix + """ + + # initialize a scoring array with no values yet + score = np.zeros((time_series.shape[0],)) + + # compute the offset + offset = (n_windows + lag) + + # iterate over all the values in the signal starting at start_idx computing the change point score + for idx in range(start_idx, time_series.shape[0], scoring_step): + + # compile the past hankel matrix (H1) + hankel_past = blg.BlockHankelRepresentation(time_series, idx - lag, window_length, n_windows, use_fast_hankel) + + # compile the future hankel matrix (H2) + hankel_future = blg.BlockHankelRepresentation(time_series, idx, window_length, n_windows, use_fast_hankel) + + # compute the score and save the returned feedback vector + score[idx-offset-scoring_step//2:idx-offset+(scoring_step+1)//2] = \ + scoring_function(np.concatenate((hankel_past, hankel_future), axis=1)) + + return score + + +def _main(): + """ + Internal quick testing function. + + :return: + """ + from time import time + # make synthetic step function + np.random.seed(123) + # synthetic (frequency change) + x0 = np.sin(2 * np.pi * 1 * np.linspace(0, 10, 3000)) + x1 = np.sin(2 * np.pi * 2 * np.linspace(0, 10, 3000)) + x2 = np.sin(2 * np.pi * 8 * np.linspace(0, 10, 3000)) + x3 = np.sin(2 * np.pi * 4 * np.linspace(0, 10, 3000)) + x = np.hstack([x0, x1, x2, x3]) + x += np.random.rand(x.size) + + # create the method + messt_recognizer = MESST(70, method='rsvd', use_fast_hankel=False) + esst_recognizer = esst.ESST(70, method='rsvd', use_fast_hankel=False) + + # compute the score + start = time() + # score1 = messt_recognizer.transform(np.concatenate((x[..., None], x[..., None]), axis=1)) + score1 = messt_recognizer.transform(x[..., None]) + print(f'Computation for {len(x)} signal values took {time()-start} s.') + + # check for similarity with the esst + score2 = esst_recognizer.transform(x) + print(score1[1000:1020]) + print(score2[1000:1020]) + print(score1[1000:1020]-score2[1000:1020]) + print(np.allclose(score1, score2)) + print(np.corrcoef(score1, score2)[0, 1]) + + +if __name__ == '__main__': + _main() diff --git a/changepoynt/algorithms/msst.py b/changepoynt/algorithms/msst.py new file mode 100644 index 0000000..a79ac52 --- /dev/null +++ b/changepoynt/algorithms/msst.py @@ -0,0 +1,287 @@ +from typing import Callable +from functools import partial + +import numpy as np + +from changepoynt.utils import block_linalg as blg +from changepoynt.utils import normalization +import changepoynt.algorithms.sst as cpsst +from changepoynt.algorithms.base_algorithm import Algorithm + + +class MSST(Algorithm): + """ + This class implements all the utility and functionality necessary to compute the Multivariate Singular Spectrum + Transformation (MSST). + + This algorithm is an extension of the SST change point detection algorithm as described in: + + [1] + Idé, Tsuyoshi, and Keisuke Inoue. + "Knowledge discovery from heterogeneous dynamic systems using change-point correlations." + Proceedings of the 2005 SIAM international conference on data mining. + Society for Industrial and Applied Mathematics, 2005. + + [2] + Idé, Tsuyoshi, and Koji Tsuda. + "Change-point detection using krylov subspace learning." + Proceedings of the 2007 SIAM International Conference on Data Mining. + Society for Industrial and Applied Mathematics, 2007. + + It is basically a merger of techniques from [1, 2] and: + [3] + Alanqary, Arwa, Abdullah Alomar, and Devavrat Shah. + "Change point detection via multivariate singular spectrum analysis." + Advances in Neural Information Processing Systems 34 (2021): 23218-23230. + + It also uses a technique called randomized singular value decomposition which is surveyed and described in + + [4] + Halko, Nathan, Per-Gunnar Martinsson, and Joel A. Tropp. + "Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions." + SIAM review 53.2 (2011): 217-288. + + For the option fast_hankel=True it uses and algorithm based on + + [5] + L. Weber and R. Lenz. + "Accelerating Singular Spectrum Transformation for Scalable Change Point Detection," + in IEEE Access, Volume 11, 2025. + doi: 10.1109/ACCESS.2025.3640386. + + There will be a parameter specifying, whether to use the implicit krylov approximation from [2]. This + significantly speeds up the computation but can reduce accuracy as the "eigensignals" are only approximated + indirectly using a krylov subspace with the possible change point signal as seed. This method also deploys + a feedback of the dominant "eigensignal" as seed into the power method (dominant eigenvector estimation) + with additive gaussian noise. + + !NOTE!: + Most computational heavy functions are implemented as standalone functions, even if they require instance variables + as this enables us to use jit compiled code provided by the numba compiler for faster calculations. + + Also, we decided to do type hinting but not type checking as it requires too much boilerplate code. We recommend + to input the specified types as the program will break in unforeseen ways otherwise. + """ + + def __init__(self, window_length: int, n_windows: int = None, lag: int = None, rank: int = 5, scale: bool = True, + method: str = 'ika', lanczos_rank: int = None, random_rank: int = None, + feedback_noise_level: float = 1e-3, scoring_step: int = 1, use_fast_hankel: bool = False) -> None: + """ + Initializing the Singular Spectrum Transformation (SST) requires setting a lot of parameters. See the parameters + explanation for some intuition into the right choices. Currently, there are two SST methods from [1] and [2] + available for use. + + :param window_length: This specifies the length of the time series (in samples), which will be used to extract + the representative "eigensignals" to be compared before and after the lag. The window length should be big + enough to cover any wanted patterns (e.g., bigger than the periodicity of periodic signals). + + :param n_windows: This specifies the number of consecutive time windows used to extract the "eigensignals" from + the given time series. It should be big enough to cover different parts of the target behavior. If one does not + specify this value, we use the rule of thumb and take as many time windows as you specified the length of the + window (rule of thumb). + + :param lag: This parameter specifies the distance of the comparison for behaviors. In easy terms it tells the + algorithms how far it should look into the future to find change in behavior to the current signal. If you do + not specify this parameter, we use a rule of thumb and look ahead half of the window length you specified to + cover the behavior. + + :param rank: This parameter specifies the amount of "eigensignals" which will be used to measure the + dissimilarity of the signal in the future behavior. As a rule of thumb, we take the five most dominant + "eigensignals" if you do not specify otherwise. + + :param scale: Due to numeric stability, we REALLY RECOMMEND scaling the signal into a restricted value range. + Per default, we use a min max scaling to ensure a restricted value range. In the presence of extreme outliers, + this could cause problems, as the signal will be squished. + + :param method: Currently, "svd" [1], "ika" [2] and "rsvd" are available. ika corresponds to IKA-SST in [2] which + will speed up the computation significantly. + + :param lanczos_rank: In order to use the implicit approximation of "eigensignals" by using "ika" [2] method, + one needs to decide the rank of the implicit approximation of each "eigensignal". As a rule of thumb, we + determine this value as being twice the amount of specified "eigensignals" to span the subspace + (parameter rank). This is also recommended by the author of IKA-SST in [2]. + + :param random_rank: In order to use the randomized singular value decomposition, one needs to provide a + randomized rank (size of second matrix dimension for the randomized matrix) as specified in [3]. The lower + this value, the faster the computation but the higher the error (as the approximation gets worse). + + :param feedback_noise_level: This specifies the amplitude of additive white gaussian noise added to the dominant + "eigensignal" of the future behavior when shifting forward. This idea is noted in [2] and initializes + the seed of the power method for dominant eigenvector estimation with the precious dominant eigenvector + plus the noise level specified here. The noise level should just be a small fraction of the value range + of the signal. + + :param scoring_step: the distance between scoring steps in samples (e.g., 2 would half the computation). + + :param use_fast_hankel: Use the O(N*logN) version for the decomposition. This is likely to have large + speed improvements for window sizes > 150 + """ + + # save the specified parameters into instance variables + self.window_length = window_length + self.n_windows = n_windows + self.lag = lag + self.rank = rank + self.scale = scale + self.method = method + self.lanczos_rank = lanczos_rank + self.random_rank = random_rank + self.noise = feedback_noise_level + self.scoring_step = scoring_step + self.use_fast_hankel = use_fast_hankel + + # set some default values when they have not been specified + if self.n_windows is None: + self.n_windows = self.window_length + if self.lag is None: + # rule of thumb + self.lag = max(self.n_windows // 3, 1) + if self.lanczos_rank is None: + # make rank even and multiply by two just as specified in [2] + self.lanczos_rank = self.rank * 2 - (self.rank & 1) + if self.random_rank is None: + # compute the rank as specified in [3] and + # https://scikit-learn.org/stable/modules/generated/sklearn.utils.extmath.randomized_svd.html + self.random_rank = min(self.rank + 10, self.window_length, self.n_windows) + + # specify the methods and their corresponding functions as lambda functions expecting only the + # 1) future hankel matrix, + # 2) the current hankel matrix and + # 3) the feedback vector (e.g. for dominant eigenvector feedback) + # all other parameter should be specified as values in the partial lambda function + self.methods = {'ika': partial(cpsst._implicit_krylov_approximation, + rank=self.rank, + lanczos_rank=self.lanczos_rank), + 'rsvd': partial(cpsst._random_singular_value_decomposition, + rank=self.rank, + randomized_rank=self.random_rank), + 'weighted': partial(cpsst._weighted_random_singular_value_decomposition, + rank=self.rank, + randomized_rank=self.random_rank), + 'symmetric': partial(cpsst._symmetric_random_singular_value_decomposition, + rank=self.rank, + randomized_rank=self.random_rank), + } + if self.method not in self.methods: + raise ValueError(f'Method {self.method} not defined. Possible methods: {list(self.methods.keys())}.') + + # set up the methods we use for the construction of the hankel matrix (either it is the fft representation + # of the other one) + if use_fast_hankel and self.method not in ["rsvd", "ika", "weighted", "symmetric"]: + raise ValueError(f'{self.method} method is not defined with use_fast_hankel=True') + + def transform(self, time_series: np.ndarray) -> np.ndarray: + """ + This function calculates the anomaly score for each sample within the time series, starting from an initial + sample being the first sample of fitting the past hankel matrix (window_length + n_windows samples) and the new + future hankel matrix (+lag). Values before that will be zero + + It also does some assertions regarding the specified parameters and whether they fit the time series. + + This function builds the interface to the jit compiled static function used to iterate over the array. + + :param time_series: 1D array containing the time series to be scored + :return: anomaly score + """ + + # check the dimensions of the input array + assert time_series.ndim > 1, "Time series needs to be an N-D array. Currently it is 1-D." + + # compute the starting point of the scoring (past and future hankel need to fit) + starting_point = self.window_length + self.n_windows + self.lag + assert starting_point < time_series.shape[0], "The time series is too short to score any points." + + # scale the time series (or just copy it if already scaled) + time_series = time_series.copy() + if self.scale: + for idx in range(time_series.shape[1]): + time_series[:, idx] = normalization.min_max_scaling(time_series[:, idx], 1, 2, inplace=True) + + # get the changepoint scorer from the different methods + scoring_function = self.methods[self.method] + + # start the scaling itself by calling the jit compiled staticmethod and return the result + score = _transform(time_series=time_series, start_idx=starting_point, window_length=self.window_length, + n_windows=self.n_windows, lag=self.lag, scoring_step=self.scoring_step, + scoring_function=scoring_function, use_fast_hankel=True) + return score + + +def _transform(time_series: np.ndarray, start_idx: int, window_length: int, n_windows: int, lag: int, scoring_step: int, + scoring_function: Callable, use_fast_hankel: bool) -> np.ndarray: + """ + Compute heavy and hopefully jit compilable score computation for the SST method. It does not do any parameter + checking and can throw cryptic errors. It's only used for internal use as a private function. + + :param time_series: 1D array containing the time series to be scored + :param start_idx: integer defining the start sample index for the score computation + :param window_length: the size of the time series window for each column of the hankel matrix + :param n_windows: number of columns in the hankel matrix + :param lag: sample distance between future and past hankel matrix + :param scoring_step: the distance between scoring steps in samples. + :param scoring_function: the function that is called every step to assess a scalar change point score + """ + + # create initial vector for ika method with feedback dominant eigenvector as proposed in [2] + # with a norm of one + x0 = np.random.rand(window_length*time_series.shape[1])[:, None] + x0 /= np.linalg.norm(x0) + + # initialize a scoring array with no values yet + score = np.zeros((time_series.shape[0],)) + + # make an offset for the data construction + offset = n_windows // 2 + lag + + # iterate over all the values in the signal starting at start_idx computing the change point score + for idx in range(start_idx, time_series.shape[0], scoring_step): + + # compile the past hankel matrix (H1) + hankel_past = blg.BlockHankelRepresentation(time_series, idx - lag, window_length, n_windows, use_fast_hankel) + + # compile the future hankel matrix (H2) + hankel_future = blg.BlockHankelRepresentation(time_series, idx, window_length, n_windows, use_fast_hankel) + + # compute the score and save the returned feedback vector + score[idx - offset - scoring_step // 2:idx - offset + (scoring_step + 1) // 2], x1 = \ + scoring_function(hankel_past, hankel_future, x0) + + # add noise to the dominant eigenvector and normalize it again + x0 = x1 + 1e-3 * np.random.rand(x0.shape[0])[:, None] + x0 /= np.linalg.norm(x0) + + return score + + +def main(): + """This function is not intended for users, but for quick testing during development.""" + from time import time + + # make synthetic step function + np.random.seed(123) + length = 300 + x = np.hstack([1 * np.ones(length) + np.random.rand(length) * 1, + 3 * np.ones(length) + np.random.rand(length) * 2, + 5 * np.ones(length) + np.random.rand(length) * 1.5]) + x += np.random.rand(x.size) + x = x[..., None] + + # create the sst method + ika_sst = MSST(31, method='ika', use_fast_hankel=True) + rsvd_sst = MSST(31, method='rsvd', use_fast_hankel=True) + weighted_sst = MSST(31, method='weighted', use_fast_hankel=True) + + # compare with original sst + orig_ika = cpsst.SST(31, method='ika', use_fast_hankel=True).transform(x[:, 0]) + print('Comparison to original score', np.corrcoef(orig_ika, ika_sst.transform(x))[0, 1]) + + # make the scoring + start = time() + ika_sst.transform(x) + print(time() - start) + rsvd_sst.transform(x) + weighted_sst.transform(x) + + +if __name__ == '__main__': + main() diff --git a/changepoynt/algorithms/sst.py b/changepoynt/algorithms/sst.py index 783cdef..b7693dc 100644 --- a/changepoynt/algorithms/sst.py +++ b/changepoynt/algorithms/sst.py @@ -143,6 +143,11 @@ def __init__(self, window_length: int, n_windows: int = None, lag: int = None, r # 2) the current hankel matrix and # 3) the feedback vector (e.g. for dominant eigenvector feedback) # all other parameter should be specified as values in the partial lambda function + # + # We recommend using: + # 1) rsvd + # 2) ika if speed is a concern + # 3) weighted if noise is low self.methods = {'ika': partial(_implicit_krylov_approximation, rank=self.rank, lanczos_rank=self.lanczos_rank), @@ -161,13 +166,16 @@ def __init__(self, window_length: int, n_windows: int = None, lag: int = None, r 'weighted': partial(_weighted_random_singular_value_decomposition, rank=self.rank, randomized_rank=self.random_rank), + 'symmetric': partial(_symmetric_random_singular_value_decomposition, + rank=self.rank, + randomized_rank=self.random_rank), } if self.method not in self.methods: raise ValueError(f'Method {self.method} not defined. Possible methods: {list(self.methods.keys())}.') # set up the methods we use for the construction of the hankel matrix (either it is the fft representation # of the other one) - if use_fast_hankel and self.method not in ["rsvd", "ika", "weighted"]: + if use_fast_hankel and self.method not in ["rsvd", "ika", "weighted", "symmetric"]: raise ValueError(f'{self.method} method is not defined with use_fast_hankel=True') self.hankel_construction = { False: lg.compile_hankel, @@ -486,6 +494,45 @@ def _weighted_random_singular_value_decomposition(hankel_past: np.ndarray, hanke return score, x0 +def _symmetric_random_singular_value_decomposition(hankel_past: np.ndarray, hankel_future: np.ndarray, x0: np.ndarray, + rank: int, randomized_rank: int): + """ + This function implements the idea of + + Idé, Tsuyoshi, and Keisuke Inoue. + "Knowledge discovery from heterogeneous dynamic systems using change-point correlations." + Proceedings of the 2005 SIAM international conference on data mining. + Society for Industrial and Applied Mathematics, 2005. + + but uses the randomized svd decomposition by + + Halko, Nathan, Per-Gunnar Martinsson, and Joel A. Tropp. + "Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions." + SIAM review 53.2 (2011): 217-288. + + !! It uses more than just the future vector. Instead, it computes a weighted score using multiple vectors + + :param hankel_past: the hankel matrix H1 before the change point + :param hankel_future: the hankel matrix H2 after the change point + :param x0: the initialization value for the power method applied to H2 to find the dominant eigenvector + :param rank: the amount of (approximated) eigenvectors as subspace of H1 + :param randomized_rank: the rank of the approximation used to construct the noise matrix + :return: the change point score, the new dominant eigenvector of H2 for the feedback into the next H2 + """ + + # compute the biggest eigenvector of the hankel matrix after the possible change point (h2) + singvecs_future, _, _ = lg.randomized_hankel_svd(hankel_future, rank, oversampling_p=randomized_rank-rank) + + # compute the eigenvectors of the past hankel matrix + singvecs_past, _, _ = lg.randomized_hankel_svd(hankel_past, rank, oversampling_p=randomized_rank-rank) + + # compute the forward score + forward_score = 1-np.sum(np.square(singvecs_past[:, :rank].T @ singvecs_future[:, 0])) + backward_score = 1-np.sum(np.square(singvecs_future[:, :rank].T @ singvecs_past[:, 0])) + score = (forward_score + backward_score) / 2 + return score, x0 + + def _naive_singular_value_decomposition(hankel_past: np.ndarray, hankel_future: np.ndarray, x0: np.ndarray, rank: int): """ diff --git a/changepoynt/utils/block_linalg.py b/changepoynt/utils/block_linalg.py new file mode 100644 index 0000000..c2e645a --- /dev/null +++ b/changepoynt/utils/block_linalg.py @@ -0,0 +1,1050 @@ +import typing +import abc + +import numpy as np +import scipy.fft as spfft + + +def block_hankel_left_matmat_fft(hankel_fft: np.ndarray, other_matrix: np.ndarray, fft_length: int, + window_length: int, window_number: int): + """ + Compute A @ H using FFT, where H is a block Hankel matrix. + + H has block structure: + + H[i, j] = B[i + j] + + + Definitions + ---------- + p, q : the shape of the Hankel matrix (p rows, q columns) + m, n : the shape of the blocks in the Hankel matrix (each block has m rows, n columns) + + Parameters + ---------- + + hankel_fft : ndarray, shape (p + q - 1, m, n) + The fft transformed (along axis 0) of the Block sequence. B[k] is an m x n matrix block. + + other_matrix : ndarray + The matrix that we multiply from the right (other_matrix @ hankel_fft). + Either: + - dense form: shape (s, p*m) + - block form: shape (p, s, m) + + fft_length : int + The length of the fft (with padding). We recommend using scipy.fft.next_fast_length + + window_length : int + The number of blocks per row of the Hankel matrix + + window_number : int + The number of blocks per column of the Hankel matrix + + Returns + ------- + result_matrix : ndarray + Outcome of other_matrix @ hankel_fft + If A was dense, returns dense shape (s, q*n). + If A was block-form, returns block shape (q, s, n). + """ + + # get the dimensions of the hankel matrix + _, m, n = hankel_fft.shape + p = window_length + q = window_number + + return_dense = False + + if other_matrix.ndim == 2: + # other_matrix is dense: shape (s, p*m) + s, cols = other_matrix.shape + + if cols != p*m: + raise ValueError("A has incompatible number of columns.") + + other_matrix_blocks = other_matrix.reshape(s, p, m).transpose(1, 0, 2) + return_dense = True + + elif other_matrix.ndim == 3: + # other_matrix is already block-form: shape (p, s, m) + p_x, s, m_x = other_matrix.shape + + if m_x != m: + raise ValueError("Block dimensions do not match.") + + other_matrix_blocks = other_matrix + + else: + raise ValueError("A must have shape (s, p*m) or (p, s, m).") + + if q <= 0: + raise ValueError("Need len(B) >= p.") + + # For left multiplication: + # + # result_matrix_j = sum_i fft_hankel_i other_matrix_{i+j} + # + # This is obtained from convolution of reversed A with B. + other_matrix_rev = other_matrix_blocks[::-1] + + other_matrix_fft = spfft.rfft(other_matrix_rev, n=fft_length, axis=0) # shape (fft_len, s, m) + + # Frequency-wise matrix-matrix multiplication: + # + # result_matrix_fft[ell] = other_matrix_fft[ell] @ fft_hankel[ell] + # + # Shapes: + # fft_hankel[ell] : (s, m) + # other_matrix_fft[ell] : (m, n) + # result_matrix_fft[ell] : (s, n) + result_matrix_fft = other_matrix_fft @ hankel_fft # shape (fft_len, s, n) + + conv = spfft.irfft(result_matrix_fft, n=fft_length, axis=0) + + # Extract the desired Hankel result + result_matrix = conv[p - 1 : p - 1 + q] # shape (q, s, n) + + result_matrix = np.real_if_close(result_matrix) + + if return_dense: + return result_matrix.transpose(1, 0, 2).reshape(s, q * n) + + return result_matrix + + +def block_hankel_right_matmat_direct(hankel: np.ndarray, other_matrix: np.ndarray): + """ + Compute H @ X without FFT, using a vectorized einsum. + + H has block structure: + + H[i, j] = B[i + j] + + Definitions + ---------- + p, q : the shape of the Hankel matrix (p rows, q columns) + m, n : the shape of the blocks in the Hankel matrix (each block has m rows, n columns) + + Parameters + ---------- + hankel : ndarray, shape (p + q - 1, m, n) + Block sequence. hankel[k] is an m x n matrix block. + + other_matrix : ndarray + Either: + - dense form: shape (q*n, r) + - block form: shape (q, n, r) + + Returns + ------- + Y : ndarray + If X was dense, returns dense shape (p*m, r). + If X was block-form, returns block shape (p, m, r). + """ + + num_blocks, m, n = hankel.shape + return_dense = False + + if other_matrix.ndim == 2: + rows, r = other_matrix.shape + + if rows % n != 0: + raise ValueError("other_matrix_blocks has incompatible number of rows.") + + q = rows // n + other_matrix_blocks = other_matrix.reshape(q, n, r) + return_dense = True + + elif other_matrix.ndim == 3: + (q, n_other_matrix, r) = other_matrix.shape + + if n_other_matrix != n: + raise ValueError("Block dimensions do not match.") + + other_matrix_blocks = other_matrix + + else: + raise ValueError("other_matrix_blocks must have shape (q*n, r) or (q, n, r).") + + p = num_blocks - q + 1 + + if p <= 0: + raise ValueError("Need len(hankel) >= q.") + + # sliding_window_view gives shape (p, m, n, q) + hankel_windows = np.lib.stride_tricks.sliding_window_view(hankel, window_shape=q, axis=0) + + # Move the window axis so that: + # b_moved[i, j] = hankel[i + j] + # b_moved has shape (p, q, m, n) + hankel_windows_moved = np.moveaxis(hankel_windows, -1, 1) + + # result_matrix[i, a, c] = sum_{j,b} b_moved[i,j,a,b] * other_matrix_blocks_blocks[j,b,c] + result_matrix = np.einsum("ijmn,jnr->imr", hankel_windows_moved, other_matrix_blocks, optimize=True) + + if return_dense: + return result_matrix.reshape(p * m, r) + + return result_matrix + + +def block_hankel_right_matmat_fft(hankel_fft: np.ndarray, other_matrix: np.ndarray, fft_length: int, + window_length: int, window_number: int): + """ + Compute H @ X using FFT, where H is a block Hankel matrix. + + H has block structure: + + H[i, j] = B[i + j] + + + Definitions + ---------- + p, q : the shape of the Hankel matrix (p rows, q columns) + m, n : the shape of the blocks in the Hankel matrix (each block has m rows, n columns) + + Parameters + ---------- + hankel_fft : ndarray, shape (p + q - 1, m, n) + The fft transformed (along axis 0) of the Block sequence. B[k] is an m x n matrix block. + + other_matrix : ndarray + Either: + - dense form: shape (q*n, r) + - block form: shape (q, n, r) + + fft_length : int + The length of the fft (with padding). We recommend using scipy.fft.next_fast_length + + window_length : int + The number of blocks per row of the Hankel matrix + + window_number : int + The number of blocks per column of the Hankel matrix + + Returns + ------- + Y : ndarray + If X was dense, returns dense shape (p*m, r). + If X was block-form, returns block shape (p, m, r). + """ + _, m, n = hankel_fft.shape + p = window_length + q = window_number + + return_dense = False + + if other_matrix.ndim == 2: + # X is dense: shape (q*n, r) + rows, r = other_matrix.shape + + if rows != q*n: + raise ValueError("X has incompatible number of rows.") + + other_matrix_blocks = other_matrix.reshape(q, n, r) + return_dense = True + + elif other_matrix.ndim == 3: + # X is already block-form: shape (q, n, r) + q_x, n_x, r = other_matrix.shape + + if q_x != q or n_x != n: + raise ValueError("Block dimensions do not match.") + + other_matrix_blocks = other_matrix + + else: + raise ValueError("X must have shape (q*n, r) or (q, n, r).") + + if p <= 0: + raise ValueError("Need len(B) >= q.") + + # Hankel multiplication: + # + # Y_i = sum_j B_{i+j} X_j + # + # This becomes convolution of B with reversed X. + other_matrix_rev = other_matrix_blocks[::-1] + other_matrix_fft = spfft.rfft(other_matrix_rev, n=fft_length, axis=0) # shape (fft_len, n, r) + + # Frequency-wise matrix-matrix multiplication: + # + # FY[ell] = FB[ell] @ FX[ell] + # + # Shapes: + # FB[ell] : (m, n) + # FX[ell] : (n, r) + # FY[ell] : (m, r) + result_matrix_fft = hankel_fft @ other_matrix_fft # shape (fft_len, m, r) + + conv = spfft.irfft(result_matrix_fft, n=fft_length, axis=0) + + # Extract the desired Hankel result + result_matrix_blocks = conv[q - 1 : q - 1 + p] # shape (p, m, r) + + result_matrix_blocks = np.real_if_close(result_matrix_blocks) + + if return_dense: + return result_matrix_blocks.reshape(p * m, r) + + return result_matrix_blocks + + +def block_hankel_left_matmat_direct(hankel: np.ndarray, other_matrix: np.ndarray): + """ + Compute A @ H without FFT, using a vectorized einsum. + + H has block structure: + + H[i, j] = B[i + j] + + Definitions + ---------- + p, q : the shape of the Hankel matrix (p rows, q columns) + m, n : the shape of the blocks in the Hankel matrix (each block has m rows, n columns) + + Parameters + ---------- + hankel : ndarray, shape (p + q - 1, m, n) + Block sequence. hankel[k] is an m x n matrix block. + + other_matrix : ndarray + Matrix to multiply by + Either: + - dense form: shape (s, p*m) + - block form: shape (p, s, m) + + Returns + ------- + Z : ndarray + If A was dense, returns dense shape (s, q*n). + If A was block-form, returns block shape (q, s, n). + """ + + num_blocks, m, n = hankel.shape + return_dense = False + + if other_matrix.ndim == 2: + s, cols = other_matrix.shape + + if cols % m != 0: + raise ValueError("A has incompatible number of columns.") + + p = cols // m + other_matrix_blocks = other_matrix.reshape(s, p, m).transpose(1, 0, 2) + return_dense = True + + elif other_matrix.ndim == 3: + p, s, m_A = other_matrix.shape + + if m_A != m: + raise ValueError("Block dimensions do not match.") + + other_matrix_blocks = other_matrix + + else: + raise ValueError("A must have shape (s, p*m) or (p, s, m).") + + q = num_blocks - p + 1 + + if q <= 0: + raise ValueError("Need len(B) >= p.") + + # sliding_window_view gives shape (q, m, n, p) + hankel_windows = np.lib.stride_tricks.sliding_window_view(hankel, window_shape=p, axis=0) + + # Move the window axis so that: + # BH[j, i] = B[j + i] + # BH has shape (q, p, m, n) + hankel_windows_moved = np.moveaxis(hankel_windows, -1, 1) + + # Z[j, s, n] = sum_{i,m} A_blocks[i,s,m] * BH[j,i,m,n] + result_matrix_blocks = np.einsum("ism,jimn->jsn", other_matrix_blocks, hankel_windows_moved, optimize=True) + + if return_dense: + return result_matrix_blocks.transpose(1, 0, 2).reshape(s, q * n) + + return result_matrix_blocks + + +def get_block_hankel_representation(time_series: np.ndarray, end_index: int, window_length: int, window_number: int) -> [np.ndarray, int]: + + # make some checks for the signal + if time_series.ndim != 2: + raise ValueError(f'Input time series has to have two dimensions. Currently: {time_series.ndim=} != 2.') + + # extract the block hankel matrix representation + # we expect it to have shape (window_number+window_length-1, m, n) where m is the number of time series + # in the multivariate time series and n is always one in our case + size = window_number+window_length-1 + return time_series[end_index-(window_number+window_length-1):end_index, :][..., None], size + + +class BlockHankel: + """ + This class is the base class for BlockHankel matrices. + It mainly serves as the interface to capture numpy matrix multiplications (@) using the functions + __matmul__ (if a BlockHankel matrix is on the left side of the @ operator) + and __array_ufunc__ (if a BlockHankel is on the right side of the @ operator) + + If these functions are called, the input is given to the abstract methods + multiply_other_from_left + and + multiply_other_from_right + which we expect the subclasses to implement in detail + """ + + # denote that we will have a shape + shape: tuple[int, ...] + + @abc.abstractmethod + def materialize(self) -> np.ndarray: + raise NotImplementedError + + @abc.abstractmethod + def multiply_other_from_right(self, other_maxtrix: np.ndarray) -> np.ndarray: + """Subclasses should implement this method.""" + raise NotImplementedError + + @abc.abstractmethod + def multiply_other_from_left(self, other_maxtrix: np.ndarray) -> np.ndarray: + """Subclasses should implement this method.""" + raise NotImplementedError + + @abc.abstractmethod + def copy(self, deep: bool = False) -> "BlockHankel": + """Subclasses should implement this method.""" + raise NotImplementedError + + @property + @abc.abstractmethod + def T(self,) -> "BlockHankel": + """Subclasses should implement this method. The method should transpose the matrix.""" + raise NotImplementedError + + def __matmul__(self, other) -> typing.Union[np.ndarray, "BlockHankel"]: + """ + This function is called by numpy if we use an instance of the current class on the left side of + a matrix multiplication operator (@) + + :param other: the other matrix + :return: + """ + # right side matmul + if isinstance(other, np.ndarray): + return self.multiply_other_from_right(other) + if isinstance(other, BlockHankel): + return BlockHankelProductRepresentation(self, other) + else: + raise ValueError('At this point matmul is only with other matrices.') + + def __array_ufunc__(self, ufunc, method, *args, **kwargs) -> typing.Union[np.ndarray, "BlockHankel"]: + """ + This function is called by numpy, if an instance of the current class on the left side of + an operator where the left side is any numpy object. + + Basically, numpy first checks if the left side of an operator has defined an operation for + both objects on the right and left. If this fails, it checks the object on the right. + + :param ufunc: the function that is called on the objects + :param method: what type of invocation it was + :param args: the arguments of the function + :param kwargs: the keyword arguments of the function + :return: + """ + if method != '__call__' or len(kwargs) > 0 or len(args) != 2: + return NotImplemented + if ufunc == np.matmul: + # left side matmul + if isinstance(args[0], np.ndarray) and args[1] is self: + return self.multiply_other_from_left(args[0]) + else: + return NotImplemented + else: + return NotImplemented + + def __array_function__(self, func, types, args, kwargs) -> "MultilevelBlockHankelRepresentation": + """ + We currently only support concatenation. + https://numpy.org/neps/nep-0018-array-function-protocol.html + """ + if func != np.concatenate: + return NotImplemented + if not all(issubclass(t, BlockHankel) for t in types) or not issubclass(self.__class__, BlockHankel): + return NotImplemented + + # check the keyword arguments + axis = kwargs.get('axis', 0) + if axis not in [0, 1]: + raise ValueError('Concatenation only in the first two dimensions.') + + # check that we only received one arg + if len(args) > 1: + raise ValueError('We only accept one non keyword argument.') + + # create the list of matrices + return MultilevelBlockHankelRepresentation(args[0], axis) + + +class BlockHankelRepresentation(BlockHankel): + """ + This matrix represents a Block Hankel matrix + for large matrices and fast_hankel=True it makes matrix multiplication faster: + + See our paper for more details: + Efficient Hankel Matrix Decomposition for Changepoint Detection + Lucas Weber, Richard Lenz 2024 + """ + + def __init__(self, time_series: np.ndarray, end_index: int, window_length: int, window_number: int, + fast_hankel: bool): + + # save the parameters that are necessary for later computations + self.window_length = window_length + self.window_number = window_number + self.fast_hankel = fast_hankel + + # get the hankel matrix block representation + hankel, size = get_block_hankel_representation(time_series, end_index, window_length, window_number) + self.size = size + + # create the representation and save it into the class + if self.fast_hankel: + + # get the next fast fft length for efficient fft + self.fft_length = spfft.next_fast_len(hankel.shape[0], real=True) + + # make the fft along the individual time series axis + hankel = spfft.rfft(hankel, n=self.fft_length, axis=0) + + # save the hankel matrix and the shape + self.hankel = hankel + self.shape = (self.window_length*time_series.shape[1], self.window_number) + + def copy(self, deep: bool = False) -> "BlockHankelRepresentation": + new = self.__class__.__new__(self.__class__) + + new.window_length = self.window_length + new.window_number = self.window_number + new.fast_hankel = self.fast_hankel + new.shape = self.shape + new.size = self.size + + if hasattr(self, "fft_length"): + new.fft_length = self.fft_length + + if deep: + new.hankel = self.hankel.copy() + else: + new.hankel = self.hankel + + return new + + def multiply_other_from_right(self, other_matrix: np.ndarray) -> np.ndarray: + # check the dimensions + if other_matrix.shape[0] != self.shape[1]: + raise ValueError(f'matmul: Right matrix has a mismatch in its core dimension 0 ' + f'(size {other_matrix.shape[0]=} is different from {self.shape[1]=})') + + # make the product based on whether we use the fast Hankel option + if self.fast_hankel: + return block_hankel_right_matmat_fft(self.hankel, other_matrix, self.fft_length, self.window_length, self.window_number) + else: + return block_hankel_right_matmat_direct(self.hankel, other_matrix) + + def multiply_other_from_left(self, other_matrix: np.ndarray) -> np.ndarray: + + # check the dimensions + if other_matrix.shape[1] != self.shape[0]: + raise ValueError(f'matmul: Left matrix has a mismatch in its core dimension 0 ' + f'(size {other_matrix.shape[1]=} is different from {self.shape[0]=})') + + # make the product based on whether we use the fast Hankel option + if self.fast_hankel: + return block_hankel_left_matmat_fft(self.hankel, other_matrix, self.fft_length, self.window_length, self.window_number) + else: + return block_hankel_left_matmat_direct(self.hankel, other_matrix) + + @property + def T(self) -> 'BlockHankelRepresentation': + + # make a copy of the object + new = self.copy(deep=False) + + # switch the shape + new.shape = tuple(reversed(new.shape)) + + # switch the number of windows and window length + new.window_length, new.window_number = new.window_number, new.window_length + + # switch the last dimensions of the hankel matrix + new.hankel = new.hankel.transpose((0, 2, 1)) + return new + + def materialize(self) -> np.ndarray: + """ + Explicitly build the full block Hankel matrix H for test multiplications. + + Definitions + ---------- + p, q : the shape of the Hankel matrix (p rows, q columns) + m, n : the shape of the blocks in the Hankel matrix (each block has m rows, n columns) + + B has shape (p + q - 1, m, n). + H has shape (p*m, q*n). + """ + hankel = self.hankel + if self.fast_hankel: + hankel = np.real_if_close(spfft.irfft(self.hankel, n=self.fft_length, axis=0))[:self.size,...] + _, m, n = hankel.shape + + hankel_matrix_full = np.zeros(self.shape) + for i in range(self.window_length): + for j in range(self.window_number): + hankel_matrix_full[i * m:(i + 1) * m, j * n:(j + 1) * n] = hankel[i + j] + + return hankel_matrix_full + + +class MultilevelBlockHankelRepresentation(BlockHankel): + """ + This class implements multiplication with a multilevel Block Hankel matrix, i.e., a matrix which can be divided into + multiple parts that each have Hankel structure. + + We expect to receive a list of: + - Hankel matrices as BlockHankelRepresentation objects + - The concatenation axis (currently only 0 and 1 are supported) + """ + + def __init__(self, matrix_tuples: list[BlockHankelRepresentation], axis: int): + + # check the input + self.shape, self.individual_shapes = self.check_input(matrix_tuples, axis) + + # extract the matrices + self.matrices = matrix_tuples + self.axis = axis + + @staticmethod + def check_input(input_list: list[BlockHankelRepresentation], axis: int) -> (int, int): + + # check whether the list is empty + if len(input_list) < 1: + raise ValueError('You did not input any Block Hankel matrices.') + + # check that the axis is either 0 or 1 + if axis not in (0, 1): + raise ValueError(f'Concatenation only in the first two dimensions. Currently specified {axis=}.') + + # save the shapes + shapes = [(0, 0)]*len(input_list) + + # check that we received only BlockHankelRepresentation matrices and the shapes match + for idx, ele in enumerate(input_list): + + # check for the type of the class + if not isinstance(ele, BlockHankelRepresentation): + raise TypeError(f'We can only concatenate BlockHankelRepresentations in this class. Got {type(ele)}.') + + # save the dimension that we have to check + shapes[idx] = ele.shape + + # decide which dimension we have to check for uniqueness + # e.g., if we concatenate along axis 0, all axis 1 sizes have to match + checkdim = 0 if axis == 1 else 1 + if len(set(ele[checkdim] for ele in shapes)) != 1: + raise ValueError(f'We cannot concatenate BlockHankelRepresentations in this class due to different shapes in dimension {checkdim}: {shapes}.') + + if axis == 1: + shape = (shapes[0][0], sum(ele[1] for ele in shapes)) + elif axis == 0: + shape = (sum(ele[0] for ele in shapes), shapes[0][1]) + else: + raise ValueError(f'{axis=} is invalid.') + return shape, shapes + + def copy(self, deep: bool = False) -> "MultilevelBlockHankelRepresentation": + new = self.__class__.__new__(self.__class__) + + if deep: + new.matrices = [ele.copy(deep=deep) for ele in self.matrices] + else: + new.matrices = self.matrices + new.axis = self.axis + new.shape = self.shape + + if deep: + new.individual_shapes = [ele.copy() for ele in self.individual_shapes] + else: + new.individual_shapes = self.individual_shapes + return new + + @property + def T(self) -> 'MultilevelBlockHankelRepresentation': + """ + Compute the transposed. + :return: + """ + new = self.copy() + # switch the shape + new.shape = tuple(reversed(new.shape)) + + # go through all of our matrices and transpose them + new.matrices = [matrix.T for matrix in new.matrices] + + # update the shapes + new.individual_shapes = [ele.shape for ele in new.matrices] + + # switch the axis + new.axis = 1 if new.axis == 0 else 0 + + return new + + def multiply_other_from_right(self, other: np.ndarray): + + # we need to create a target matrix + result = np.zeros((self.shape[0], other.shape[1])) + + # check the matrix shapes + if self.shape[1] != other.shape[0]: + raise ValueError(f'Shape missmatch: {self.shape=} multiplied with {other.shape=} is not possible.') + + # choose depending on the axis of concatenation + starting_point = 0 + if self.axis == 0: # matrices are stacked ontop of each other + for idx, matrix in enumerate(self.matrices): + end_point = starting_point+self.individual_shapes[idx][0] + result[starting_point:end_point, :] = matrix @ other + starting_point = end_point + if self.axis == 1: # matrices are stacked next to each other + for idx, matrix in enumerate(self.matrices): + end_point = starting_point + starting_point+self.individual_shapes[idx][1] + result += matrix @ other[starting_point:end_point, :] + starting_point = end_point + + return result + + def multiply_other_from_left(self, other: np.ndarray): + + # we need to create a target matrix + result = np.zeros((other.shape[0], self.shape[1])) + + # check the matrix shapes + if self.shape[0] != other.shape[1]: + raise ValueError(f'Shape missmatch: {other.shape=} multiplied with {self.shape=} is not possible.') + + # choose depending on the axis of concatenation + starting_point = 0 + if self.axis == 0: # matrices are stacked ontop of each other + for idx, matrix in enumerate(self.matrices): + end_point = starting_point + starting_point + self.individual_shapes[idx][0] + result += other[:, starting_point:end_point] @ matrix + starting_point = end_point + if self.axis == 1: # matrices are stacked next to each other + for idx, matrix in enumerate(self.matrices): + end_point = starting_point + self.individual_shapes[idx][1] + result[:, starting_point:end_point] = other @ matrix + starting_point = end_point + return result + + def materialize(self) -> np.ndarray: + """ + Explicitly build the full multilevel block Hankel matrix H for test multiplications. + """ + + return np.concatenate(tuple(mat.materialize()for mat in self.matrices), axis=self.axis) + + +class BlockHankelProductRepresentation(BlockHankel): + """ + This class represents a matrix-matrix product of two hankel matrices. + H @ H + Instead of computing the product and materializing it, we keep both matrices so we have hankel property for + products. + + See Section III-G in: + L. Weber and R. Lenz, + "Accelerating Singular Spectrum Transformation for Scalable Change Point Detection," + in IEEE Access, vol. 13, pp. 213556-213577, 2025, doi: 10.1109/ACCESS.2025.3640386. + for details. + """ + + def __init__(self, first_hankel: BlockHankel, second_hankel: BlockHankel): + + # save the matrices for keeping + self.first_hankel = first_hankel + self.second_hankel = second_hankel + + # check that the dimensions match + if self.first_hankel.shape[1] != self.second_hankel.shape[0]: + raise ValueError(f'matmul: Operant shape missmatch: {self.first_hankel.shape} @ {self.second_hankel.shape} missmatch in inner dimensions.') + + # save the shape + self.shape = (self.first_hankel.shape[0], self.second_hankel.shape[1]) + + def copy(self, deep: bool = False) -> "BlockHankelProductRepresentation": + new = self.__class__.__new__(self.__class__) + + # copy the matrices + new.first_hankel = self.first_hankel.copy(deep=deep) + new.second_hankel = self.second_hankel.copy(deep=deep) + + # copy the shape + new.shape = self.shape + return new + + @property + def T(self): + + # make a copy of myself + new = self.copy() + + # switch matrix order and transpose the matrices + new.first_hankel, new.second_hankel = new.second_hankel.T, new.first_hankel.T + + # update the shape + new.shape = tuple(reversed(new.shape)) + return new + + def multiply_other_from_right(self, other: np.ndarray) -> np.ndarray: + return self.first_hankel @ (self.second_hankel @ other) + + def multiply_other_from_left(self, other: np.ndarray) -> np.ndarray: + return (other @ self.first_hankel) @ self.second_hankel + + def materialize(self) -> np.ndarray: + # we need to materialize the second otherwise, we will again only create a + # BlockHankelProductRepresentation (the current class) + # + # materializing the second makes this a product of a BlockHankel with a np.ndarray + # instead of BlockHankel with BlockHankel + return self.first_hankel @ self.second_hankel.materialize() + + +def check_matrix_operations(hankel_typed_matrix: BlockHankel, materialized_hankel: np.ndarray, other_size: int = 10): + + # print the checks + print() + print(f"## Starting checks for matrix of type {hankel_typed_matrix.__class__.__name__}. ##") + print('################################################') + + # check the shape + assert hankel_typed_matrix.shape == materialized_hankel.shape, f'Shapes are off: {hankel_typed_matrix.shape=} != {materialized_hankel.shape=}.' + assert materialized_hankel.ndim == 2, 'Materialized matrix is strange.' + assert np.allclose(hankel_typed_matrix.materialize(), materialized_hankel), 'Hankel matrix is off.' + + # create the right matrix + right_matrix = np.arange(0, materialized_hankel.shape[1]*other_size) + right_matrix = right_matrix.reshape(materialized_hankel.shape[1], other_size) + + # check for right multiplication + naive_matmul = materialized_hankel @ right_matrix + hankel_typed_matmul = hankel_typed_matrix @ right_matrix + is_close = np.allclose(naive_matmul, hankel_typed_matmul) + print('Right multiplication valid?', is_close) + assert is_close, 'Right multiplication failed.' + + # create the left matrix + left_matrix = np.arange(0, materialized_hankel.shape[0] * other_size) + left_matrix = left_matrix.reshape(other_size, materialized_hankel.shape[0]) + + # check for left multiplication + naive_matmul = left_matrix @ materialized_hankel + hankel_typed_matmul = left_matrix @ hankel_typed_matrix + is_close = np.allclose(naive_matmul, hankel_typed_matmul) + print('Left multiplication valid?', is_close) + assert is_close, 'Left multiplication failed.' + + # check for transpose + assert np.allclose(hankel_typed_matrix.T.materialize(), materialized_hankel.T), 'Transpose does not work.' + + # check for right transposed multiplication + naive_matmul = materialized_hankel.T @ left_matrix.T + hankel_typed_matmul = hankel_typed_matrix.T @ left_matrix.T + is_close = np.allclose(naive_matmul, hankel_typed_matmul) + print('Transposed right multiplication valid?', is_close) + assert is_close, 'Transposed right multiplication failed.' + + # check for left transposed multiplication + naive_matmul = right_matrix.T @ materialized_hankel.T + hankel_typed_matmul = right_matrix.T @ hankel_typed_matrix.T + is_close = np.allclose(naive_matmul, hankel_typed_matmul) + print('Transposed left multiplication valid?', is_close) + assert is_close, 'Transposed left multiplication failed.' + print('################################################') + print() + + +def main(): + import timeit + # ------------------------- + # Example + # ------------------------- + + p = 70 # block rows of H + q = 60 # block columns of H + m = 2 # rows per block + s = 15 # rows of left multiplier A + + # make a signal of the correct shape + signal = np.arange(0, (p+q-1)*m) + signal = signal.reshape(p+q-1, m) + + # make the Hankel block matrix representation + hankel_direct = BlockHankelRepresentation(signal, p+q-1, p, q, fast_hankel=False) + hankel_fft = BlockHankelRepresentation(signal, p+q-1, p, q, fast_hankel=True) + hankel_naive = hankel_direct.materialize() + + # check whether the materialized hankel matrices are equal + print('Materialized Hankel matrices equal?', np.allclose(hankel_direct.materialize(), hankel_fft.materialize())) + + # make some checks using the functions + check_matrix_operations(hankel_fft, hankel_naive, s) + check_matrix_operations(hankel_direct, hankel_naive, s) + + # create right matrix + right_matrix = np.arange(0, q * s) + right_matrix = right_matrix.reshape(q, s) + + # make a print to start multiplication from the right + print() + print('Starting multiplication timing from the RIGHT-------') + + run_num = 100 + fft_time = timeit.timeit(lambda: BlockHankelRepresentation(signal, p+q-1, p, q, fast_hankel=True) @ right_matrix, number=run_num) / run_num * 1000 + print('FFT right multiplication took:', fft_time) + direct_time = timeit.timeit(lambda: BlockHankelRepresentation(signal, p+q-1, p, q, fast_hankel=False) @ right_matrix, number=run_num) / run_num * 1000 + print('Direct right multiplication took:', direct_time) + direct_time = timeit.timeit(lambda: BlockHankelRepresentation(signal, p+q-1, p, q, fast_hankel=False).materialize() @ right_matrix, number=run_num) / run_num * 1000 + print('Naive right multiplication took:', direct_time) + direct_time = timeit.timeit(lambda: hankel_naive @ right_matrix, number=run_num) / run_num * 1000 + print('Naive (already constructed, unfair) right multiplication took:', direct_time) + + # make a print to start multiplication from the left + print() + print('Starting multiplication timing from the LEFT-------') + + # create left matrix + left_matrix = np.arange(0, p * m * s) + left_matrix = left_matrix.reshape(s, p*m) + + run_num = 100 + fft_time = timeit.timeit(lambda: left_matrix @ BlockHankelRepresentation(signal, p+q-1, p, q, fast_hankel=True), number=run_num) / run_num * 1000 + print('FFT left multiplication took:', fft_time) + direct_time = timeit.timeit(lambda: left_matrix @ BlockHankelRepresentation(signal, p+q-1, p, q, fast_hankel=False), number=run_num) / run_num * 1000 + print('Direct left multiplication took:', direct_time) + direct_time = timeit.timeit(lambda: left_matrix @ BlockHankelRepresentation(signal, p+q-1, p, q, fast_hankel=False).materialize(), number=run_num) / run_num * 1000 + print('Naive left multiplication took:', direct_time) + direct_time = timeit.timeit(lambda: left_matrix @ hankel_naive,number=run_num) / run_num * 1000 + print('Naive (already constructed, unfair) left multiplication took:', direct_time) + + # make a print to start working on concatenated hankel matrices + print() + print('Starting working on concatenated matrices-------') + + naive_multilevel = np.concatenate((hankel_naive, hankel_naive), axis=1) + direct_multilevel = np.concatenate((hankel_direct, hankel_direct), axis=1) + fft_multilevel = np.concatenate((hankel_fft, hankel_fft), axis=1) + + # make some checks using the functions + check_matrix_operations(fft_multilevel, naive_multilevel, s) + check_matrix_operations(direct_multilevel, naive_multilevel, s) + + # create the right matrix + right_matrix = np.arange(0, naive_multilevel.shape[1] * s) + right_matrix = right_matrix.reshape(naive_multilevel.shape[1], s) + + # create the right matrix + left_matrix = np.arange(0, naive_multilevel.shape[0] * s) + left_matrix = left_matrix.reshape(s, naive_multilevel.shape[0]) + + run_num = 100 + fft_time = timeit.timeit( + lambda: np.concatenate((BlockHankelRepresentation(signal, p+q-1, p, q, fast_hankel=True), + BlockHankelRepresentation(signal, p+q-1, p, q, fast_hankel=True)), axis=1) @ right_matrix, + number=run_num) / run_num * 1000 + print('FFT right multiplication took:', fft_time) + direct_time = timeit.timeit( + lambda: np.concatenate((BlockHankelRepresentation(signal, p+q-1, p, q, fast_hankel=False), + BlockHankelRepresentation(signal, p+q-1, p, q, fast_hankel=False)), axis=1) @ right_matrix, + number=run_num) / run_num * 1000 + print('Direct right multiplication took:', direct_time) + direct_time = timeit.timeit( + lambda: np.concatenate((hankel_direct.materialize(), hankel_direct.materialize()), axis=1) @ right_matrix, + number=run_num) / run_num * 1000 + print('Naive right multiplication took:', direct_time) + direct_time = timeit.timeit( + lambda: np.concatenate((hankel_naive, hankel_naive), axis=1) @ right_matrix, + number=run_num) / run_num * 1000 + print('Naive right (already constructed, unfair) multiplication took:', direct_time) + + run_num = 100 + fft_time = timeit.timeit( + lambda: left_matrix @ np.concatenate((BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=True), + BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=True)), + axis=1), + number=run_num) / run_num * 1000 + print('FFT left multiplication took:', fft_time) + direct_time = timeit.timeit( + lambda: left_matrix @ np.concatenate((BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=False), + BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=False)), + axis=1), + number=run_num) / run_num * 1000 + print('Direct left multiplication took:', direct_time) + direct_time = timeit.timeit( + lambda: left_matrix @ np.concatenate((hankel_direct.materialize(), hankel_direct.materialize()), axis=1), + number=run_num) / run_num * 1000 + print('Naive left multiplication took:', direct_time) + direct_time = timeit.timeit( + lambda: left_matrix @ np.concatenate((hankel_naive, hankel_naive), axis=1), + number=run_num) / run_num * 1000 + print('Naive left (already constructed, unfair) multiplication took:', direct_time) + + # make a print to start working on product hankel matrices + print() + print('Starting working on product Hankel matrices-------') + + hankel_product_direct = hankel_direct @ hankel_direct.T + hankel_product_fft = hankel_fft @ hankel_fft.T + hankel_product_naive = hankel_naive @ hankel_naive.T + + check_matrix_operations(hankel_product_direct, hankel_product_naive) + check_matrix_operations(hankel_product_fft, hankel_product_naive) + + # create the right matrix + right_matrix = np.arange(0, hankel_product_naive.shape[1] * s) + right_matrix = right_matrix.reshape(hankel_product_naive.shape[1], s) + + # create the left matrix + left_matrix = np.arange(0, hankel_product_naive.shape[0] * s) + left_matrix = left_matrix.reshape(s, hankel_product_naive.shape[0]) + + run_num = 100 + fft_time = timeit.timeit(lambda: left_matrix @ (hankel_fft @ hankel_fft.T), number=run_num) / run_num * 1000 + print('FFT left multiplication took:', fft_time) + direct_time = timeit.timeit(lambda: left_matrix @ (hankel_direct @ hankel_direct.T), + number=run_num) / run_num * 1000 + print('Direct left multiplication took:', direct_time) + direct_time = timeit.timeit(lambda: left_matrix @ (hankel_direct.materialize() @ hankel_direct.materialize().T), + number=run_num) / run_num * 1000 + print('Naive left multiplication took:', direct_time) + direct_time = timeit.timeit(lambda:left_matrix @ (hankel_naive @ hankel_naive.T), number=run_num) / run_num * 1000 + print('Naive left (already constructed, unfair) multiplication took:', direct_time) + + run_num = 100 + fft_time = timeit.timeit(lambda: (hankel_fft @ hankel_fft.T) @ right_matrix, number=run_num) / run_num * 1000 + print('FFT right multiplication took:', fft_time) + direct_time = timeit.timeit(lambda: (hankel_direct @ hankel_direct.T) @ right_matrix, + number=run_num) / run_num * 1000 + print('Direct right multiplication took:', direct_time) + direct_time = timeit.timeit(lambda: (hankel_direct.materialize() @ hankel_direct.materialize().T) @ right_matrix, + number=run_num) / run_num * 1000 + print('Naive right multiplication took:', direct_time) + direct_time = timeit.timeit(lambda: (hankel_naive @ hankel_naive.T) @ right_matrix, number=run_num) / run_num * 1000 + print('Naive right (already constructed, unfair) multiplication took:', direct_time) + + + +if __name__ == '__main__': + main() diff --git a/tests/test_block_linalg.py b/tests/test_block_linalg.py new file mode 100644 index 0000000..1c56641 --- /dev/null +++ b/tests/test_block_linalg.py @@ -0,0 +1,142 @@ +import numpy as np +import pytest + +import changepoynt.utils.block_linalg as blg + + +def _setup_signal(): + p = 70 # block rows of H + q = 60 # block columns of H + m = 2 # rows per block + s = 15 # rows of left multiplier A + + # make a signal of the correct shape + signal = np.arange(0, (p + q - 1) * m) + signal = signal.reshape(p + q - 1, m) + + return signal, p, q, m, s + + + +def test_block_hankel_direct(): + # make the signal + signal, p, q, m, s = _setup_signal() + + # make the Hankel block matrix representation + hankel_direct = blg.BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=False) + hankel_naive = hankel_direct.materialize() + + # run the checks + blg.check_matrix_operations(hankel_direct, hankel_naive, s) + + +def test_block_hankel_fft(): + + # make the signal + signal, p, q, m, s = _setup_signal() + + # make the Hankel block matrix representation + hankel_direct = blg.BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=False) + hankel_naive = hankel_direct.materialize() + hankel_fft = blg.BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=True) + + # run the checks + blg.check_matrix_operations(hankel_fft, hankel_naive, s) + + +def test_wrong_signal_shape(): + # make the signal + signal, p, q, m, s = _setup_signal() + + # make the Hankel block matrix representation + with pytest.raises(ValueError): + hankel_direct = blg.BlockHankelRepresentation(signal[:, 0], p + q - 1, p, q, fast_hankel=False) + + +def test_multilevel_block_hankel(): + + # make the signal + signal, p, q, m, s = _setup_signal() + geng = np.random.default_rng(42) + signal2 = signal + geng.random(size=signal.shape) + + # make the Hankel block matrix representation + hankel_direct = blg.BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=False) + hankel_direct2 = blg.BlockHankelRepresentation(signal2, p + q - 1, p, q, fast_hankel=False) + hankel_naive = hankel_direct.materialize() + hankel_naive2 = hankel_direct2.materialize() + hankel_fft = blg.BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=True) + hankel_fft2 = blg.BlockHankelRepresentation(signal2, p + q - 1, p, q, fast_hankel=True) + + for axis in (0, 1): + + naive_multilevel = np.concatenate((hankel_naive, hankel_naive2), axis=axis) + direct_multilevel = np.concatenate((hankel_direct, hankel_direct2), axis=axis) + fft_multilevel = np.concatenate((hankel_fft, hankel_fft2), axis=axis) + + # make some checks using the functions + blg.check_matrix_operations(fft_multilevel, naive_multilevel, s) + blg.check_matrix_operations(direct_multilevel, naive_multilevel, s) + +def test_wrong_axis_multilevel_block_hankel(): + # make the signal + signal, p, q, m, s = _setup_signal() + geng = np.random.default_rng(42) + signal2 = signal + geng.random(size=signal.shape) + + # make the Hankel block matrix representation + hankel_direct = blg.BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=False) + hankel_direct2 = blg.BlockHankelRepresentation(signal2, p + q - 1, p, q, fast_hankel=False) + + # check whether we stop wrong axis + with pytest.raises(ValueError): + direct_multilevel = np.concatenate((hankel_direct, hankel_direct2), axis=2) + + +def test_block_hankel_product(): + + # make the signal + signal, p, q, m, s = _setup_signal() + geng = np.random.default_rng(42) + signal2 = signal + geng.random(size=signal.shape) + + # make the Hankel block matrix representation + hankel_direct = blg.BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=False) + hankel_direct2 = blg.BlockHankelRepresentation(signal2, p + q - 1, p, q, fast_hankel=False) + hankel_naive = hankel_direct.materialize() + hankel_naive2 = hankel_direct2.materialize() + hankel_fft = blg.BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=True) + hankel_fft2 = blg.BlockHankelRepresentation(signal2, p + q - 1, p, q, fast_hankel=True) + + # compute the products + naive_product = hankel_naive @ hankel_naive2.T + direct_product = hankel_direct @ hankel_direct2.T + fft_product = hankel_fft @ hankel_fft2.T + + # make some checks using the functions + blg.check_matrix_operations(fft_product, naive_product, s) + blg.check_matrix_operations(direct_product, naive_product, s) + + +def test_block_hankel_product_wrong_shape(): + + # make the signal + signal, p, q, m, s = _setup_signal() + geng = np.random.default_rng(42) + signal2 = signal + geng.random(size=signal.shape) + + # make the Hankel block matrix representation + hankel_direct = blg.BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=False) + hankel_direct2 = blg.BlockHankelRepresentation(signal2, p + q - 1, p, q, fast_hankel=False) + hankel_naive = hankel_direct.materialize() + hankel_naive2 = hankel_direct2.materialize() + hankel_fft = blg.BlockHankelRepresentation(signal, p + q - 1, p, q, fast_hankel=True) + hankel_fft2 = blg.BlockHankelRepresentation(signal2, p + q - 1, p, q, fast_hankel=True) + + # compute the products + with pytest.raises(ValueError): + naive_product = hankel_naive @ hankel_naive2 + with pytest.raises(ValueError): + direct_product = hankel_direct @ hankel_direct2 + with pytest.raises(ValueError): + fft_product = hankel_fft @ hankel_fft2 \ No newline at end of file diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 1f68712..3084a31 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -8,9 +8,13 @@ class TestLinearAlgebra: def setup_method(self): # set a random seed np.random.seed(3455) + rng = np.random.default_rng(3455) # create a random square matrix - self.A = np.random.rand(50, 50) + # self.A = np.arange(50*50).reshape(50, 50)+ np.random.rand(50, 50) + U = rng.normal(size=(50, 3)) + V = rng.normal(size=(50, 3)) + self.A = U @ np.diag(np.array([5, 4, 3])) @ V.T + rng.normal(size=(50, 50))*1e-3 # create random starting vector for the power method # and normalize it @@ -79,13 +83,23 @@ def test_rayleigh_ritz_svd(self): def test_randomized_svd(self): # set a threshold for the amount of eigenvectors we want to have and compute them - k = 3 - k_eigvals, k_eigvecs = lg.facebook_randomized_svd(self.A, 50) + k = 2 + k_eigvals, k_eigvecs = lg.facebook_randomized_svd(self.A, 20) # compare the values np.testing.assert_almost_equal(k_eigvals[:k], self.singvals[:k]) np.testing.assert_almost_equal(np.abs(k_eigvecs[:, :k]), np.abs(self.singvecs[:, :k])) + def test_own_randomized_svd(self): + + # set a threshold for the amount of eigenvectors we want to have and compute them + k = 2 + k_left_singvec, k_singvals, _ = lg.randomized_hankel_svd(self.A, k, 2, oversampling_p=18) + + # compare the values + np.testing.assert_almost_equal(k_singvals[:k], self.singvals[:k]) + np.testing.assert_almost_equal(np.abs(k_left_singvec[:, :k]), np.abs(self.singvecs[:, :k])) + def test_right_hankel_product(self): # make the multiplications res = self.hankel_matrix_fft @ self.other_matrix diff --git a/tests/test_messt.py b/tests/test_messt.py new file mode 100644 index 0000000..e29eba9 --- /dev/null +++ b/tests/test_messt.py @@ -0,0 +1,130 @@ +import numpy as np +import pytest + +from changepoynt.algorithms.messt import MESST +import changepoynt.algorithms.esst as cpesst + + +def _make_change_signal( + n_per_segment: int = 320, + mean_before: int = 48, + mean_after: int = 14, + noise: float = 4, + seed: int = 1234, +): + rng = np.random.default_rng(seed) + left = np.ones(n_per_segment) * mean_before + right = np.ones(n_per_segment) * mean_after + signal = np.concatenate([left, right]) + signal += noise * rng.standard_normal(signal.shape[0]) + return signal, n_per_segment + + +def _outside_region(score: np.ndarray, center: int, half_width: int, valid_start: int) -> np.ndarray: + left = score[valid_start:max(valid_start, center - half_width)] + right = score[min(center + half_width, score.shape[0]):] + if left.size and right.size: + return np.concatenate([left, right]) + if left.size: + return left + if right.size: + return right + raise AssertionError("Need some samples outside the change region for comparison.") + + +def test_messt_unknown_method_raises_value_error(): + with pytest.raises(ValueError): + MESST(window_length=40, method="does-not-exist") + + +def test_messt_smoketest(): + allowed_methods = MESST(window_length=40).methods + signal, _ = _make_change_signal() + signal = np.concatenate([signal[..., None], signal[..., None]], axis=1) + for method in allowed_methods: + transf = MESST(window_length=40, method=method) + transf.transform(signal) + + +def test_messt_rejects_fast_hankel_for_fbrsvd(): + with pytest.raises(ValueError): + MESST(window_length=40, method="fbrsvd", use_fast_hankel=True) + + +def test_messt_rejects_1d_input(): + signal, _ = _make_change_signal() + detector = MESST(window_length=40, method="rsvd") + with pytest.raises(AssertionError): + detector.transform(signal) + + +def test_messt_rejects_too_short_signal(): + detector = MESST(window_length=40, n_windows=20, lag=20, method="rsvd") + too_short = np.linspace(0.0, 1.0, 80) + with pytest.raises(AssertionError): + detector.transform(too_short) + + +def test_messt_score_is_zero_before_first_possible_output(): + signal, _ = _make_change_signal() + detector = MESST(window_length=40, n_windows=20, lag=20, method="rsvd") + np.random.seed(7) + score = detector.transform(signal[..., None]) + + expected_first_scored_idx = detector.window_length + np.testing.assert_allclose(score[:expected_first_scored_idx], 0.0) + + +def test_messt_detects_frequency_change_near_boundary(): + signal, change_idx = _make_change_signal() + detector = MESST(window_length=48, n_windows=24, lag=24, rank=2, method="rsvd") + + np.random.seed(11) + score = detector.transform(signal[..., None]) + + neighborhood = score[change_idx - 60: change_idx + 60] + outside = _outside_region(score, center=change_idx, half_width=120, valid_start=detector.window_length) + + assert neighborhood.size > 0 + assert outside.size > 0 + assert np.isfinite(score).all() + assert neighborhood.max() > np.percentile(outside, 97) + + +def test_messt_fast_hankel_tracks_reference_implementation(): + signal, _ = _make_change_signal() + slow = MESST(window_length=40, n_windows=20, lag=20, rank=2, method="rsvd", use_fast_hankel=False) + fast = MESST(window_length=40, n_windows=20, lag=20, rank=2, method="rsvd", use_fast_hankel=True) + + np.random.seed(31) + slow_score = slow.transform(signal[..., None]) + np.random.seed(31) + fast_score = fast.transform(signal[..., None]) + + valid_start = slow.window_length + corr = np.corrcoef(slow_score[valid_start:], fast_score[valid_start:])[0, 1] + assert np.isfinite(corr) + assert corr > 0.95 + + +def test_messt_tracks_esst(): + signal, _ = _make_change_signal() + ws = 30 + slow = MESST(window_length=ws, n_windows=ws, lag=ws//2, rank=5, method="rsvd", use_fast_hankel=False) + fast = MESST(window_length=ws, n_windows=ws, lag=ws//2, rank=5, method="rsvd", use_fast_hankel=True) + orig_esst = cpesst.ESST(window_length=ws, n_windows=ws, lag=ws//2, rank=5, method="rsvd", use_fast_hankel=False) + + np.random.seed(31) + slow_score = slow.transform(signal[..., None]) + np.random.seed(31) + fast_score = fast.transform(signal[..., None]) + np.random.seed(31) + esst_score = orig_esst.transform(signal) + + valid_start = slow.window_length + corr = np.corrcoef(slow_score[valid_start:], esst_score[valid_start:])[0, 1] + assert np.isfinite(corr) + assert corr > 0.95 + corr = np.corrcoef(fast_score[valid_start:], esst_score[valid_start:])[0, 1] + assert np.isfinite(corr) + assert corr > 0.95 \ No newline at end of file diff --git a/tests/test_msst.py b/tests/test_msst.py new file mode 100644 index 0000000..f2449c8 --- /dev/null +++ b/tests/test_msst.py @@ -0,0 +1,130 @@ +import numpy as np +import pytest + +from changepoynt.algorithms.msst import MSST +import changepoynt.algorithms.sst as cpsst + + +def _make_change_signal( + n_per_segment: int = 320, + mean_before: int = 48, + mean_after: int = 14, + noise: float = 4, + seed: int = 1234, +): + rng = np.random.default_rng(seed) + left = np.ones(n_per_segment) * mean_before + right = np.ones(n_per_segment) * mean_after + signal = np.concatenate([left, right]) + signal += noise * rng.standard_normal(signal.shape[0]) + return signal, n_per_segment + + +def _outside_region(score: np.ndarray, center: int, half_width: int, valid_start: int) -> np.ndarray: + left = score[valid_start:max(valid_start, center - half_width)] + right = score[min(center + half_width, score.shape[0]):] + if left.size and right.size: + return np.concatenate([left, right]) + if left.size: + return left + if right.size: + return right + raise AssertionError("Need some samples outside the change region for comparison.") + + +def test_msst_unknown_method_raises_value_error(): + with pytest.raises(ValueError): + MSST(window_length=40, method="does-not-exist") + + +def test_msst_smoketest(): + allowed_methods = MSST(window_length=40).methods + signal, _ = _make_change_signal() + signal = np.concatenate([signal[..., None], signal[..., None]], axis=1) + for method in allowed_methods: + transf = MSST(window_length=40, method=method) + transf.transform(signal) + + + +def test_msst_rejects_fast_hankel_for_fbrsvd(): + with pytest.raises(ValueError): + MSST(window_length=40, method="fbrsvd", use_fast_hankel=True) + + +def test_msst_rejects_1d_input(): + signal, _ = _make_change_signal() + detector = MSST(window_length=40, method="rsvd") + with pytest.raises(AssertionError): + detector.transform(signal) + + +def test_msst_rejects_too_short_signal(): + detector = MSST(window_length=40, n_windows=20, lag=20, method="rsvd") + too_short = np.linspace(0.0, 1.0, 80) + with pytest.raises(AssertionError): + detector.transform(too_short) + + +def test_msst_score_is_zero_before_first_possible_output(): + signal, _ = _make_change_signal() + detector = MSST(window_length=40, n_windows=20, lag=20, method="rsvd") + np.random.seed(7) + score = detector.transform(signal[..., None]) + + expected_first_scored_idx = detector.window_length + np.testing.assert_allclose(score[:expected_first_scored_idx], 0.0) + + +def test_msst_detects_frequency_change_near_boundary(): + signal, change_idx = _make_change_signal() + detector = MSST(window_length=48, n_windows=24, lag=24, rank=2, method="rsvd") + + np.random.seed(11) + score = detector.transform(signal[..., None]) + + neighborhood = score[change_idx - 60: change_idx + 60] + outside = _outside_region(score, center=change_idx, half_width=120, valid_start=detector.window_length) + + assert neighborhood.size > 0 + assert outside.size > 0 + assert np.isfinite(score).all() + assert neighborhood.max() > np.percentile(outside, 98) + + +def test_msst_fast_hankel_tracks_reference_implementation(): + signal, _ = _make_change_signal() + slow = MSST(window_length=40, n_windows=20, lag=20, rank=2, method="rsvd", use_fast_hankel=False) + fast = MSST(window_length=40, n_windows=20, lag=20, rank=2, method="rsvd", use_fast_hankel=True) + + np.random.seed(31) + slow_score = slow.transform(signal[..., None]) + np.random.seed(31) + fast_score = fast.transform(signal[..., None]) + + valid_start = slow.window_length + corr = np.corrcoef(slow_score[valid_start:], fast_score[valid_start:])[0, 1] + assert np.isfinite(corr) + assert corr > 0.95 + + +def test_msst_tracks_esst(): + signal, _ = _make_change_signal() + slow = MSST(window_length=40, n_windows=20, lag=20, rank=2, method="rsvd", use_fast_hankel=False) + fast = MSST(window_length=40, n_windows=20, lag=20, rank=2, method="rsvd", use_fast_hankel=True) + orig_esst = cpsst.SST(window_length=40, n_windows=20, lag=20, rank=2, method="rsvd", use_fast_hankel=False) + + np.random.seed(31) + slow_score = slow.transform(signal[..., None]) + np.random.seed(31) + fast_score = fast.transform(signal[..., None]) + np.random.seed(31) + esst_score = orig_esst.transform(signal) + + valid_start = slow.window_length + corr = np.corrcoef(slow_score[valid_start:], esst_score[valid_start:])[0, 1] + assert np.isfinite(corr) + assert corr > 0.95 + corr = np.corrcoef(fast_score[valid_start:], esst_score[valid_start:])[0, 1] + assert np.isfinite(corr) + assert corr > 0.95 \ No newline at end of file