From 344651c32db2b60ec42ec6fed0007933156c55ba Mon Sep 17 00:00:00 2001 From: vloison Date: Fri, 4 Jul 2025 15:37:36 +0200 Subject: [PATCH 01/41] homogeneize fadin and unhap init arguments --- fadin/solver.py | 143 +++++++++++++++++++------------ fadin/utils/compute_constants.py | 4 +- 2 files changed, 91 insertions(+), 56 deletions(-) diff --git a/fadin/solver.py b/fadin/solver.py index 8767cf5..9cdef46 100644 --- a/fadin/solver.py +++ b/fadin/solver.py @@ -88,6 +88,12 @@ class FaDIn(object): optim : `str` in ``{'RMSprop' | 'Adam' | 'GD'}``, default='RMSprop' The algorithms used to optimize the Hawkes processes parameters. + params_optim : dict, {'lr', ...}, default=dict() + Learning rate and parameters of the chosen optimization algorithm. + Will be passed as arguments to the `torch.optimizer` constructor chosen + via the `optim` parameter. + If 'lr' is not given, it is set to 1e-3. + step_size : `float`, `default=1e-3` Learning rate of the chosen optimization algorithm. @@ -100,9 +106,6 @@ class FaDIn(object): ztzG is approximated with Toeplitz matrix not taking into account edge effects. - log : `boolean`, `default=False` - Record the loss values during the optimization. - grad_kernel : `None` or `callable`, default=None If kernel in ``{'raised_cosine'| 'truncated_gaussian' | 'truncated_exponential'}`` the gradient function is implemented. @@ -113,21 +116,39 @@ class FaDIn(object): criterion is below it). If not reached the solver does 'max_iter' iterations. + random_state : `int`, `RandomState` instance or `None`, `default=None` + Set the torch seed to 'random_state'. + If set to `None`, torch seed will be set to 0. + Attributes ---------- - param_baseline : `tensor`, shape (n_dim) - Baseline parameter of the Hawkes process. + params_intens : `list` of `tensor` + List of parameters of the Hawkes process at the end of the fit. + The list contains the following tensors in this order: + - `baseline`: `tensor`, shape (n_dim,): Baseline parameter. + - `alpha`: `tensor`, shape (n_dim, n_dim): Weight parameter. + - `kernel`: `list` of `tensor`, shape (n_dim, n_dim): + Kernel parameters. The size of the list varies depending + the number of parameters. The shape of each tensor is + `(n_dim, n_dim)`. + + param_baseline : `tensor`, shape (max_iter, n_dim) + Baseline parameter of the Hawkes process for each fit iteration. - param_alpha : `tensor`, shape (n_dim, n_dim) - Weight parameter of the Hawkes process. + param_baseline_noise : `tensor`, shape (max_iter, n_dim) + Baseline parameter of the Hawkes process for each fit iteration. + + param_alpha : `tensor`, shape (max_iter, n_dim, n_dim) + Weight parameter of the Hawkes process for each fit iteration. param_kernel : `list` of `tensor` - list containing tensor array of kernels parameters. + list containing tensor array of kernels parameters for each fit + iteration. The size of the list varies depending the number of parameters. The shape of each tensor is `(n_dim, n_dim)`. v_loss : `tensor`, shape (n_iter) - If `log=True`, compute the loss accross iterations. + loss accross iterations. If no early stopping, `n_iter` is equal to `max_iter`. """ compute_gradient = staticmethod(compute_gradient_fadin) @@ -136,8 +157,7 @@ class FaDIn(object): def __init__(self, n_dim, kernel, init='random', optim_mask=None, kernel_length=1, delta=0.01, optim='RMSprop', params_optim=dict(), max_iter=2000, ztzG_approx=True, - log=False, grad_kernel=None, - tol=10e-5, random_state=None): + grad_kernel=None, tol=10e-5, random_state=None): # Discretization parameters self.delta = delta @@ -149,7 +169,6 @@ def __init__(self, n_dim, kernel, init='random', optim_mask=None, self.kernel = kernel self.solver = optim self.max_iter = max_iter - self.log = log self.tol = tol self.n_dim = n_dim self.kernel_model = DiscreteKernelFiniteSupport( @@ -347,18 +366,25 @@ class UNHaP(object): - 'kernel': `list` of tensors of shape (n_dim, n_dim): Initial kernel parameters. - baseline_mask : `tensor` of shape (n_dim,), or `None` - Tensor of same shape as the baseline vector, with values in (0, 1). - `baseline` coordinates where `baseline_mask` is equal to 0 - will stay constant equal to zero and not be optimized. - If set to `None`, all coordinates of baseline will be optimized. - - alpha_mask : `tensor` of shape (n_dim, n_dim), or `None` - Tensor of same shape as the `alpha` tensor, with values in (0, 1). - `alpha` coordinates and kernel parameters where `alpha_mask` = 0 - will not be optimized. - If set to `None`, all coordinates of alpha will be optimized, - and all kernel parameters will be optimized if optimize_kernel=`True`. + optim_mask: `dict` of `tensor` or `None`, default=`None`. + Dictionary containing the masks for the optimization of the parameters + of the Hawkes process. If set to `None`, all parameters are optimized. + The dictionary must contain the following keys: + - 'baseline': `tensor` of shape (n_dim,), or `None`. + Tensor of same shape as the baseline vector, with values in (0, 1). + `baseline` coordinates where then tensor is equal to 0 + will not be optimized. If set to `None`, all coordinates of + baseline will be optimized. + - 'baseline_noise': `tensor` of shape (n_dim,), or `None`. + Tensor of same shape as the noise baseline vector, with values in + (0, 1). `baseline_noise` coordinates where then tensor is equal + to 0 will not be optimized. If set to `None`, all coordinates of + baseline_noise will be optimized. + - 'alpha': `tensor` of shape (n_dim, n_dim), or `None`. + Tensor of same shape as the `alpha` tensor, with values in (0, 1). + `alpha` coordinates and kernel parameters where `alpha_mask` = 0 + will not be optimized. If set to `None`, all coordinates of alpha + and kernel parameters will be optimized. kernel_length : `float`, `default=1.` Length of kernels in the Hawkes process. @@ -369,68 +395,79 @@ class UNHaP(object): optim : `str` in ``{'RMSprop' | 'Adam' | 'GD'}``, default='RMSprop' The algorithm used to optimize the Hawkes processes parameters. - params_optim : dict, {'lr', ...} + params_optim : dict, {'lr', ...}, default=dict() Learning rate and parameters of the chosen optimization algorithm. + Will be passed as arguments to the `torch.optimizer` constructor chosen + via the `optim` parameter. + If 'lr' is not given, it is set to 1e-3. max_iter : `int`, `default=2000` Maximum number of iterations during fit. + batch_rho : `int` + Number of FaDIn iterations between latent variables rho updates. + ztzG_approx : `boolean`, `default=True` If ztzG_approx is false, compute the true ztzG precomputation constant that is the computational bottleneck of FaDIn. if ztzG_approx is true, ztzG is approximated with Toeplitz matrix not taking into account edge effects. - batch_rho : `int` - Number of FaDIn iterations between latent variables rho updates. - tol : `float`, `default=1e-5` The tolerance of the solver (iterations stop when the stopping criterion is below it). If not reached the solver does 'max_iter' iterations. density_hawkes : `str`, `default=linear` - Density of the marks of the Hawkes process, in + Density of the marks of the Hawkes process events, in ``{'linear' | 'uniform'}``. density_noise : `str`, `default=reverse_linear` - Density of the marks of the Hawkes process, in + Density of the marks of the spurious events, in ``{'reverse_linear' | 'uniform'}``. - moment_matching: boolean, `default=False` - If set to False, baseline, alpha and kernel parameters are randomly - chosen. If set to True, baseline, alpha and kernel parameters are - chosen using the smart init strategy. - The smart init strategy is only implemented - for truncated gaussian and raised_cosine kernels. - random_state : `int`, `RandomState` instance or `None`, `default=None` Set the torch seed to 'random_state'. If set to `None`, torch seed will be set to 0. Attributes ---------- - param_baseline : `tensor`, shape (n_dim) - Baseline parameter of the Hawkes process. - - param_baseline_noise : `tensor`, shape (n_dim) - Baseline parameter of the Hawkes process. - - param_alpha : `tensor`, shape (n_dim, n_dim) - Weight parameter of the Hawkes process. + params_intens : `list` of `tensor` + List of parameters of the Hawkes process at the end of the fit. + The list contains the following tensors in this order: + - `baseline`: `tensor`, shape (n_dim,): Baseline parameter. + - `baseline_noise`: `tensor`, shape (n_dim,): Baseline noise parameter. + - `alpha`: `tensor`, shape (n_dim, n_dim): Weight parameter. + - `kernel`: `list` of `tensor`, shape (n_dim, n_dim): + Kernel parameters. The size of the list varies depending + the number of parameters. The shape of each tensor is + `(n_dim, n_dim)`. + + param_baseline : `tensor`, shape (max_iter, n_dim) + Baseline parameter of the Hawkes process for each fit iteration. + + param_baseline_noise : `tensor`, shape (max_iter, n_dim) + Baseline parameter of the Hawkes process for each fit iteration. + + param_alpha : `tensor`, shape (max_iter, n_dim, n_dim) + Weight parameter of the Hawkes process for each fit iteration. param_kernel : `list` of `tensor` - list containing tensor array of kernels parameters. + list containing tensor array of kernels parameters for each fit + iteration. The size of the list varies depending the number of parameters. The shape of each tensor is `(n_dim, n_dim)`. + v_loss : `tensor`, shape (n_iter) + loss accross iterations. + If no early stopping, `n_iter` is equal to `max_iter`. + """ def __init__(self, n_dim, kernel, init='random', optim_mask=None, kernel_length=1, delta=0.01, optim='RMSprop', params_optim=dict(), max_iter=2000, batch_rho=100, ztzG_approx=True, tol=10e-5, density_hawkes='linear', - density_noise='uniform', moment_matching=False, - random_state=None): + density_noise='uniform', random_state=None): # Set discretization parameters self.delta = delta @@ -447,7 +484,6 @@ def __init__(self, n_dim, kernel, init='random', optim_mask=None, self.density_hawkes = density_hawkes self.density_noise = density_noise - self.moment_matching = moment_matching # Set model parameters self.n_dim = n_dim @@ -572,11 +608,10 @@ def fit(self, events, end_time): # add rho parameter at the end of params self.params_mixture = [self.rho] - if self.moment_matching: - # Smart initialization of solver parameters - self.params_intens = init_hawkes_params_unhap( - self, self.init, marked_events, n_ground_events, end_time - ) + # Smart initialization of solver parameters + self.params_intens = init_hawkes_params_unhap( + self, self.init, marked_events, n_ground_events, end_time + ) self.opt_intens, self.opt_mixture = optimizer_unhap( [self.params_intens, self.params_mixture], diff --git a/fadin/utils/compute_constants.py b/fadin/utils/compute_constants.py index a3b024f..8aaa2aa 100644 --- a/fadin/utils/compute_constants.py +++ b/fadin/utils/compute_constants.py @@ -166,7 +166,7 @@ def compute_marked_quantities( marks_grid_hawkes = events_grid else: raise NotImplementedError( - "this density is not implemented \ + f"events mark density {density_hawkes} is not implemented \ must be in linear | uniform " ) @@ -178,7 +178,7 @@ def compute_marked_quantities( marks_grid_noise = events_grid else: raise NotImplementedError( - "this density is not implemented \ + f"noise mark density {density_noise} is not implemented \ must be in reverse_linear | uniform " ) From 6eac89d57d0acc42ee1156b49a7f72660578fa47 Mon Sep 17 00:00:00 2001 From: vloison Date: Fri, 4 Jul 2025 16:17:53 +0200 Subject: [PATCH 02/41] fix compute_gradient --- fadin/loss_and_gradient.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/fadin/loss_and_gradient.py b/fadin/loss_and_gradient.py index 9620108..2fa487d 100644 --- a/fadin/loss_and_gradient.py +++ b/fadin/loss_and_gradient.py @@ -35,14 +35,17 @@ def compute_gradient_fadin(solver, events_grid, discretization, discretization ) - if solver.log: - solver.v_loss[i] = \ - discrete_l2_loss_precomputation(solver.zG, solver.zN, solver.ztzG, - solver.params_intens[0], - solver.params_intens[1], - kernel, n_events, - solver.delta, - end_time).detach() + solver.v_loss[i] = \ + discrete_l2_loss_precomputation( + solver.zG, + solver.zN, + solver.ztzG, + solver.params_intens[0], + solver.params_intens[1], + kernel, n_events, + solver.delta, + end_time + ).detach() # Update baseline gradient solver.params_intens[0].grad = get_grad_baseline( solver.zG, From 76f811aa68b472a9d13c1f6c45e20896801fec76 Mon Sep 17 00:00:00 2001 From: vloison Date: Fri, 4 Jul 2025 17:38:47 +0200 Subject: [PATCH 03/41] fix unhap moment matching --- examples/plot_unhap.py | 2 +- fadin/init.py | 13 +++++++------ fadin/solver.py | 27 +++++++++------------------ 3 files changed, 17 insertions(+), 25 deletions(-) diff --git a/examples/plot_unhap.py b/examples/plot_unhap.py index d8bdf60..0295bda 100644 --- a/examples/plot_unhap.py +++ b/examples/plot_unhap.py @@ -119,6 +119,7 @@ def simulate_data(baseline, baseline_noise, alpha, end_time, seed=0): n_dim=1, kernel="truncated_gaussian", kernel_length=1.0, + init='moment_matching_mean', delta=delta, optim="RMSprop", params_optim={"lr": 1e-3}, @@ -126,7 +127,6 @@ def simulate_data(baseline, baseline_noise, alpha, end_time, seed=0): batch_rho=batch_rho, density_hawkes="linear", density_noise="uniform", - moment_matching=True, ) solver.fit(ev, end_time) diff --git a/fadin/init.py b/fadin/init.py index 271e1d6..7d96bc0 100644 --- a/fadin/init.py +++ b/fadin/init.py @@ -220,8 +220,8 @@ def momentmatching_fadin(solver, events, n_ground_events, end_time, return baseline, alpha, kernel_params_init -def init_hawkes_params_unhap(solver, init, events, n_ground_events, - end_time): +def init_hawkes_params_unhap(solver, init, events, events_grid, + n_ground_events, end_time): """ Computes the initial Hawkes parameters for the UNHaP solver. @@ -240,6 +240,7 @@ def init_hawkes_params_unhap(solver, init, events, n_ground_events, 'alpha' and 'kernel'. events: list of array of size number of timestamps, list size is self.n_dim. + events_grid: TODO n_ground_events : torch.tensor Number of ground events for each dimension end_time: float @@ -268,6 +269,7 @@ def init_hawkes_params_unhap(solver, init, events, n_ground_events, baseline, baseline_noise, alpha, kernel_params_init = momentmatching_unhap( solver, events, + events_grid, n_ground_events, end_time, mm_mode @@ -312,10 +314,9 @@ def momentmatching_unhap(solver, events, events_grid, n_ground_events, alpha = torch.ones(solver.n_dim, solver.n_dim) / (solver.n_dim + 2) # Kernel parameters init - if solver.optimize_kernel: - kernel_params_init = momentmatching_kernel( - solver, events, events_grid, n_ground_events, plot_delta, mode - ) + kernel_params_init = momentmatching_kernel( + solver, events, n_ground_events, plot_delta, mode + ) return baseline, baseline_noise, alpha, kernel_params_init diff --git a/fadin/solver.py b/fadin/solver.py index 9cdef46..997de4c 100644 --- a/fadin/solver.py +++ b/fadin/solver.py @@ -481,12 +481,19 @@ def __init__(self, n_dim, kernel, init='random', optim_mask=None, self.max_iter = max_iter self.tol = tol self.batch_rho = batch_rho + self.n_dim = n_dim + self.kernel_model = DiscreteKernelFiniteSupport( + self.delta, + self.n_dim, + kernel, + self.kernel_length, + 0 + ) self.density_hawkes = density_hawkes self.density_noise = density_noise # Set model parameters - self.n_dim = n_dim self.baseline = torch.rand(self.n_dim) @@ -529,22 +536,6 @@ def __init__(self, n_dim, kernel, init='random', optim_mask=None, f"Got {init}." ) self.init = init - self.params_intens = init_hawkes_params_unhap( - self, - self.init, - kernel, - self.kernel_length, - self.n_dim - ) - - self.kernel_model = DiscreteKernelFiniteSupport( - self.delta, - self.n_dim, - kernel, - self.kernel_length, - 0 - ) - self.kernel = kernel # If the learning rate is not given, fix it to 1e-3 if 'lr' not in params_optim.keys(): @@ -610,7 +601,7 @@ def fit(self, events, end_time): # Smart initialization of solver parameters self.params_intens = init_hawkes_params_unhap( - self, self.init, marked_events, n_ground_events, end_time + self, self.init, marked_events, events_grid, n_ground_events, end_time ) self.opt_intens, self.opt_mixture = optimizer_unhap( From 979626bf8949c1dadaf31119e6f610d03d715f91 Mon Sep 17 00:00:00 2001 From: vloison Date: Mon, 7 Jul 2025 10:52:15 +0200 Subject: [PATCH 04/41] update readme --- README.md | 28 ++++++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 5490ea4..91cc771 100644 --- a/README.md +++ b/README.md @@ -1,17 +1,27 @@ # FaDIn: Fast Discretized Inference For Hawkes Processes with General Parametric Kernels +# UNHaP: Unmixing Noise from Hawkes Processes ![build](https://img.shields.io/github/actions/workflow/status/GuillaumeStaermanML/FaDIn/unit_tests.yml?event=push&style=for-the-badge) ![python version](https://img.shields.io/badge/python-3.7_|_3.8_|_3.9_|_3.10_|_3.11-blue?style=for-the-badge) ![license](https://img.shields.io/github/license/GuillaumeStaermanML/FaDIn?style=for-the-badge) ![code style](https://img.shields.io/badge/code_style-black-black?style=for-the-badge) -This Package implements FaDIn. +This package implements FaDIn and UNHaP. FaDIn and UNHaP are solvers for inference of Hawkes Processes with finite-support kernels on simulated or real data, with the following features: +- Computation time is low compared to other methods. +- Compatible in an univariate setting as well as a multivariate setting. +- Classical kernels (exponential truncated gaussian, raised cosine) are implemented. The user can add their own kernel for inference. +- Masking: if only a few Hawkes Parameters need to be inferred, the user can mask the other parameters. +- Smart initialization of parameters before optimization: the user can choose between `random` (purely random), `moment_matching_max` (moment matching with maximum mode) and `moment_matching_mean` (moment matching with mean mode). + + +[FaDIn](https://proceedings.mlr.press/v202/staerman23a/staerman23a.pdf) does classical Hawkes inference with gradient descent. +[UNHaP](https://raw.githubusercontent.com/mlresearch/v258/main/assets/loison25a/loison25a.pdf) does Hawkes inference where the Hawkes Process is marked and mixed with a noisy Poisson process. + ## Installation **To install this package, make sure you have an up-to-date version of** `pip`. - ### From PyPI (coming soon) In a dedicated Python env, run: @@ -41,6 +51,12 @@ pip install -e ".[dev]" pre-commit install ``` +## Short examples +A few illustrative examples are provided in the `examples` folder of this repository: +- `plot_univariate_fadin`: simulate an univariate unmarked Hawkes process, infer Hawkes Process parameters using FaDIn, and plot inferred kernel. +- `plot_multivariate_fadin`: same as `plot_univariate_fadin` but in the multivariate case. +- `plot_unhap`: simulate an univariate marked Hawkes process and a marked Poisson process, infer Hawkes Process parameters using UNHaP, ald plot inferred kernels. + ## Citing this work If this package was useful to you, please cite it in your work: @@ -54,4 +70,12 @@ If this package was useful to you, please cite it in your work: year={2023}, organization={PMLR} } + +@improceedings{loison2025unhap, +title={UNHaP: Unmixing Noise from Hawkes Process}, +author={Loison, Virginie and Staerman, Guillaume and Moreau, Thomas}, +booktitle={International Conference on Artificial Intelligence and Statistics}, +pages={1342--1350}, +year={2025} +} ``` From 7df13a2b77d539ad22c8de421886113195995008 Mon Sep 17 00:00:00 2001 From: vloison Date: Mon, 7 Jul 2025 11:04:08 +0200 Subject: [PATCH 05/41] small readme update --- README.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 91cc771..c5e9925 100644 --- a/README.md +++ b/README.md @@ -21,7 +21,9 @@ This package implements FaDIn and UNHaP. FaDIn and UNHaP are solvers for inferen ## Installation **To install this package, make sure you have an up-to-date version of** `pip`. - +```bash +python3 -m pip install --upgrade pip +``` ### From PyPI (coming soon) In a dedicated Python env, run: From c0eda903143f73ebb79910f2ce01c5eb0fdce8ed Mon Sep 17 00:00:00 2001 From: vloison Date: Tue, 8 Jul 2025 15:34:44 +0200 Subject: [PATCH 06/41] add moment matchig test --- fadin/tests/test_init.py | 130 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 130 insertions(+) create mode 100644 fadin/tests/test_init.py diff --git a/fadin/tests/test_init.py b/fadin/tests/test_init.py new file mode 100644 index 0000000..58d2a88 --- /dev/null +++ b/fadin/tests/test_init.py @@ -0,0 +1,130 @@ +import pytest +import torch +import numpy as np +from fadin.init import ( + momentmatching_kernel, + init_hawkes_params_fadin, + momentmatching_fadin, + init_hawkes_params_unhap, + momentmatching_unhap, + random_params_fadin, + random_params_unhap, +) +# from fadin.kernels import truncated_gaussian +from fadin.utils.functions import identity, linear_zero_one +from fadin.utils.functions import reverse_linear_zero_one, truncated_gaussian + +from fadin.utils.utils_simu import simu_marked_hawkes_cluster, custom_density +from fadin.utils.utils_simu import simu_multi_poisson +from fadin.solver import FaDIn, UNHaP + +baseline = np.array([0.3]) +baseline_noise = np.array([0.05]) +alpha = np.array([[1.45]]) +mu = np.array([[0.4]]) +sigma = np.array([[0.1]]) + +delta = 0.01 +end_time = 1000 +seed = 0 +max_iter = 20 +batch_rho = 200 + + +# Create the simulating function +def simulate_marked_data(baseline, baseline_noise, alpha, end_time, seed=0): + n_dim = len(baseline) + + marks_kernel = identity + marks_density = linear_zero_one + time_kernel = truncated_gaussian + + params_marks_density = dict() + # params_marks_density = dict(scale=1) + params_marks_kernel = dict(slope=1.2) + params_time_kernel = dict(mu=mu, sigma=sigma) + + marked_events, _ = simu_marked_hawkes_cluster( + end_time, + baseline, + alpha, + time_kernel, + marks_kernel, + marks_density, + params_marks_kernel=params_marks_kernel, + params_marks_density=params_marks_density, + time_kernel_length=None, + marks_kernel_length=None, + params_time_kernel=params_time_kernel, + random_state=seed, + ) + + noisy_events_ = simu_multi_poisson(end_time, [baseline_noise]) + + random_marks = [np.random.rand(noisy_events_[i].shape[0]) for i in range(n_dim)] + noisy_marks = [ + custom_density( + reverse_linear_zero_one, + dict(), + size=noisy_events_[i].shape[0], + kernel_length=1.0, + ) + for i in range(n_dim) + ] + noisy_events = [ + np.concatenate( + (noisy_events_[i].reshape(-1, 1), random_marks[i].reshape(-1, 1)), axis=1 + ) + for i in range(n_dim) + ] + + events = [ + np.concatenate((noisy_events[i], marked_events[i]), axis=0) + for i in range(n_dim) + ] + + events_cat = [events[i][events[i][:, 0].argsort()] for i in range(n_dim)] + + labels = [ + np.zeros(marked_events[i].shape[0] + noisy_events_[i].shape[0]) + for i in range(n_dim) + ] + labels[0][-marked_events[0].shape[0]:] = 1.0 + true_rho = [labels[i][events[i][:, 0].argsort()] for i in range(n_dim)] + + return events_cat, noisy_marks, true_rho + + +ev, noisy_marks, true_rho = simulate_marked_data( + baseline, baseline_noise.item(), alpha, end_time, seed=0 +) + + +def test_unhap_init(): + unhap_mmmax = UNHaP( + n_dim=1, + kernel="truncated_gaussian", + init="moment_matching_max", + max_iter=max_iter + ) + unhap_mmmax.fit(ev, end_time) + + unhap_mmtimean = UNHaP( + n_dim=1, + kernel="truncated_gaussian", + init="moment_matching_mean", + max_iter=max_iter + ) + unhap_mmtimean.fit(ev, end_time) + + unhap_random = UNHaP( + n_dim=1, + kernel="truncated_gaussian", + init="random", + max_iter=max_iter + ) + unhap_random.fit(ev, end_time) + assert unhap_mmmax is not None, "UNHaP moment matching max failed" + assert unhap_mmtimean is not None, "UNHaP moment matching mean failed" + assert unhap_random is not None, "UNHaP random initialization failed" + return None From 40692171eedd96c375a0e7881413f9acf1c8135c Mon Sep 17 00:00:00 2001 From: vloison Date: Tue, 8 Jul 2025 15:34:59 +0200 Subject: [PATCH 07/41] minor changes --- examples/plot_unhap.py | 6 +++--- fadin/tests/test_mask.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/plot_unhap.py b/examples/plot_unhap.py index 0295bda..69a83ee 100644 --- a/examples/plot_unhap.py +++ b/examples/plot_unhap.py @@ -45,7 +45,7 @@ # %% Create the simulating function -def simulate_data(baseline, baseline_noise, alpha, end_time, seed=0): +def simulate_marked_data(baseline, baseline_noise, alpha, end_time, seed=0): n_dim = len(baseline) marks_kernel = identity @@ -102,7 +102,7 @@ def simulate_data(baseline, baseline_noise, alpha, end_time, seed=0): np.zeros(marked_events[i].shape[0] + noisy_events_[i].shape[0]) for i in range(n_dim) ] - labels[0][-marked_events[0].shape[0] :] = 1.0 + labels[0][-marked_events[0].shape[0]:] = 1.0 true_rho = [labels[i][events[i][:, 0].argsort()] for i in range(n_dim)] # put the mark to one to test the impact of the marks # events_cat[0][:, 1] = 1. @@ -110,7 +110,7 @@ def simulate_data(baseline, baseline_noise, alpha, end_time, seed=0): return events_cat, noisy_marks, true_rho -ev, noisy_marks, true_rho = simulate_data( +ev, noisy_marks, true_rho = simulate_marked_data( baseline, baseline_noise.item(), alpha, end_time, seed=0 ) # %% Apply UNHAP diff --git a/fadin/tests/test_mask.py b/fadin/tests/test_mask.py index 5dfb23a..86a0346 100644 --- a/fadin/tests/test_mask.py +++ b/fadin/tests/test_mask.py @@ -127,6 +127,6 @@ def test_rc_mask(): ztzG_approx=ztzG_approx, random_state=simu_random_state ) - assert torch.allclose(rc_bl, torch.Tensor([0., 0.])) - assert torch.allclose(rc_alpha * torch.Tensor([[1., 1.], [0., 1.]]), - torch.zeros(2, 2)) + assert torch.allclose(rc_bl, torch.Tensor([0., 0.])) + assert torch.allclose(rc_alpha * torch.Tensor([[1., 1.], [0., 1.]]), + torch.zeros(2, 2)) From 0b8a32caedf114feb5c60408e75d27aaa679e711 Mon Sep 17 00:00:00 2001 From: vloison Date: Tue, 8 Jul 2025 15:38:41 +0200 Subject: [PATCH 08/41] fix linter --- fadin/tests/test_init.py | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/fadin/tests/test_init.py b/fadin/tests/test_init.py index 58d2a88..8972e2e 100644 --- a/fadin/tests/test_init.py +++ b/fadin/tests/test_init.py @@ -1,22 +1,11 @@ -import pytest -import torch import numpy as np -from fadin.init import ( - momentmatching_kernel, - init_hawkes_params_fadin, - momentmatching_fadin, - init_hawkes_params_unhap, - momentmatching_unhap, - random_params_fadin, - random_params_unhap, -) # from fadin.kernels import truncated_gaussian from fadin.utils.functions import identity, linear_zero_one from fadin.utils.functions import reverse_linear_zero_one, truncated_gaussian from fadin.utils.utils_simu import simu_marked_hawkes_cluster, custom_density from fadin.utils.utils_simu import simu_multi_poisson -from fadin.solver import FaDIn, UNHaP +from fadin.solver import UNHaP baseline = np.array([0.3]) baseline_noise = np.array([0.05]) From c93e84a97747e4fa65ccbb07c60d1b326293d8a6 Mon Sep 17 00:00:00 2001 From: vloison Date: Tue, 8 Jul 2025 17:57:54 +0200 Subject: [PATCH 09/41] add unhap ecg experiments --- experiments/unhap/ecg_gait/run_ecg.py | 449 ++++++++++++++ experiments/unhap/ecg_gait/utils_ecg.py | 785 ++++++++++++++++++++++++ 2 files changed, 1234 insertions(+) create mode 100644 experiments/unhap/ecg_gait/run_ecg.py create mode 100644 experiments/unhap/ecg_gait/utils_ecg.py diff --git a/experiments/unhap/ecg_gait/run_ecg.py b/experiments/unhap/ecg_gait/run_ecg.py new file mode 100644 index 0000000..8f0b7eb --- /dev/null +++ b/experiments/unhap/ecg_gait/run_ecg.py @@ -0,0 +1,449 @@ +"""Apply Hawkes processes solvers to ECG data. +Done following vitaldb example notebook: +https://github.com/vitaldb/examples/blob/master/vitaldb_python_library.ipynb . + +Iteratively runs a set of chosen solvers on 20 5-minute ECG intervals. +Computes the mean absolute error and mean relative absolute error between +the heart rate estimated by each solver and the ground truth heart rate. + +""" +import time +from pathlib import Path + + +import numpy as np +import pandas as pd +from joblib import Memory + +from pyhrv.hrv import hrv +import neurokit2 as nk + +# %% Imports +from fadin.solver import FaDIn, UNHaP + +from utils_ecg import dataload, RESAMP_FREQ, ECG_CH_NAME +from utils_ecg import peak_events, plot_signal_rho +from utils_ecg import process_cdl_actis, fadin_fit, cdl + +memory = Memory(location='__cache__', verbose=10) + +# Segment time intervals for ECG data from vitaldb +SLOTS = { + 0: ['https://api.vitaldb.net/0010.vital', 300, 305], + 2: ['https://api.vitaldb.net/0001.vital', 35, 40], + 1: ['https://api.vitaldb.net/0001.vital', 20, 25], + 3: ['https://api.vitaldb.net/0001.vital', 55, 60], + 4: ['https://api.vitaldb.net/0002.vital', 55, 60], + 5: ['https://api.vitaldb.net/0002.vital', 120, 125], # 125, 130 + 6: ['https://api.vitaldb.net/0002.vital', 150, 155], + 7: ['https://api.vitaldb.net/0003.vital', 45, 50], + 8: ['https://api.vitaldb.net/0004.vital', 55, 60], + 9: ['https://api.vitaldb.net/0004.vital', 105, 110], + 10: ['https://api.vitaldb.net/0004.vital', 120, 125], + 12: ['https://api.vitaldb.net/0006.vital', 10, 15], + 13: ['https://api.vitaldb.net/0006.vital', 25, 30], + 14: ['https://api.vitaldb.net/0007.vital', 45, 50], + 15: ['https://api.vitaldb.net/0007.vital', 50, 55], + 16: ['https://api.vitaldb.net/0007.vital', 115, 120], + 17: ['https://api.vitaldb.net/0010.vital', 120, 125], + 18: ['https://api.vitaldb.net/0010.vital', 175, 180], + 19: ['https://api.vitaldb.net/0010.vital', 230, 235], +} + +# FaDIn solver hyperparameters +KERNEL_LENGTH = 1.5 +DELTA = 0.01 +LR = 1e-3 +L = int(1 / DELTA) +KERNEL_TYPE = 'truncated_gaussian' +N_ITERATIONS_FADIN = 50_000 +RANDOM_STATE = 40 # Only used if model=FaDIn or smart_init=False + +# Dictionary learning hyperparameters +N_ATOMS = 1 +N_TIMES_ATOMS = int(RESAMP_FREQ) # int(0.7 * RESAMP_FREQ) +N_ITERATIONS = 40 +SORT_ATOMS = True # Return the atoms in decreasing signal explanability +WINDOW = True # Whether the atoms are forced to be = 0 on their edges +RANDOM_SEED = 0 +REG = 2.5 # reg parameter of BatchCDL function 1.5, 3 + +# Module choice +MULTI_START = False +CDL_OR_PEAKS = 'cdl' # 'cdl' or 'peaks' +RUN_BASELINES = False +RUN_SOLVER = True + +# Change following line to choose which solver models will be used. +MODEL_LIST = [UNHaP, 'StocUNHaP'] +BATCH_RHO = 1000 +INIT = 'moment_matching_max' # 'moment_matching_mean', 'moment_matching_max', 'random' + +# CDL activations processing parameters: +MARKED_EVENTS = True +ACTI_THRESH = 0. +SPARSIFY_ACTIS = False + +# Path for saving figure. If set to None, figures are not saved. +SAVEFIG = None # 'ECG_plots/august_2024/' +RESPATH = Path('results') / 'ECG' / CDL_OR_PEAKS +RESPATH.mkdir(parents=True, exist_ok=True) + +if MULTI_START: + # indexes = np.squeeze( + # np.repeat(list(SLOTS), 4).reshape((1, -1), order='F') + # ) + # print('indexes', indexes) + df_runs = pd.DataFrame( + columns=['slot', 'init_mode', 'log_likelihood', 'abs_error', + 'rel_abs_error', 'bl_noise', 'bl_hawkes', 'alpha_hawkes', + 'mean', 'std'], + ) + + +# Functions + +# Initialize Dataframes to save errors +ABS_ERROR = pd.DataFrame( + data=100., + index=SLOTS.keys(), + columns=[str(m) for m in MODEL_LIST] + ['pyHRV'] + ['Neurokit'] +) +REL_ABS_ERROR = pd.DataFrame( + data=1., + index=SLOTS.keys(), + columns=[str(m) for m in MODEL_LIST] + ['pyHRV'] + ['Neurokit'] +) +RUNTIME = pd.DataFrame( + data=1., + index=SLOTS.keys(), + columns=[str(m) for m in MODEL_LIST] + ['pyHRV'] + ['Neurokit'] +) +############################################################################## +# PYHRV AND NEUROKIT +############################################################################## +if RUN_BASELINES: + for KEY in SLOTS.keys(): + FILEURL = SLOTS[KEY][0] + TIME_START = SLOTS[KEY][1] # minutes + TIME_STOP = SLOTS[KEY][2] # minutes + print('File', FILEURL) + print(TIME_START, '-', TIME_STOP, 'time interval') + + # Load data + T, clean_df, cleandf_hr = dataload( + FILEURL, + time_start=TIME_START, + time_stop=TIME_STOP, + verbose=False, + plotfig=False + ) + print('Ground truth heart rate:', np.mean(cleandf_hr)) + # pyHRV + pyhrv_start = time.time() + hrv_stats = hrv( + signal=clean_df[ECG_CH_NAME].to_numpy(), sampling_rate=RESAMP_FREQ + ) + mean_hr_pyhrv = hrv_stats['hr_mean'] + print('pyHRV Mean Heart rate:', mean_hr_pyhrv) + + ae_pyhrv = np.abs(mean_hr_pyhrv - np.mean(cleandf_hr)) + rae_pyhrv = ae_pyhrv / np.mean(cleandf_hr) + ABS_ERROR.loc[KEY, 'pyHRV'] = ae_pyhrv + REL_ABS_ERROR.loc[KEY, 'pyHRV'] = rae_pyhrv + pyhrv_time = time.time() - pyhrv_start + print('pyHRV time:', pyhrv_time) + RUNTIME.loc[KEY, 'pyHRV'] = pyhrv_time + # Neurokit + nk_start = time.time() + signals, info = nk.ecg_process( + clean_df[ECG_CH_NAME].to_numpy(dtype=np.float32), + sampling_rate=RESAMP_FREQ, + method='neurokit' + ) + mean_hr_nk = signals['ECG_Rate'].mean() + print('Neurokit Mean Heart rate:', mean_hr_nk) + + ae_nk = np.abs(mean_hr_nk - np.mean(cleandf_hr)) + rae_nk = ae_nk / np.mean(cleandf_hr) + ABS_ERROR.loc[KEY, 'Neurokit'] = ae_nk + REL_ABS_ERROR.loc[KEY, 'Neurokit'] = rae_nk + nk_time = time.time() - nk_start + print('Neurokit time:', nk_time) + RUNTIME.loc[KEY, 'Neurokit'] = nk_time + + ABS_ERROR.to_csv(f'{RESPATH}abs_error_pyHRV_neurokit.csv') + REL_ABS_ERROR.to_csv(f'{RESPATH}rel_abs_error_pyHRV_neurokit.csv') + RUNTIME.to_csv(f'{RESPATH}runtime_pyHRV_neurokit.csv') + + print('pyHRV mean absolute error:', ABS_ERROR['pyHRV'].mean()) + print('pyHRV absolute error std:', ABS_ERROR['pyHRV'].std()) + print('pyHRV mean relative absolute error:', REL_ABS_ERROR['pyHRV'].mean()) + print('pyHRV relative absolute error std:', REL_ABS_ERROR['pyHRV'].std()) + print('Neurokit mean absolute error:', ABS_ERROR['Neurokit'].mean()) + print('Neurokit absolute error std:', ABS_ERROR['Neurokit'].std()) + print( + 'Neurokit mean relative absolute error:', + REL_ABS_ERROR['Neurokit'].mean() + ) + print( + 'Neurokit relative absolute error std:', + REL_ABS_ERROR['Neurokit'].std() + ) + print('PyHRV mean runtime:', RUNTIME['pyHRV'].mean()) + print('PyHRV std runtime:', RUNTIME['pyHRV'].std()) + print('Neurokit mean runtime:', RUNTIME['Neurokit'].mean()) + print('Neurokit std runtime:', RUNTIME['Neurokit'].std()) +############################################################################## +# EVENT DETECTION + UNHAP +############################################################################## +for KEY in SLOTS.keys(): + FILEURL = SLOTS[KEY][0] + TIME_START = SLOTS[KEY][1] # minutes + TIME_STOP = SLOTS[KEY][2] # minutes + print('File', FILEURL) + print(TIME_START, '-', TIME_STOP, 'time interval') + + # Load data + T, clean_df, cleandf_hr = dataload( + FILEURL, + time_start=TIME_START, + time_stop=TIME_STOP, + verbose=False, + plotfig=False + ) + model_start = time.time() + # Naive event detection + # events_t = naive_events( + # clean_df, time_start=TIME_START, thresh=0.5, plotfig=False + # ) + + if CDL_OR_PEAKS == 'peaks': + # Event detection using peaks + acti_events = peak_events(clean_df, RESAMP_FREQ, marked=MARKED_EVENTS) + + if CDL_OR_PEAKS == 'cdl': + # Event detection using convolutional dictionary learning + cdl_actis = cdl(clean_df, REG, KEY, plotfig=False, save_fig=SAVEFIG) + + # Process activations + acti_events = process_cdl_actis( + cdl_actis, + freq=RESAMP_FREQ, + reg=REG, + threshold=ACTI_THRESH, + sparsify=SPARSIFY_ACTIS, + marked=MARKED_EVENTS, + save_fig=SAVEFIG + ) + + # Run solver on events + for model in MODEL_LIST: + if not RUN_SOLVER: + break + if MULTI_START: + df_runs.to_csv(f'results/ECG/{MODEL_LIST}df_runs_ecg.csv') + assert model != FaDIn or not MARKED_EVENTS, ( + "FaDIn cannot be fed marked events. Change MARKED_EVENTS to False." + ) + print(str(model)) + if model in [UNHaP, 'StocUNHaP']: + stoc_classif = model == 'StocUNHaP' + if not MULTI_START: + s, ae, rae = fadin_fit( + cleandf_hr, + TIME_START, + TIME_STOP, + model=UNHaP, + events=acti_events, + T=T, + figtitle=f'{model} on ECG CDL raw actis', + figname=f'{FILEURL[-8:-6]}{model}_raw_{TIME_START}-' + f'{TIME_STOP}min_lmbd{REG}_smartinit', + verbose=True, + init=INIT, + batch_rho=BATCH_RHO, + stoc_classif=stoc_classif, + ) + if MULTI_START: + # Run multiple starts: + # moment matching mean, moment matching max, and random init. + s, ae, rae, ll = fadin_fit( + cleandf_hr, + TIME_START, + TIME_STOP, + model=model, + events=acti_events, + T=T, + figtitle=f'{model} on ECG CDL raw actis', + figname=f'{FILEURL[-8:-6]}{model}_raw_{TIME_START}-' + f'{TIME_STOP}min_lmbd{REG}_smartinit', + verbose=True, + init='moment_matching_mean', + batch_rho=BATCH_RHO, + return_ll=True + ) + results_mmmean = pd.DataFrame.from_dict({ + 'slot': [KEY], + 'init_mode': ['moment_matching_mean'], + 'log_likelihood': [ll], + 'abs_error': [ae], + 'rel_abs_error': [rae], + 'bl_noise': [s.baseline_noise.item()], + 'bl_hawkes': [s.baseline.item()], + 'alpha_hawkes': [s.alpha.item()], + 'mean': [s.kernel_params_fixed[0].item()], + 'std': [s.kernel_params_fixed[1].item()] + } + ) + + s, ae, rae, ll = fadin_fit( + cleandf_hr, + TIME_START, + TIME_STOP, + model=model, + events=acti_events, + T=T, + figtitle=f'{model} on ECG CDL raw actis', + figname=f'{FILEURL[-8:-6]}{model}_raw_{TIME_START}-' + f'{TIME_STOP}min_lmbd{REG}_smartinit', + verbose=True, + init='moment_matching_max', + batch_rho=BATCH_RHO, + return_ll=True + ) + results_mmmax = pd.DataFrame.from_dict({ + 'slot': [KEY], + 'init_mode': [INIT], + 'log_likelihood': [ll], + 'abs_error': [ae], + 'rel_abs_error': [rae], + 'bl_noise': [s.baseline_noise.item()], + 'bl_hawkes': [s.baseline.item()], + 'alpha_hawkes': [s.alpha.item()], + 'mean': [s.kernel_params_fixed[0].item()], + 'std': [s.kernel_params_fixed[1].item()] + } + ) + + s, ae, rae, ll = fadin_fit( + cleandf_hr, + TIME_START, + TIME_STOP, + model=model, + events=acti_events, + T=T, + figtitle=f'{model} on ECG CDL raw actis', + figname=f'{FILEURL[-8:-6]}{model}_raw_{TIME_START}-' + f'{TIME_STOP}min_lmbd{REG}_smartinit', + verbose=True, + init='random', + random_state=1, + batch_rho=BATCH_RHO, + return_ll=True + ) + results_random1 = pd.DataFrame(data={ + 'slot': [KEY], + 'init_mode': ['random1'], + 'log_likelihood': [ll], + 'abs_error': [ae], + 'rel_abs_error': [rae], + 'bl_noise': [s.baseline_noise.item()], + 'bl_hawkes': [s.baseline.item()], + 'alpha_hawkes': [s.alpha.item()], + 'mean': [s.kernel_params_fixed[0].item()], + 'std': [s.kernel_params_fixed[1].item()] + } + ) + s, ae, rae, ll = fadin_fit( + cleandf_hr, + TIME_START, + TIME_STOP, + model=model, + events=acti_events, + T=T, + figtitle=f'{model} on ECG CDL raw actis', + figname=f'{FILEURL[-8:-6]}{model}_raw_{TIME_START}-' + f'{TIME_STOP}min_lmbd{REG}_smartinit', + verbose=True, + init='random', + random_state=12, + batch_rho=BATCH_RHO, + return_ll=True + ) + results_random2 = pd.DataFrame(data={ + 'slot': [KEY], + 'init_mode': ['random2'], + 'log_likelihood': [ll], + 'abs_error': [ae], + 'rel_abs_error': [rae], + 'bl_noise': [s.baseline_noise.item()], + 'bl_hawkes': [s.baseline.item()], + 'alpha_hawkes': [s.alpha.item()], + 'mean': [s.kernel_params_fixed[0].item()], + 'std': [s.kernel_params_fixed[1].item()] + } + ) + df_runs = pd.concat([df_runs, results_mmmean, results_mmmax, + results_random1, results_random2], + ignore_index=True) + + ae = np.min( + [results_mmmean['abs_error'], + results_mmmax['abs_error'], + results_random1['abs_error'], + results_random2['abs_error']] + ) + rae = np.min( + [results_mmmean['rel_abs_error'], + results_mmmax['rel_abs_error'], + results_random1['rel_abs_error'], + results_random2['rel_abs_error']] + ) + + else: + # FaDIn + s, ae, rae = fadin_fit( + cleandf_hr, + TIME_START, + TIME_STOP, + model=model, + events=acti_events, + T=T, + figtitle=f'{model} on ECG CDL raw actis', + figname=f'{FILEURL[-8:-6]}{model}_raw_{TIME_START}-' + f'{TIME_STOP}min_lmbd{REG}_smartinit', + verbose=True, + init='random', + ) + model_time = time.time() + print('Model time:', model_time) + ABS_ERROR.loc[KEY, str(model)] = ae + REL_ABS_ERROR.loc[KEY, str(model)] = rae + RUNTIME.loc[KEY, str(model)] = model_time + if model == UNHaP and KEY == 0: + plot_signal_rho( + s, begin_s=206, end_s=209, + figtitle='figsignalrhozoom2', save_fig=SAVEFIG + ) + # plt.show(block=True) + # if KEY == 0: + # plot_signal_cdl(cdl_actis, clean_df, begin_s=206, end_s=209, + # figtitle='figsignalactizoom2', save_fig=SAVEFIG) + # plt.show(block=True) + +mae = ABS_ERROR.mean(axis=0) +mrae = REL_ABS_ERROR.mean(axis=0) +std_ae = ABS_ERROR.std(axis=0) +std_rae = REL_ABS_ERROR.std(axis=0) +ABS_ERROR.to_csv(f'{RESPATH}abs_error{MODEL_LIST}.csv') +REL_ABS_ERROR.to_csv(f'{RESPATH}rel_abs_error{MODEL_LIST}.csv') +RUNTIME.to_csv(f'{RESPATH}runtime{MODEL_LIST}.csv') +print('Absolute error: \n', ABS_ERROR) +print('Relative absolute error: \n', REL_ABS_ERROR) +print('Stats on absolute error: \n', ABS_ERROR.describe()) +print('Stats on relative absolute error: \n', REL_ABS_ERROR.describe()) +print('Mean Absolute Error: \n', mae, ) +print('std of Absolute Error: \n', std_ae, '\n') +print('Mean Relative Absolute Error: \n', mrae) +print('std of Relative Absolute Error: \n', std_rae, '\n') diff --git a/experiments/unhap/ecg_gait/utils_ecg.py b/experiments/unhap/ecg_gait/utils_ecg.py new file mode 100644 index 0000000..8558c97 --- /dev/null +++ b/experiments/unhap/ecg_gait/utils_ecg.py @@ -0,0 +1,785 @@ +import vitaldb +import pandas as pd +import numpy as np +import alphacsc +import plotly.graph_objects as go +from plotly.subplots import make_subplots +import matplotlib.pyplot as plt +from scipy.signal import find_peaks +from joblib import Memory +import torch + +from fadin.solver import FaDIn, UNHaP +from fadin.loss_and_gradient import discrete_ll_loss_conv +from fadin.utils.utils import projected_grid_marked +from fadin.kernels import DiscreteKernelFiniteSupport + +memory = Memory() + +# Dataload parameters +ECG_CH_NAME = "SNUADC/ECG_II" # Name of ECG vitaldb track +HR_CH_NAME = "Solar8000/HR" # Name of heart rate vitaldb track +RESAMP_FREQ = 200 # Frequency (Hz) of data downsampling 200 + +# FaDIn solver hyperparameters +KERNEL_LENGTH = 1.5 +DELTA = 0.01 +LR = 1e-3 +L = int(1 / DELTA) +KERNEL_TYPE = "truncated_gaussian" +N_ITERATIONS_FADIN = 50_000 +RANDOM_STATE = 40 # Only used if model=FaDIn or smart_init=False + +# Dictionary learning hyperparameters +N_ATOMS = 1 +N_TIMES_ATOMS = int(RESAMP_FREQ) # int(0.7 * RESAMP_FREQ) +N_ITERATIONS = 40 +SORT_ATOMS = True # Return the atoms in decreasing signal explanability +WINDOW = True # Whether the atoms are forced to be = 0 on their edges +RANDOM_SEED = 0 + +SLOTS = { + 0: ["https://api.vitaldb.net/0010.vital", 300, 305], + 2: ["https://api.vitaldb.net/0001.vital", 35, 40], + 1: ["https://api.vitaldb.net/0001.vital", 20, 25], + 3: ["https://api.vitaldb.net/0001.vital", 55, 60], + 4: ["https://api.vitaldb.net/0002.vital", 55, 60], + 5: ["https://api.vitaldb.net/0002.vital", 120, 125], # 125, 130 + 6: ["https://api.vitaldb.net/0002.vital", 150, 155], + 7: ["https://api.vitaldb.net/0003.vital", 45, 50], + 8: ["https://api.vitaldb.net/0004.vital", 55, 60], + 9: ["https://api.vitaldb.net/0004.vital", 105, 110], + 10: ["https://api.vitaldb.net/0004.vital", 120, 125], + 12: ["https://api.vitaldb.net/0006.vital", 10, 15], + 13: ["https://api.vitaldb.net/0006.vital", 25, 30], + 14: ["https://api.vitaldb.net/0007.vital", 45, 50], + 15: ["https://api.vitaldb.net/0007.vital", 50, 55], + 16: ["https://api.vitaldb.net/0007.vital", 115, 120], + 17: ["https://api.vitaldb.net/0010.vital", 120, 125], + 18: ["https://api.vitaldb.net/0010.vital", 175, 180], + 19: ["https://api.vitaldb.net/0010.vital", 230, 235], +} + + +def plot_signal_rho(solver, begin_s, end_s, figtitle="", save_fig=None): + rho = solver.param_rho[0] + print("rho", rho) + print("shape rho", rho.shape) + begin_rho = int(begin_s / solver.delta) + end_rho = int(end_s / solver.delta) + subrho = rho[begin_rho:end_rho] + figrho, axrho = plt.subplots(1, 1, figsize=(4, 3), squeeze=True) + # Rho + axrho.stem( + subrho, + label="rho", + markerfmt="go", + linefmt="g-", + basefmt="", + ) + axrho.tick_params(axis="both", which="major", labelsize="xx-large") + axrho.legend(loc="best", fontsize="xx-large") + axrho.set_xlabel("Time (seconds)", size="xx-large") + axrho.set_title("rho", size="xx-large") + if save_fig is not None: + figrho.savefig(f"{save_fig}signal-rho{figtitle}stem.svg") + figrho.tight_layout() + figrho.show() + + +def plot_signal_cdl( + z_hat, clean_df, begin_s=None, end_s=None, figtitle="", save_fig=None +): + if begin_s is None: + begin_s = 0 + if end_s is None: + end_s = len(clean_df) // RESAMP_FREQ + subdf = clean_df.iloc[ + begin_s * RESAMP_FREQ: end_s * RESAMP_FREQ + ] + subzhat = z_hat[0][0, begin_s * RESAMP_FREQ: end_s * RESAMP_FREQ] + fig1, ax1 = plt.subplots(1, 1, figsize=(4, 3), squeeze=True) + fig2, ax2 = plt.subplots(1, 1, figsize=(4, 4), squeeze=True) + # Original signal + ax2.plot( + subdf["Time"], + subdf[ECG_CH_NAME], + label="ECG", + linewidth=3, + color="black" + ) + # CDL activations + pos_index = np.argwhere(subzhat > 0).squeeze() + ax1.stem( + subdf.iloc[pos_index]["Time"], + subzhat[pos_index], + markerfmt="bo", + linefmt="b-", + basefmt="", + label="CDL activations", + ) + ax1.hlines( + y=0, + xmin=subdf.iloc[0]["Time"], + xmax=subdf.iloc[len(subdf) - 1]["Time"] + ) + ax1.tick_params(axis="both", which="major", labelsize="xx-large") + ax1.legend(loc="best", fontsize="xx-large") + ax1.set_xlabel("Time (seconds)", size="xx-large") + ax1.set_title("Dictionary activations events", size="xx-large") + fig1.savefig(f"{save_fig}signal-actis{figtitle}stem.svg") + fig1.tight_layout() + fig1.show() + + ax2.tick_params(axis="both", which="major", labelsize="xx-large") + ax2.legend(loc="best", fontsize="xx-large") + ax2.set_xlabel("Time (seconds)", size="xx-large") + ax2.set_ylabel("mV", size="xx-large") + ax2.set_title("ECG signal", size="xx-large") + fig2.savefig(f"{save_fig}signal{figtitle}.svg") + fig2.tight_layout() + fig2.show() + + +def dataload(fileurl, time_start, time_stop, verbose=False, plotfig=False): + """Load vitaldb recording, extract [time_start, time_stop] of ECG and + heart rate channels. + Parameters + ---------- + + fileurl: str + URL of file loaded by `vitaldb.VitalFile`. + + time_start: float + Beginning of time interval to load, in minutes. + If `time_start` and `time_end` are set to None, the entire recording + is loaded. + + time_end: float + End of time interval to load, in minutes. + If `time_start` and `time_end` are set to None, the entire recording + is loaded. + + verbose: bool, default=`False` + + plot: bool, default=`False` + If set to True, the loaded data is plotted. + + Returns: + T, clean_df, cleandf_hr + T (float) is the time interval length in seconds + clean_df (pandas DataFrame) is a DataFrame with time index, and 2 columns: + ECG and heart rate. + cleandf_hr (pandas DataFrame) is a DataFrame with time index and 1 column: + heart rate with no NaN value. + """ + + # Read vitaldb file + vf = vitaldb.VitalFile(fileurl) # 'https://api.vitaldb.net/0001.vital' + # Load data in DataFrame format + df = vf.to_pandas( + [ECG_CH_NAME, HR_CH_NAME], + interval=1 / RESAMP_FREQ, + return_datetime=True + ) + # Select clean time interval for event detection + if time_start is not None and time_stop is not None: + T = int(60 * (time_stop - time_start)) + clean_df = df.iloc[ + int(60 * RESAMP_FREQ * time_start): int( + RESAMP_FREQ * (60 * time_start + T) + ) + ] + else: + T = int(len(df) / RESAMP_FREQ) + clean_df = df + cleandf_hr = clean_df.loc[pd.notna(clean_df[HR_CH_NAME])][HR_CH_NAME] + + if verbose: + # Print file metadata + print(f"Recording duration: {vf.dtend - vf.dtstart:.3f} seconds") + print(f"{len(vf.trks)} tracks") + print("Recording tracks", vf.get_track_names()) + + # ECG metadata + track = vf.trks[ECG_CH_NAME] + print("ECG Channel name", track.name) + print(f"Heart rate track type: {track.type} (1:wav, 2:num, 5:str)") + print("ECG sampling rate", track.srate, "Hz") + print("Channel unit", track.unit) + + # Heart rate medatada + track_hr = vf.trks[HR_CH_NAME] + print("Heart rate channel name:", track_hr.name) + print(f"Heart rate track type: {track_hr.type} (1:wav, 2:num, 5:str)") + print("Heart rate unit", track_hr.unit) + + if plotfig: + subdf_hr = df.loc[pd.notna(df[HR_CH_NAME])] # for visual plot + fig_hr = make_subplots(rows=2, cols=1, shared_xaxes=True) + fig_hr.add_trace( + go.Scattergl( + x=subdf_hr["Time"], + y=subdf_hr[HR_CH_NAME], + name="Heart rate", + ), + row=1, + col=1, + ) + fig_hr.add_trace( + go.Scattergl( + x=clean_df["Time"], + y=clean_df[ECG_CH_NAME], + name="ECG" + ), + row=2, + col=1, + ) + fig_hr.update_layout(title=fileurl) + fig_hr.show(renderer="browser") + + return T, clean_df, cleandf_hr + + +def peak_events(clean_df, freq, marked): + """Detect heartbeats on ECG using scipy.find_peaks. + + Each peak detected in the ECG signal is considered as a event. + The mark of each event is its peak amplitude. + + Parameters + ---------- + clean_df: pandas DataFrame + Signal DataFrame. index=timestamps, column ECG_NAME contain ECG values. + freq: float + Sampling frequency of the ECG signal. + marked: boolean + If set to True, output is a numpy array containing + detected events timestamps and their associated marks. + If set to False, output is a numpy array containing + detected events timestamps. + + Returns + ------- + events_t: numpy array + Timestamps of detected events in seconds. + Activation events on which TPP solver will be fitted. + If `marked=True`, events are listed as their timestamps. + If `marked=False`, events are listed as [timestamp, mark]. + + """ + # Detect ECG peaks + peaks_ind, props = find_peaks(clean_df[ECG_CH_NAME], height=0.1) + + if marked: + marks = props["peak_heights"] + return [ + np.array( + [[peaks_ind[i] / freq, marks[i]] + for i in range(len(peaks_ind))] + ) + ] + else: + # Convert index to time in seconds + events_t = np.reshape(peaks_ind / freq, (1, len(peaks_ind))) + return events_t + + +def naive_events(clean_df, time_start, thresh=0.5, plotfig=False): + """Detect heartbeats on ECG using a naive thresholding method. + Each time block where signal > thresh is considered as a heartbeat. + + Parameters + ---------- + clean_df: pandas DataFrame + Signal DataFrame. index=timestamps, column ECG_NAME contain ECG values. + + time_start: float + Beginning time of events time interval, in minutes. + + thresh: float, default=0.5 + Threshold value for heartbeat detection, in mV. + + plot: boolean, default=False + If set to True, ECG and detected events are plotted. + + Returns + ------- + events_t: numpy array + Timestamps of detected events in seconds. + """ + # Detect ECG peaks, ie where ECG > threshold + clean_df["discretized"] = np.int32( + np.where(clean_df[ECG_CH_NAME] > thresh, True, False) + ) + # Mark beginning of ECG peaks for FaDIn events + # When contiguous indexes are < threshold, only keep first index of block + clean_df["events"] = np.roll( + np.int32( + np.where( + np.diff(clean_df["discretized"], append=np.array([0])) == 1.0, + True, + False, + ) + ), + 1, + ) + # Retrieve timestamps of ECG peaks + events_ind = clean_df.loc[clean_df["events"] == 1].index.to_numpy() + events_t = np.reshape( + events_ind / RESAMP_FREQ - 60 * time_start, (1, len(events_ind)) + ) # seconds + + if plotfig: + # Plot ECG and detected events + fig2 = go.Figure( + go.Scattergl( + x=clean_df["Time"], + y=clean_df[ECG_CH_NAME], + name=ECG_CH_NAME + ) + ) + fig2.add_trace( + go.Scattergl( + x=clean_df["Time"], + y=clean_df["discretized"], + name="discretized" + ) + ) + fig2.add_trace( + go.Scattergl( + x=clean_df["Time"], + y=clean_df["events"], + name="events", + line_color="red", + ) + ) + fig2.show(renderer="browser") + + return events_t + + +def process_cdl_actis( + z_hat, + freq, + reg, + threshold=0.0, + sparsify=False, + marked=False, + plothist=False, + save_fig=None, +): + """Turns activations into arguments for FaDIn solvers. + Events with activation < `threshold` will be discarded. + If sparsify=`True`, consecutive activation blocks will be discarded and the + first activation of each activation block will be kept. + + Parameters + ---------- + z_hat: numpy array + activation vector of learned dictionary. + + freq: `float` or `int` + sampling frequency of z_hat. + + reg: float + Regularization parameter of CDL. + + threshold: `float` (default 0) + Threshold to apply to activations. + Only activations > threshold will be kept. + + sparsify: `bool` (default `False`) + How to handle blocks of consecutive non-zero activations. + if set to `True`, only the first activation of each block is kept. + If set to `False`, blocks are left untouched. + + marked: `bool` (default `False`) + If set to `False`, output is a numpy array containing + processed activations timestamps, suitable for FaDIn. + If set to `True`, output is a numpy array containing + processed activations timestamps and their associated marks, + suitable for MarkedFaDIn. + + plothist: `bool` (default `False`) + Whether activation blocks length histogram is plotted. The length + in this histogram are computed after thresholding and before + sparsifying. + + save_fig: string or None, default=None. + Path where to save histogram plot. + If set to None, histogram is not saved. + + Returns + ------- + events: numpy array of dimension (1, n_events) + Activation events on which FaDIn solver will be fitted. + If `marked=True`, events are listed as their timestamps. + If `marked=False`, events are listed as [timestamp, mark]. + """ + z = z_hat[0][0].copy() + # Normalize activations by max absolute value to have marks in [0, 1] + z = np.abs(z) / np.abs(z).max() + + # Threshold activations + actis_bool = np.int32(z > threshold) + not_actis = np.where(actis_bool == 0)[0] + acti_blocks_len = np.diff(not_actis) - 1 + figactis, axactis = plt.subplots() + axactis.hist( + acti_blocks_len, + bins=0.5 + np.arange(np.max(acti_blocks_len)) + ) + axactis.set_title( + f"CDL acti blocks lengths, $\\lambda$={reg}, thresh={threshold}" + ) + axactis.set_xlabel("Block length") + axactis.set_xlim([0, 10]) + if save_fig is not None: + figactis.savefig( + f"{save_fig}CDL_blocks_lengths_lmbd{reg}_thresh{threshold}.svg" + ) + if plothist: + plt.show() + + if sparsify: + # Only keep first activation of each successive activation block + actis_bool = np.diff(actis_bool, append=np.array([0])) + actis_bool = np.roll(actis_bool, 1) + # Recover index of remaining activations + events_ind = np.where(actis_bool == 1)[0] + + if marked: + return [np.array([[i / freq, z[i]] for i in events_ind])] + else: + # Convert index to time in seconds + events_t = np.reshape(events_ind / freq, (1, len(events_ind))) + return events_t + + +@memory.cache +def cdl(clean_df, reg, key, plotfig=False, save_fig=None): + """ + Learn convolutional dictionary on ECG data. A solver can then + estimate MHP parameters using dictionary activations as events. + + Parameters + ---------- + clean_df: pandas DataFrame + Signal DataFrame. index=timestamps, column ECG_NAME contain ECG values. + + reg: float + Regularization parameter of CDL. + + plotfig: boolean, default=False + If True, learned temporal atom, CDL reconstructed signal, activations, + and original signal are plotted. + + save_fig: string or None, default=None. + Path where to save figure if plotfig=true. If set to None, figure will + not be saved. + + Returns + ------- + cd.z_hat_: numpy array + Activation vector of learned dictionary. + """ + # Reshape timeseries for alphaCSC + cd_timeseries = ( + clean_df[ECG_CH_NAME].to_numpy().reshape(1, 1, -1).astype(np.float32) + ) + + # CDL solver + + cd = alphacsc.BatchCDL( + n_atoms=N_ATOMS, + n_times_atom=N_TIMES_ATOMS, + rank1=True, + n_iter=N_ITERATIONS, + n_jobs=1, + reg=reg, + sort_atoms=SORT_ATOMS, + window=WINDOW, + random_state=RANDOM_SEED, + lmbd_max="fixed", + verbose=1, + ) + cd.fit(cd_timeseries) + + if plotfig: + # Display learned atoms + fig, axes = plt.subplots( + N_ATOMS, 1, figsize=(4, 4 * N_ATOMS), squeeze=False + ) + for k in range(N_ATOMS): + axes[k, 0].plot( + np.arange(len(cd.v_hat_[k])) / RESAMP_FREQ, + cd.v_hat_[k], + linewidth=3, + color="blue", + ) + axes[N_ATOMS - 1, 0].set_xlabel("Time(seconds)", size="xx-large") + axes[0, 0].set_title("Temporal atom", size="xx-large") + axes[0, 0].tick_params( + axis="both", which="major", labelsize="xx-large") + axes[0, 0].set_ylabel("mV", size="xx-large") + if save_fig is not None: + fig.savefig(f"{save_fig}temp_atom_lmbd{reg}_key{key}.svg") + fig.show() + + # Plot activations parallelled with data + # Reconstructed signal + rec_sig = alphacsc.utils.convolution.construct_X_multi( + cd.z_hat_, cd.D_hat_, n_channels=1 + ) + z_hat_thresh = cd.z_hat_.copy() + z_hat_thresh[z_hat_thresh < 0.5] = 0 + rec_sig_thresh = alphacsc.utils.convolution.construct_X_multi( + z_hat_thresh, cd.D_hat_, n_channels=1 + ) + # Plots + # Original signal + fig3 = go.Figure( + go.Scattergl( + x=clean_df["Time"], + y=clean_df[ECG_CH_NAME], + name=ECG_CH_NAME, + opacity=0.4, + ) + ) + # CDL activations + fig3.add_trace( + go.Scattergl( + x=clean_df["Time"], + y=cd.z_hat_[0][0], + name="CDL activations", + line_color="green", + opacity=0.4, + ) + ) + # CDL reconstructed signal + fig3.add_trace( + go.Scattergl( + x=clean_df["Time"], + y=rec_sig[0][0], + name="CDL reconstructed signal", + line_color="blue", + opacity=0.4, + line={"dash": "dash"}, + ) + ) + # CDL reconstructed signal from thresholded activations + fig3.add_trace( + go.Scattergl( + x=clean_df["Time"], + y=rec_sig_thresh[0][0], + name="CDL threshold reconstructed signal", + line_color="red", + opacity=0.4, + line={"dash": "dash"}, + ) + ) + fig3.show(renderer="browser") + return cd.z_hat_ + + +def evaluate_intensity(baseline, alpha, mean, sigma, delta, events_grid): + L = int(1 / delta) + TG = DiscreteKernelFiniteSupport( + delta, n_dim=1, kernel="truncated_gaussian", lower=0, upper=1 + ) + + intens = TG.intensity_eval( + torch.tensor(baseline), + torch.tensor(alpha), + [torch.Tensor(mean), torch.Tensor(sigma)], + events_grid, + torch.linspace(0, 1, L), + ) + return intens + + +@memory.cache +def fadin_fit( + cleandf_hr, + TIME_START, + TIME_STOP, + model, + events, + T, + kernel_type=KERNEL_TYPE, + kernel_length=KERNEL_LENGTH, + lr=LR, + delta=DELTA, + max_iter=N_ITERATIONS_FADIN, + random_state=RANDOM_STATE, + init="moment_matching_max", + figtitle="", + figname="0", + verbose=False, + return_ll=False, + **fadin_init, +): + """ + Builds FaDIn solver and estimates its parameters on events. + Parameters: + ----------- + cleandf_hr: pandas DataFrame + DataFrame with time index and 1 column: heart rate with no NaN value. + TIME_START: float + Beginning time of events time interval, in minutes. + TIME_STOP: float + End time of events time interval, in minutes. + model: class instance + Type of FaDIn solver. Supported values are FaDIn, MarkedFaDIn, + JointFaDIn and JointFaDInDenoising. + + events: list + Events on which FaDIn solver will be fitted. List of length n_dim. + Each list element is a numpy array containing events for one dimension. + If the solver is non-marked, events are listed as their timestamps. + If the solver is marked, events are listed as [timestamp, mark]. + + kernel_type: str + Type of kernel for FaDIn solver. Supported values are + 'truncated_gaussian', 'raised_cosine', and 'truncated_gaussian'. + + lr: float + Solver learning rate. + + delta: float + Time delta for solver time discretization. + + T: int + Length of events time interval. + + max_iter: int + Max number of iterations of solver. + + random_state: int + Random state for solver parameters intialization. + Only used if `moment_matching=False`. + + moment_matching: boolean, default=`False` + `moment_matching` parameter of solver. + + figtitle: str + Title of kernel figure. + + figname: str + Name of file for saving kernel figure, if SAVEFIG is not None. + + verbose: boolean, `default=False` + + return_ll: boolean, `default=False` + If `True`, computes the log-likelihood of the solver on the events. + + **fadin_init: + Additional parameters of solver constructor. + """ + # Solver initialization + if model == FaDIn: + solver = model( + n_dim=1, + kernel=kernel_type, + kernel_length=kernel_length, + params_optim={"lr": lr}, + delta=delta, + max_iter=max_iter, + random_state=random_state, + **fadin_init, + ) + else: + solver = model( + n_dim=1, + kernel=kernel_type, + kernel_length=kernel_length, + params_optim={"lr": lr}, + delta=delta, + max_iter=max_iter, + random_state=random_state, + init=init, + **fadin_init, + ) + # Fit + solver.fit(events, T) + + # Print estimated solver parameters + estimated_baseline = solver.baseline + estimated_alpha = solver.alpha + param_kernel = [ + solver.param_kernel[0][-1].item(), + solver.param_kernel[1][-1].item(), + ] # [mean, std] + + if param_kernel[0] == 0: + abs_error = np.nan + else: + abs_error = np.abs(60 / param_kernel[0] - np.mean(cleandf_hr)) + rel_abs_error = abs_error / np.mean(cleandf_hr) + + if return_ll: + # Compute solver intensity function on events + if model in [UNHaP]: + intens_bl = estimated_baseline + solver.baseline_noise + else: + intens_bl = estimated_baseline + _, events_grid = projected_grid_marked(events, delta, L * T + 1) + intens = evaluate_intensity( + [intens_bl], + solver.alpha, + solver.param_kernel[0][-1], + solver.param_kernel[1][-1], + delta, + events_grid, + ) + # Compute log-likelihood of intensity function on events + solver_ll = discrete_ll_loss_conv(intens, events_grid, delta, T) + if verbose: + print("Data time interval:", TIME_START, " - ", TIME_STOP, "minutes") + # Print estimated parameters + print( + "Estimated baseline is:", + torch.round(estimated_baseline, decimals=4).item() + ) + if model in [UNHaP]: + print( + "Estimated noise baseline is:", + torch.round(solver.baseline_noise, decimals=4).item(), + ) + print( + "Estimated alpha is:", + torch.round(estimated_alpha, decimals=4).item() + ) + print( + "Estimated mean of gaussian kernel:", + np.round(param_kernel[0], decimals=4) + ) + print( + "Estimated std of gaussian kernel:", + np.round(param_kernel[1], decimals=4) + ) + # Print monitor QRS interval for comparison + print( + "Mean monitored QRS interval:", + np.round(np.mean(60 / cleandf_hr), decimals=4), + "s", + ) + print( + "std of monitored QRS interval:", + np.round(np.std(60 / cleandf_hr), decimals=4), + "s", + ) + print("Heart Rate Absolute Error", + np.round(abs_error, decimals=4), + "beats/min") + print( + "Heart Rate Relative Absolute Error", + np.round(rel_abs_error, decimals=4), + "\n", + ) + + if not return_ll: + return solver, abs_error, rel_abs_error + if return_ll: + return solver, abs_error, rel_abs_error, solver_ll From 77046f495fa8c124b8f0b559798f3d71ff1e0736 Mon Sep 17 00:00:00 2001 From: vloison Date: Tue, 8 Jul 2025 17:59:20 +0200 Subject: [PATCH 10/41] add utils --- fadin/loss_and_gradient.py | 20 ++++++++++++++++++++ fadin/utils/utils.py | 30 ++++++++++++++++++++++++++++++ 2 files changed, 50 insertions(+) diff --git a/fadin/loss_and_gradient.py b/fadin/loss_and_gradient.py index 2fa487d..decf6d9 100644 --- a/fadin/loss_and_gradient.py +++ b/fadin/loss_and_gradient.py @@ -151,6 +151,26 @@ def discrete_l2_loss_precomputation(zG, zN, ztzG, baseline, alpha, kernel, return loss_precomp / n_events.sum() +def discrete_ll_loss_conv(intensity, events_grid, delta, end_time): + """Compute the LL discrete loss using convolutions. + + Parameters + ---------- + intensity : tensor, shape (n_dim, n_grid) + Values of the intensity function evaluated on the grid. + + events_grid : tensor, shape (n_dim, n_grid) + Events projected on the pre-defined grid. + + delta : float + Step size of the discretization grid. + """ + mask = events_grid > 0 + intens = torch.log(intensity[mask]) + return (intensity.sum(1) * delta - + intens.sum()).sum() / end_time + + def squared_compensator_1(baseline): """Compute the value of the first term of the discrete l2 loss using precomputations diff --git a/fadin/utils/utils.py b/fadin/utils/utils.py index ff6a1a1..d8f8fb1 100644 --- a/fadin/utils/utils.py +++ b/fadin/utils/utils.py @@ -86,6 +86,36 @@ def projected_grid(events, grid_step, n_grid): return events_grid +def projected_grid_marked(events, delta, n_grid): + """Project the marked events on the defined grid with their associated mark + + Return a torch tensor with values of the mark on the discretized grid. + If several events are projected on the same step of the grid, + marks are added together. + """ + n_dim = len(events) + size_discret = int(1 / delta) + marks_grid = torch.zeros(n_dim, n_grid) + events_grid = torch.zeros(n_dim, n_grid) + for i in range(n_dim): + temp = np.round(events[i][:, 0] / delta) * delta + temp2 = np.round(temp * size_discret) + idx, a, _, count = np.unique( + temp2, return_counts=True, return_index=True, return_inverse=True) + + marks_grid[i, idx.astype(int)] += torch.tensor(events[i][a, 1]) + events_grid[i, idx.astype(int)] += torch.tensor(count) + # sum marked values when more than one events are projected on the + # same element of the grid + u = np.where(count > 1)[0] + for j in range(len(u)): + id = u[j] + marks_grid[i, int(idx[id])] += torch.tensor( + events[i][a[u[j]]+1:a[u[j]]+count[u[j]], 1]).sum() + + return marks_grid, events_grid + + def smooth_projection_marked(marked_events, delta, n_grid): """Project events on the grid and remove duplica in both events grid and the original events lists. Return also the marks on the grid.""" From d6b29b8b60243ad88db95c3fea0163a56052b2dd Mon Sep 17 00:00:00 2001 From: vloison Date: Tue, 8 Jul 2025 17:59:36 +0200 Subject: [PATCH 11/41] add stocunhap option to unhap --- fadin/solver.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/fadin/solver.py b/fadin/solver.py index 997de4c..3846d48 100644 --- a/fadin/solver.py +++ b/fadin/solver.py @@ -467,7 +467,8 @@ def __init__(self, n_dim, kernel, init='random', optim_mask=None, kernel_length=1, delta=0.01, optim='RMSprop', params_optim=dict(), max_iter=2000, batch_rho=100, ztzG_approx=True, tol=10e-5, density_hawkes='linear', - density_noise='uniform', random_state=None): + density_noise='uniform', stoc_classif=False, + random_state=None): # Set discretization parameters self.delta = delta @@ -492,6 +493,7 @@ def __init__(self, n_dim, kernel, init='random', optim_mask=None, self.density_hawkes = density_hawkes self.density_noise = density_noise + self.stoc_classif = stoc_classif # Set model parameters @@ -716,8 +718,17 @@ def fit(self, events, end_time): precomputations = compute_constants_unhap( z_tilde, marks_grid, events_grid, self.param_rho, marked_quantities[1], self.L) - - self.param_rho = torch.round(self.params_mixture[0].detach()) + if self.stoc_classif is False: + # vanilla UNHaP + self.param_rho = torch.round(self.params_mixture[0].detach()) + if self.stoc_classif is True: + # StocUNHaP + random_rho = torch.rand(self.n_dim, n_grid) + self.param_rho = torch.where( + random_rho < self.params_mixture[0].data, + 1., + 0. + ) # Save and clip parameters self.params_intens[0].data = self.params_intens[0].data.clip(0) * \ From 891598dce80e072fa7803ea4bf925902d39a0b25 Mon Sep 17 00:00:00 2001 From: Virginie Loison <48970001+vloison@users.noreply.github.com> Date: Wed, 9 Jul 2025 14:36:09 +0200 Subject: [PATCH 12/41] Update README.md Co-authored-by: Thomas Moreau --- README.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/README.md b/README.md index c5e9925..dc04821 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,4 @@ -# FaDIn: Fast Discretized Inference For Hawkes Processes with General Parametric Kernels -# UNHaP: Unmixing Noise from Hawkes Processes +# FaDIn: a tool box for fast and robust inference for parametric point processes ![build](https://img.shields.io/github/actions/workflow/status/GuillaumeStaermanML/FaDIn/unit_tests.yml?event=push&style=for-the-badge) ![python version](https://img.shields.io/badge/python-3.7_|_3.8_|_3.9_|_3.10_|_3.11-blue?style=for-the-badge) From 3cb77b08fe44a9684187c644cfd393a59bb224c8 Mon Sep 17 00:00:00 2001 From: Virginie Loison <48970001+vloison@users.noreply.github.com> Date: Wed, 9 Jul 2025 14:36:21 +0200 Subject: [PATCH 13/41] Update README.md Co-authored-by: Thomas Moreau --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index dc04821..a5adc58 100644 --- a/README.md +++ b/README.md @@ -53,7 +53,7 @@ pre-commit install ``` ## Short examples -A few illustrative examples are provided in the `examples` folder of this repository: +A few illustrative examples are provided in the `examples` folder of this repository, in particular: - `plot_univariate_fadin`: simulate an univariate unmarked Hawkes process, infer Hawkes Process parameters using FaDIn, and plot inferred kernel. - `plot_multivariate_fadin`: same as `plot_univariate_fadin` but in the multivariate case. - `plot_unhap`: simulate an univariate marked Hawkes process and a marked Poisson process, infer Hawkes Process parameters using UNHaP, ald plot inferred kernels. From da89433fb4df5da3aeb1261046fb96b8666f5bf1 Mon Sep 17 00:00:00 2001 From: Virginie Loison <48970001+vloison@users.noreply.github.com> Date: Wed, 9 Jul 2025 14:36:36 +0200 Subject: [PATCH 14/41] Update fadin/loss_and_gradient.py Co-authored-by: Thomas Moreau --- fadin/loss_and_gradient.py | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/fadin/loss_and_gradient.py b/fadin/loss_and_gradient.py index decf6d9..11b853f 100644 --- a/fadin/loss_and_gradient.py +++ b/fadin/loss_and_gradient.py @@ -35,17 +35,16 @@ def compute_gradient_fadin(solver, events_grid, discretization, discretization ) - solver.v_loss[i] = \ - discrete_l2_loss_precomputation( - solver.zG, - solver.zN, - solver.ztzG, - solver.params_intens[0], - solver.params_intens[1], - kernel, n_events, - solver.delta, - end_time - ).detach() + solver.v_loss[i] = discrete_l2_loss_precomputation( + solver.zG, + solver.zN, + solver.ztzG, + solver.params_intens[0], + solver.params_intens[1], + kernel, n_events, + solver.delta, + end_time + ).detach() # Update baseline gradient solver.params_intens[0].grad = get_grad_baseline( solver.zG, From f414269d74c16d196805d65d415743cae7e27ce8 Mon Sep 17 00:00:00 2001 From: Virginie Loison <48970001+vloison@users.noreply.github.com> Date: Wed, 9 Jul 2025 14:38:01 +0200 Subject: [PATCH 15/41] Update fadin/solver.py Co-authored-by: Thomas Moreau --- fadin/solver.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fadin/solver.py b/fadin/solver.py index 3846d48..49cd2fa 100644 --- a/fadin/solver.py +++ b/fadin/solver.py @@ -718,10 +718,10 @@ def fit(self, events, end_time): precomputations = compute_constants_unhap( z_tilde, marks_grid, events_grid, self.param_rho, marked_quantities[1], self.L) - if self.stoc_classif is False: + if not self.stoc_classif: # vanilla UNHaP self.param_rho = torch.round(self.params_mixture[0].detach()) - if self.stoc_classif is True: + else: # StocUNHaP random_rho = torch.rand(self.n_dim, n_grid) self.param_rho = torch.where( From a69dcf4fbe22ae99aca83a4a50b467e16365cc31 Mon Sep 17 00:00:00 2001 From: vloison Date: Wed, 9 Jul 2025 15:12:46 +0200 Subject: [PATCH 16/41] add simulate_marked_data to utils_simu --- examples/plot_unhap.py | 94 ++++------------------ fadin/tests/test_init.py | 105 +++++-------------------- fadin/utils/utils_simu.py | 160 +++++++++++++++++++++++++++++++++----- 3 files changed, 175 insertions(+), 184 deletions(-) diff --git a/examples/plot_unhap.py b/examples/plot_unhap.py index 69a83ee..c1a5c43 100644 --- a/examples/plot_unhap.py +++ b/examples/plot_unhap.py @@ -20,15 +20,12 @@ import numpy as np import matplotlib.pyplot as plt -from fadin.utils.utils_simu import simu_marked_hawkes_cluster, custom_density -from fadin.utils.utils_simu import simu_multi_poisson +from fadin.utils.utils_simu import simulate_marked_data from fadin.solver import UNHaP -from fadin.utils.functions import identity, linear_zero_one -from fadin.utils.functions import reverse_linear_zero_one, truncated_gaussian from fadin.utils.vis import plot -# %% Fixing the parameter of the simulation setting +# %% Fix the simulation and solver parameters baseline = np.array([0.3]) baseline_noise = np.array([0.05]) @@ -37,83 +34,17 @@ sigma = np.array([[0.1]]) delta = 0.01 -end_time = 10000 +end_time = 1000 seed = 0 -max_iter = 20000 +max_iter = 2000 batch_rho = 200 -# %% Create the simulating function - - -def simulate_marked_data(baseline, baseline_noise, alpha, end_time, seed=0): - n_dim = len(baseline) - - marks_kernel = identity - marks_density = linear_zero_one - time_kernel = truncated_gaussian - - params_marks_density = dict() - # params_marks_density = dict(scale=1) - params_marks_kernel = dict(slope=1.2) - params_time_kernel = dict(mu=mu, sigma=sigma) - - marked_events, _ = simu_marked_hawkes_cluster( - end_time, - baseline, - alpha, - time_kernel, - marks_kernel, - marks_density, - params_marks_kernel=params_marks_kernel, - params_marks_density=params_marks_density, - time_kernel_length=None, - marks_kernel_length=None, - params_time_kernel=params_time_kernel, - random_state=seed, - ) - - noisy_events_ = simu_multi_poisson(end_time, [baseline_noise]) - - random_marks = [np.random.rand(noisy_events_[i].shape[0]) for i in range(n_dim)] - noisy_marks = [ - custom_density( - reverse_linear_zero_one, - dict(), - size=noisy_events_[i].shape[0], - kernel_length=1.0, - ) - for i in range(n_dim) - ] - noisy_events = [ - np.concatenate( - (noisy_events_[i].reshape(-1, 1), random_marks[i].reshape(-1, 1)), axis=1 - ) - for i in range(n_dim) - ] - - events = [ - np.concatenate((noisy_events[i], marked_events[i]), axis=0) - for i in range(n_dim) - ] - - events_cat = [events[i][events[i][:, 0].argsort()] for i in range(n_dim)] - - labels = [ - np.zeros(marked_events[i].shape[0] + noisy_events_[i].shape[0]) - for i in range(n_dim) - ] - labels[0][-marked_events[0].shape[0]:] = 1.0 - true_rho = [labels[i][events[i][:, 0].argsort()] for i in range(n_dim)] - # put the mark to one to test the impact of the marks - # events_cat[0][:, 1] = 1. - - return events_cat, noisy_marks, true_rho - +# %% Simulate Hawkes Process with truncated Gaussian kernel and Poisson noise ev, noisy_marks, true_rho = simulate_marked_data( - baseline, baseline_noise.item(), alpha, end_time, seed=0 + baseline, baseline_noise.item(), alpha, end_time, mu, sigma, seed=0 ) -# %% Apply UNHAP +# %% Initiate and fit UNHAP to the simulated events solver = UNHaP( n_dim=1, @@ -148,9 +79,14 @@ def simulate_marked_data(baseline, baseline_noise, alpha, end_time, seed=0): sum_error = error_baseline + error_baseline_noise + error_alpha + error_mu + error_sigma error_params = np.sqrt(sum_error) -print("L2 square errors of the vector of parameters is:", error_params) +print("L2 square error of the vector of parameters is:", error_params) # %% Plot estimated parameters -fig, axs = plot(solver, plotfig=False, bl_noise=True, title="UNHaP fit", savefig=None) +fig, axs = plot( + solver, + plotfig=False, + bl_noise=True, + title="UNHaP fit", + savefig=None +) plt.show(block=True) -# %% diff --git a/fadin/tests/test_init.py b/fadin/tests/test_init.py index 8972e2e..06e0106 100644 --- a/fadin/tests/test_init.py +++ b/fadin/tests/test_init.py @@ -1,100 +1,31 @@ import numpy as np -# from fadin.kernels import truncated_gaussian -from fadin.utils.functions import identity, linear_zero_one -from fadin.utils.functions import reverse_linear_zero_one, truncated_gaussian -from fadin.utils.utils_simu import simu_marked_hawkes_cluster, custom_density -from fadin.utils.utils_simu import simu_multi_poisson +from fadin.utils.utils_simu import simulate_marked_data from fadin.solver import UNHaP -baseline = np.array([0.3]) -baseline_noise = np.array([0.05]) -alpha = np.array([[1.45]]) -mu = np.array([[0.4]]) -sigma = np.array([[0.1]]) -delta = 0.01 -end_time = 1000 -seed = 0 -max_iter = 20 -batch_rho = 200 - - -# Create the simulating function -def simulate_marked_data(baseline, baseline_noise, alpha, end_time, seed=0): - n_dim = len(baseline) +def test_unhap_init(): - marks_kernel = identity - marks_density = linear_zero_one - time_kernel = truncated_gaussian + baseline = np.array([0.3]) + baseline_noise = np.array([0.05]) + alpha = np.array([[1.45]]) + mu = np.array([[0.4]]) + sigma = np.array([[0.1]]) - params_marks_density = dict() - # params_marks_density = dict(scale=1) - params_marks_kernel = dict(slope=1.2) - params_time_kernel = dict(mu=mu, sigma=sigma) + end_time = 1000 + seed = 0 + max_iter = 20 + batch_rho = 5 - marked_events, _ = simu_marked_hawkes_cluster( - end_time, - baseline, - alpha, - time_kernel, - marks_kernel, - marks_density, - params_marks_kernel=params_marks_kernel, - params_marks_density=params_marks_density, - time_kernel_length=None, - marks_kernel_length=None, - params_time_kernel=params_time_kernel, - random_state=seed, + ev, noisy_marks, true_rho = simulate_marked_data( + baseline, baseline_noise.item(), alpha, end_time, mu, sigma, seed=seed ) - - noisy_events_ = simu_multi_poisson(end_time, [baseline_noise]) - - random_marks = [np.random.rand(noisy_events_[i].shape[0]) for i in range(n_dim)] - noisy_marks = [ - custom_density( - reverse_linear_zero_one, - dict(), - size=noisy_events_[i].shape[0], - kernel_length=1.0, - ) - for i in range(n_dim) - ] - noisy_events = [ - np.concatenate( - (noisy_events_[i].reshape(-1, 1), random_marks[i].reshape(-1, 1)), axis=1 - ) - for i in range(n_dim) - ] - - events = [ - np.concatenate((noisy_events[i], marked_events[i]), axis=0) - for i in range(n_dim) - ] - - events_cat = [events[i][events[i][:, 0].argsort()] for i in range(n_dim)] - - labels = [ - np.zeros(marked_events[i].shape[0] + noisy_events_[i].shape[0]) - for i in range(n_dim) - ] - labels[0][-marked_events[0].shape[0]:] = 1.0 - true_rho = [labels[i][events[i][:, 0].argsort()] for i in range(n_dim)] - - return events_cat, noisy_marks, true_rho - - -ev, noisy_marks, true_rho = simulate_marked_data( - baseline, baseline_noise.item(), alpha, end_time, seed=0 -) - - -def test_unhap_init(): unhap_mmmax = UNHaP( n_dim=1, kernel="truncated_gaussian", init="moment_matching_max", - max_iter=max_iter + max_iter=max_iter, + batch_rho=batch_rho ) unhap_mmmax.fit(ev, end_time) @@ -102,7 +33,8 @@ def test_unhap_init(): n_dim=1, kernel="truncated_gaussian", init="moment_matching_mean", - max_iter=max_iter + max_iter=max_iter, + batch_rho=batch_rho ) unhap_mmtimean.fit(ev, end_time) @@ -110,7 +42,8 @@ def test_unhap_init(): n_dim=1, kernel="truncated_gaussian", init="random", - max_iter=max_iter + max_iter=max_iter, + batch_rho=batch_rho ) unhap_random.fit(ev, end_time) assert unhap_mmmax is not None, "UNHaP moment matching max failed" diff --git a/fadin/utils/utils_simu.py b/fadin/utils/utils_simu.py index 6aae9e9..cd2dad0 100644 --- a/fadin/utils/utils_simu.py +++ b/fadin/utils/utils_simu.py @@ -4,10 +4,16 @@ from scipy.integrate import simpson import scipy.stats as stats +from fadin.utils.functions import identity, linear_zero_one +from fadin.utils.functions import reverse_linear_zero_one, truncated_gaussian + def find_max(intensity_function, duration): """Find the maximum intensity of a function.""" - res = minimize_scalar(lambda x: -intensity_function(x), bounds=(0, duration)) + res = minimize_scalar( + lambda x: -intensity_function(x), + bounds=(0, duration) + ) return -res.fun @@ -29,7 +35,8 @@ def check_random_state(seed): def kernel_norm(x, kernel, kernel_length, params_kernel=dict()): - "Normalize the given kernel on a finite support between 0 and kernel_length" + """Normalize the given kernel on a finite support between 0 and + kernel_length""" cdf_normalization = kernel.cdf(kernel_length, **params_kernel) \ - kernel.cdf(0, **params_kernel) return kernel.pdf(x, **params_kernel) / cdf_normalization @@ -37,7 +44,8 @@ def kernel_norm(x, kernel, kernel_length, params_kernel=dict()): def compute_intensity(events, s, baseline, alpha, kernel, kernel_length, params_kernel=dict()): - "Compute the intensity function at events s giving the history of events" + """Compute the intensity function at events s giving the history of events + """ diff = [] n_dim = len(events) for i in range(n_dim): @@ -49,7 +57,9 @@ def compute_intensity(events, s, baseline, alpha, kernel, for j in range(n_dim): # param du kernel ij to extend contrib_dims[i] += alpha[i, j] *\ - kernel_norm(diff[j], kernel, kernel_length, **params_kernel).sum() + kernel_norm( + diff[j], kernel, kernel_length, **params_kernel + ).sum() intens = baseline + contrib_dims return np.array(intens) @@ -119,7 +129,8 @@ def _ppf(self, x): def simu_poisson(end_time, intensity, upper_bound=None, random_state=None): """ Simulate univariate Poisson processes on [0, end_time] with the Ogata's modified thinning algorithm. - If the intensity is a numerical value, simulate a Homegenous Poisson Process, + If the intensity is a numerical value, simulate a Homegenous Poisson + Process. If the intensity is a function, simulate an Inhomogenous Poisson Process. Parameters @@ -167,10 +178,13 @@ def simu_poisson(end_time, intensity, upper_bound=None, random_state=None): return events -def simu_multi_poisson(end_time, intensity, upper_bound=None, random_state=None): +def simu_multi_poisson(end_time, intensity, upper_bound=None, + random_state=None): """Simulate multivariate Poisson processes on [0, end_time] with - the Ogata's modified thinning algorithm by superposition of univariate processes. - If the intensity is a numerical value, simulate a Homegenous Poisson Process, + the Ogata's modified thinning algorithm by superposition of univariate + processes. + If the intensity is a numerical value, simulate a Homegenous Poisson + Process. If the intensity is a function, simulate an Inhomogenous Poisson Process. Parameters @@ -227,15 +241,18 @@ def simu_multi_poisson(end_time, intensity, upper_bound=None, random_state=None) def simu_hawkes_cluster(end_time, baseline, alpha, kernel, params_kernel=dict(), kernel_length=None, upper_bound=None, random_state=None): - """ Simulate a multivariate Hawkes process following an immigration-birth procedure. + """ Simulate a multivariate Hawkes process following an immigration-birth + procedure. Edge effects may be reduced according to the second references below. References: - Møller, J., & Rasmussen, J. G. (2006). Approximate simulation of Hawkes processes. + Møller, J., & Rasmussen, J. G. (2006). Approximate simulation of Hawkes + processes. Methodology and Computing in Applied Probability, 8, 53-64. - Møller, J., & Rasmussen, J. G. (2005). Perfect simulation of Hawkes processes. + Møller, J., & Rasmussen, J. G. (2005). Perfect simulation of Hawkes + processes. Advances in applied probability, 37(3), 629-646. Parameters @@ -251,7 +268,8 @@ def simu_hawkes_cluster(end_time, baseline, alpha, kernel, kernel: str or callable The choice of the kernel for the simulation. - String kernel available are probability distribution from scipy.stats module. + String kernel available are probability distribution from scipy.stats + module. A custom kernel can be implemented with the form kernel(x, **params). params_kernel: dict @@ -299,7 +317,11 @@ def simu_hawkes_cluster(end_time, baseline, alpha, kernel, nij = Dk[i][j].sum() C[i][j] = np.repeat(Ck[j], repeats=Dk[i][j]) Eij = custom_density( - kernel, params_kernel, size=nij, kernel_length=kernel_length) + kernel, + params_kernel, + size=nij, + kernel_length=kernel_length + ) Fij = C[i][j] + Eij Fi.append(Fij) s += Fij.shape[0] @@ -407,8 +429,10 @@ def simu_hawkes_thinning(end_time, baseline, alpha, kernel, def simu_marked_hawkes_cluster(end_time, baseline, alpha, time_kernel, marks_kernel, marks_density, params_time_kernel=dict(), - params_marks_kernel=dict(), params_marks_density=dict(), - time_kernel_length=None, marks_kernel_length=None, + params_marks_kernel=dict(), + params_marks_density=dict(), + time_kernel_length=None, + marks_kernel_length=None, random_state=None): """ Simulate a multivariate marked Hawkes process following an immigration-birth procedure. @@ -510,7 +534,9 @@ def simu_marked_hawkes_cluster(end_time, baseline, alpha, time_kernel, marks_ker nij = induced_ev[i][j].sum() - time_ev[i][j] = np.repeat(time_ev_k[j], repeats=induced_ev[i][j]) + time_ev[i][j] = np.repeat( + time_ev_k[j], repeats=induced_ev[i][j] + ) inter_arrival_ij = custom_density( time_kernel, params_time_kernel, size=nij, @@ -518,15 +544,12 @@ def simu_marked_hawkes_cluster(end_time, baseline, alpha, time_kernel, marks_ker sim_ev_ij = time_ev[i][j] + inter_arrival_ij sim_ev_i.append(sim_ev_ij) - # sim_marks_ij = mark_distribution.rvs(**params_marks_density, size=nij) # custom distrib on mark density sim_marks_ij = custom_density( marks_density, params_marks_density, size=nij, kernel_length=marks_kernel_length) - # sim_marks_ij = marksij / marksij.max() - assert (sim_marks_ij > 1).sum() == 0 sim_marks_i.append(sim_marks_ij) sim_labels_i.append(gen_labels_ij) @@ -565,3 +588,102 @@ def simu_marked_hawkes_cluster(end_time, baseline, alpha, time_kernel, marks_ker it += 1 return events, labels + + +def simulate_marked_data(baseline, baseline_noise, alpha, end_time, mu, + sigma, seed=0): + """ Simulate a marked Hawkes process with noise. + + The marked Hawkes process has a truncated_gaussian kernel. The marks + densities are `linear_zero_ones` for the Hawkes process and + `reverse_linear_zero_ones` for the noisy Poisson process. + Parameters + ---------- + baseline : array of float of size (n_dim,) + Baseline parameter of the Hawkes process. + baseline_noise : float + Baseline parameter of the noisy Poisson process. + alpha : array of float of size (n_dim, n_dim) + Weight parameter associated to the kernel function. + end_time : int | float + Duration of the event time segment. + mu : float + Mean of the truncated Gaussian kernel. + sigma : float + Standard deviation of the truncated Gaussian kernel. + seed : int, default=0 + Random seed for reproducibility. + + Returns + ------- + events_cat : list of arrays + The timestamps and the marks of the point process' events. + Timestamps are first coordinate. + noisy_marks : list of arrays + The marks of the noisy Poisson process. + true_rho : list of arrays + The labels of the events, where 0 indicates noise and 1 indicates + marked Hawkes events. + """ + n_dim = len(baseline) + + marks_kernel = identity + marks_density = linear_zero_one + time_kernel = truncated_gaussian + + params_marks_density = dict() + params_marks_kernel = dict(slope=1.2) + params_time_kernel = dict(mu=mu, sigma=sigma) + + marked_events, _ = simu_marked_hawkes_cluster( + end_time, + baseline, + alpha, + time_kernel, + marks_kernel, + marks_density, + params_marks_kernel=params_marks_kernel, + params_marks_density=params_marks_density, + time_kernel_length=None, + marks_kernel_length=None, + params_time_kernel=params_time_kernel, + random_state=seed, + ) + + noisy_events_ = simu_multi_poisson(end_time, [baseline_noise]) + + random_marks = [ + np.random.rand(noisy_events_[i].shape[0]) for i in range(n_dim) + ] + noisy_marks = [ + custom_density( + reverse_linear_zero_one, + dict(), + size=noisy_events_[i].shape[0], + kernel_length=1.0, + ) + for i in range(n_dim) + ] + noisy_events = [ + np.concatenate( + (noisy_events_[i].reshape(-1, 1), random_marks[i].reshape(-1, 1)), + axis=1 + ) + for i in range(n_dim) + ] + + events = [ + np.concatenate((noisy_events[i], marked_events[i]), axis=0) + for i in range(n_dim) + ] + + events_cat = [events[i][events[i][:, 0].argsort()] for i in range(n_dim)] + + labels = [ + np.zeros(marked_events[i].shape[0] + noisy_events_[i].shape[0]) + for i in range(n_dim) + ] + labels[0][-marked_events[0].shape[0]:] = 1.0 + true_rho = [labels[i][events[i][:, 0].argsort()] for i in range(n_dim)] + + return events_cat, noisy_marks, true_rho From 9a1eab6575e9073611a0b16946292be24b071ee8 Mon Sep 17 00:00:00 2001 From: vloison Date: Wed, 9 Jul 2025 15:46:00 +0200 Subject: [PATCH 17/41] update docstrings, add _ to solver attributes. --- fadin/init.py | 4 +- fadin/solver.py | 194 +++++++++++++++++++++++++++--------------------- 2 files changed, 113 insertions(+), 85 deletions(-) diff --git a/fadin/init.py b/fadin/init.py index 7d96bc0..25588e1 100644 --- a/fadin/init.py +++ b/fadin/init.py @@ -240,7 +240,9 @@ def init_hawkes_params_unhap(solver, init, events, events_grid, 'alpha' and 'kernel'. events: list of array of size number of timestamps, list size is self.n_dim. - events_grid: TODO + events_grid: list of array of size number of timestamps, + list size is self.n_dim. + Events projected on the time grid. n_ground_events : torch.tensor Number of ground events for each dimension end_time: float diff --git a/fadin/solver.py b/fadin/solver.py index 49cd2fa..9d4c04c 100644 --- a/fadin/solver.py +++ b/fadin/solver.py @@ -132,22 +132,22 @@ class FaDIn(object): the number of parameters. The shape of each tensor is `(n_dim, n_dim)`. - param_baseline : `tensor`, shape (max_iter, n_dim) + param_baseline_ : `tensor`, shape (max_iter, n_dim) Baseline parameter of the Hawkes process for each fit iteration. - param_baseline_noise : `tensor`, shape (max_iter, n_dim) + param_baseline_noise_ : `tensor`, shape (max_iter, n_dim) Baseline parameter of the Hawkes process for each fit iteration. - param_alpha : `tensor`, shape (max_iter, n_dim, n_dim) + param_alpha_ : `tensor`, shape (max_iter, n_dim, n_dim) Weight parameter of the Hawkes process for each fit iteration. - param_kernel : `list` of `tensor` + param_kernel_ : `list` of `tensor` list containing tensor array of kernels parameters for each fit iteration. The size of the list varies depending the number of parameters. The shape of each tensor is `(n_dim, n_dim)`. - v_loss : `tensor`, shape (n_iter) + v_loss_ : `tensor`, shape (n_iter) loss accross iterations. If no early stopping, `n_iter` is equal to `max_iter`. """ @@ -247,7 +247,7 @@ def fit(self, events, end_time): print('number of events is:', n_ground_events) n_ground_events = torch.tensor(n_ground_events) # Initialize Hawkes parameters - self.params_intens = init_hawkes_params_fadin( + self.params_intens_ = init_hawkes_params_fadin( self, self.init, events, @@ -257,7 +257,7 @@ def fit(self, events, end_time): # Initialize optimizer self.opt = optimizer_fadin( - self.params_intens, + self.params_intens_, self.params_solver, solver=self.solver ) @@ -276,18 +276,23 @@ def fit(self, events, end_time): #################################################### self.v_loss = torch.zeros(self.max_iter) - self.param_baseline = torch.zeros(self.max_iter + 1, self.n_dim) - self.param_alpha = torch.zeros(self.max_iter + 1, - self.n_dim, - self.n_dim) - self.param_kernel = torch.zeros(self.n_kernel_params, - self.max_iter + 1, - self.n_dim, self.n_dim) - self.param_baseline[0] = self.params_intens[0].detach() - self.param_alpha[0] = self.params_intens[1].detach() + self.param_baseline_ = torch.zeros(self.max_iter + 1, self.n_dim) + self.param_alpha_ = torch.zeros( + self.max_iter + 1, + self.n_dim, + self.n_dim + ) + self.param_kernel_ = torch.zeros( + self.n_kernel_params, + self.max_iter + 1, + self.n_dim, + self.n_dim + ) + self.param_baseline_[0] = self.params_intens_[0].detach() + self.param_alpha_[0] = self.params_intens_[1].detach() for i in range(self.n_kernel_params): - self.param_kernel[i, 0] = self.params_intens[2 + i].detach() + self.param_kernel_[i, 0] = self.params_intens_[2 + i].detach() #################################################### start = time.time() @@ -303,34 +308,34 @@ def fit(self, events, end_time): self.opt.step() # Save parameters - self.params_intens[0].data = self.params_intens[0].data.clip(0) * \ + self.params_intens_[0].data = self.params_intens_[0].data.clip(0) * \ self.baseline_mask - self.params_intens[1].data = self.params_intens[1].data.clip(0) * \ + self.params_intens_[1].data = self.params_intens_[1].data.clip(0) * \ self.alpha_mask - self.param_baseline[i + 1] = self.params_intens[0].detach() - self.param_alpha[i + 1] = self.params_intens[1].detach() + self.param_baseline_[i + 1] = self.params_intens_[0].detach() + self.param_alpha_[i + 1] = self.params_intens_[1].detach() for j in range(self.n_kernel_params): - self.params_intens[2 + j].data = \ - self.params_intens[2 + j].data.clip(0) - self.param_kernel[j, i + 1] = \ - self.params_intens[2 + j].detach() + self.params_intens_[2 + j].data = \ + self.params_intens_[2 + j].data.clip(0) + self.param_kernel_[j, i + 1] = \ + self.params_intens_[2 + j].detach() # Early stopping if i % 100 == 0: - error_b = torch.abs(self.param_baseline[i + 1] - - self.param_baseline[i]).max() - error_al = torch.abs(self.param_alpha[i + 1] - - self.param_alpha[i]).max() - error_k = torch.abs(self.param_kernel[0, i + 1] - - self.param_kernel[0, i]).max() + error_b = torch.abs(self.param_baseline_[i + 1] - + self.param_baseline_[i]).max() + error_al = torch.abs(self.param_alpha_[i + 1] - + self.param_alpha_[i]).max() + error_k = torch.abs(self.param_kernel_[0, i + 1] - + self.param_kernel_[0, i]).max() if error_b < self.tol and error_al < self.tol \ and error_k < self.tol: print('early stopping at iteration:', i) - self.param_baseline = self.param_baseline[:i + 1] - self.param_alpha = self.param_alpha[:i + 1] + self.param_baseline_ = self.param_baseline_[:i + 1] + self.param_alpha_ = self.param_alpha_[:i + 1] for j in range(self.n_kernel_params): - self.param_kernel[j] = self.param_kernel[j, i + 1] + self.param_kernel_[j] = self.param_kernel_[j, i + 1] break print('iterations in ', time.time() - start) @@ -341,7 +346,8 @@ class UNHaP(object): """Define the UNHaP framework for estimated mixture of Hawkes and Poisson processes. - Unhap minimizes the discretized L2 mixture loss of Hawkes and Poisson processes. + Unhap minimizes the discretized L2 mixture loss of a mixture of Hawkes and + Poisson processes. Parameters ---------- @@ -426,13 +432,19 @@ class UNHaP(object): Density of the marks of the spurious events, in ``{'reverse_linear' | 'uniform'}``. + stoc_classif: `boolean`, `default=False`. + Whether to add stochasticity in the classification step of the + optimization. + If set to `True`, the rho parameter is updated stochastically. + If set to `False`, the rho parameter is updated deterministically. + random_state : `int`, `RandomState` instance or `None`, `default=None` Set the torch seed to 'random_state'. If set to `None`, torch seed will be set to 0. Attributes ---------- - params_intens : `list` of `tensor` + params_intens_ : `list` of `tensor` List of parameters of the Hawkes process at the end of the fit. The list contains the following tensors in this order: - `baseline`: `tensor`, shape (n_dim,): Baseline parameter. @@ -443,22 +455,25 @@ class UNHaP(object): the number of parameters. The shape of each tensor is `(n_dim, n_dim)`. - param_baseline : `tensor`, shape (max_iter, n_dim) + param_baseline_ : `tensor`, shape (max_iter, n_dim) Baseline parameter of the Hawkes process for each fit iteration. - param_baseline_noise : `tensor`, shape (max_iter, n_dim) + param_baseline_noise_ : `tensor`, shape (max_iter, n_dim) Baseline parameter of the Hawkes process for each fit iteration. - param_alpha : `tensor`, shape (max_iter, n_dim, n_dim) + param_alpha_ : `tensor`, shape (max_iter, n_dim, n_dim) Weight parameter of the Hawkes process for each fit iteration. - param_kernel : `list` of `tensor` + param_kernel_ : `list` of `tensor` list containing tensor array of kernels parameters for each fit iteration. The size of the list varies depending the number of parameters. The shape of each tensor is `(n_dim, n_dim)`. - v_loss : `tensor`, shape (n_iter) + param_rho_ : `tensor`, shape (n_dim, n_grid) + Final latent variable rho of the mixture model. + + v_loss_ : `tensor`, shape (n_iter) loss accross iterations. If no early stopping, `n_iter` is equal to `max_iter`. @@ -501,7 +516,9 @@ def __init__(self, n_dim, kernel, init='random', optim_mask=None, # Set optimization masks if optim_mask is None: - optim_mask = {'baseline': None, 'baseline_noise': None, 'alpha': None} + optim_mask = { + 'baseline': None, 'baseline_noise': None, 'alpha': None + } # Set baseline optimization mask if optim_mask['baseline'] is None: self.baseline_mask = torch.ones([n_dim]) @@ -523,7 +540,6 @@ def __init__(self, n_dim, kernel, init='random', optim_mask=None, assert optim_mask['alpha'].shape == torch.Size([n_dim, n_dim]), \ "Invalid alpha_mask shape, must be (n_dim, n_dim)" self.alpha_mask = optim_mask['alpha'] - # Initialization option for Hawkes parameters s = ['random', 'moment_matching_max', 'moment_matching_mean'] if isinstance(init, str): @@ -602,17 +618,22 @@ def fit(self, events, end_time): self.params_mixture = [self.rho] # Smart initialization of solver parameters - self.params_intens = init_hawkes_params_unhap( - self, self.init, marked_events, events_grid, n_ground_events, end_time + self.params_intens_ = init_hawkes_params_unhap( + self, + self.init, + marked_events, + events_grid, + n_ground_events, + end_time ) self.opt_intens, self.opt_mixture = optimizer_unhap( - [self.params_intens, self.params_mixture], + [self.params_intens_, self.params_mixture], self.params_solver, solver=self.solver ) - self.param_rho = self.params_mixture[0].detach() + self.param_rho_ = self.params_mixture[0].detach() #################################################### # Precomputations @@ -623,7 +644,7 @@ def fit(self, events, end_time): precomputations = compute_constants_unhap(z_tilde, marks_grid, events_grid, - self.param_rho, + self.param_rho_, marked_quantities[1], self.L) print('precomput:', time.time() - start) @@ -631,22 +652,24 @@ def fit(self, events, end_time): #################################################### # save results #################################################### - self.v_loss = torch.zeros(self.max_iter) + self.v_loss_ = torch.zeros(self.max_iter) + + self.param_baseline_ = torch.zeros(self.max_iter + 1, self.n_dim) + self.param_baseline_noise_ = torch.zeros(self.max_iter + 1, self.n_dim) + self.param_alpha_ = torch.zeros( + self.max_iter + 1, self.n_dim, self.n_dim) + self.param_kernel_ = torch.zeros( + self.n_kernel_params, + self.max_iter + 1, + self.n_dim, + self.n_dim + ) - self.param_baseline = torch.zeros(self.max_iter + 1, self.n_dim) - self.param_baseline_noise = torch.zeros(self.max_iter + 1, self.n_dim) - self.param_alpha = torch.zeros(self.max_iter + 1, - self.n_dim, self.n_dim) - self.param_kernel = torch.zeros(self.n_kernel_params, - self.max_iter + 1, - self.n_dim, self.n_dim) - # self.param_rho = torch.zeros(self.max_iter + 1, self.n_dim, n_grid) - - self.param_baseline[0] = self.params_intens[0].detach() - self.param_baseline_noise[0] = self.params_intens[1].detach() - self.param_alpha[0] = self.params_intens[2].detach() + self.param_baseline_[0] = self.params_intens_[0].detach() + self.param_baseline_noise_[0] = self.params_intens_[1].detach() + self.param_alpha_[0] = self.params_intens_[2].detach() for i in range(self.n_kernel_params): - self.param_kernel[i, 0] = self.params_intens[3 + i].detach() + self.param_kernel_[i, 0] = self.params_intens_[3 + i].detach() #################################################### start = time.time() @@ -658,7 +681,7 @@ def fit(self, events, end_time): # Update kernel and grad values kernel, grad_eta = self.kernel_model.kernel_and_grad( - self.params_intens[3:3 + self.n_kernel_params], + self.params_intens_[3:3 + self.n_kernel_params], discretization ) @@ -667,25 +690,25 @@ def fit(self, events, end_time): #################################################### # Update baseline, baseline noise and alpha - self.params_intens[0].grad, self.params_intens[1].grad, \ - self.params_intens[2].grad = compute_base_gradients( + self.params_intens_[0].grad, self.params_intens_[1].grad, \ + self.params_intens_[2].grad = compute_base_gradients( precomputations, - self.params_intens, + self.params_intens_, kernel, self.delta, end_time, - self.param_rho, + self.param_rho_, marked_quantities, n_ground_events ) # Update kernel for j in range(self.n_kernel_params): - self.params_intens[3 + j].grad = \ + self.params_intens_[3 + j].grad = \ get_grad_eta_mixture( precomputations, - self.params_intens[0], - self.params_intens[2], + self.params_intens_[0], + self.params_intens_[2], kernel, grad_eta[j], self.delta, @@ -703,7 +726,7 @@ def fit(self, events, end_time): marks_grid, kernel, marked_quantities[0], - self.params_intens, + self.params_intens_, self.delta, mask_void, n_ground_events, @@ -711,40 +734,43 @@ def fit(self, events, end_time): self.opt_mixture.step() - self.params_mixture[0].data = self.params_mixture[0].data.clip(0, 1) + self.params_mixture[0].data = \ + self.params_mixture[0].data.clip(0, 1) self.params_mixture[0].data[mask_void] = 0. # round the rho z_tilde = marks_grid * torch.round(self.params_mixture[0].data) precomputations = compute_constants_unhap( z_tilde, marks_grid, events_grid, - self.param_rho, marked_quantities[1], self.L) + self.param_rho_, marked_quantities[1], self.L) if not self.stoc_classif: # vanilla UNHaP - self.param_rho = torch.round(self.params_mixture[0].detach()) + self.param_rho_ = torch.round( + self.params_mixture[0].detach() + ) else: # StocUNHaP random_rho = torch.rand(self.n_dim, n_grid) - self.param_rho = torch.where( + self.param_rho_ = torch.where( random_rho < self.params_mixture[0].data, 1., 0. ) # Save and clip parameters - self.params_intens[0].data = self.params_intens[0].data.clip(0) * \ + self.params_intens_[0].data = self.params_intens_[0].data.clip(0) * \ self.baseline_mask - self.params_intens[1].data = self.params_intens[1].data.clip(0) - self.params_intens[2].data = self.params_intens[2].data.clip(0) * \ + self.params_intens_[1].data = self.params_intens_[1].data.clip(0) + self.params_intens_[2].data = self.params_intens_[2].data.clip(0) * \ self.alpha_mask - self.param_baseline[i+1] = self.params_intens[0].detach() - self.param_alpha[i+1] = self.params_intens[2].detach() - self.param_baseline_noise[i+1] = self.params_intens[1].detach() + self.param_baseline_[i+1] = self.params_intens_[0].detach() + self.param_alpha_[i+1] = self.params_intens_[2].detach() + self.param_baseline_noise_[i+1] = self.params_intens_[1].detach() for j in range(self.n_kernel_params): - self.params_intens[3+j].data = \ - self.params_intens[3+j].data.clip(1e-3) - self.param_kernel[j, i+1] = self.params_intens[3+j].detach() + self.params_intens_[3+j].data = \ + self.params_intens_[3+j].data.clip(1e-3) + self.param_kernel_[j, i+1] = self.params_intens_[3+j].detach() print('iterations in ', time.time() - start) From a08883807dfe6392d51cc66f03960a0e4a9a5a84 Mon Sep 17 00:00:00 2001 From: vloison Date: Wed, 9 Jul 2025 18:55:42 +0200 Subject: [PATCH 18/41] update solver attributes: params_intens to private, add public baseline_, alpha_, kernel_, baseline_noise --- fadin/solver.py | 170 ++++++++++++++++++++++++++++-------------------- 1 file changed, 101 insertions(+), 69 deletions(-) diff --git a/fadin/solver.py b/fadin/solver.py index 9d4c04c..14e8ed9 100644 --- a/fadin/solver.py +++ b/fadin/solver.py @@ -122,31 +122,25 @@ class FaDIn(object): Attributes ---------- - params_intens : `list` of `tensor` - List of parameters of the Hawkes process at the end of the fit. - The list contains the following tensors in this order: - - `baseline`: `tensor`, shape (n_dim,): Baseline parameter. - - `alpha`: `tensor`, shape (n_dim, n_dim): Weight parameter. - - `kernel`: `list` of `tensor`, shape (n_dim, n_dim): - Kernel parameters. The size of the list varies depending - the number of parameters. The shape of each tensor is - `(n_dim, n_dim)`. + + baseline_ : `tensor`, shape (n_dim,) + Final baseline parameter of the Hawkes process after fitting. + alpha_ : `tensor`, shape (n_dim, n_dim) + Final weight parameter of the Hawkes process kernel after fitting. + kernel_ : `list` of `tensor` + Final kernels parameters values after fitting. param_baseline_ : `tensor`, shape (max_iter, n_dim) Baseline parameter of the Hawkes process for each fit iteration. - param_baseline_noise_ : `tensor`, shape (max_iter, n_dim) Baseline parameter of the Hawkes process for each fit iteration. - param_alpha_ : `tensor`, shape (max_iter, n_dim, n_dim) Weight parameter of the Hawkes process for each fit iteration. - param_kernel_ : `list` of `tensor` list containing tensor array of kernels parameters for each fit iteration. The size of the list varies depending the number of parameters. The shape of each tensor is `(n_dim, n_dim)`. - v_loss_ : `tensor`, shape (n_iter) loss accross iterations. If no early stopping, `n_iter` is equal to `max_iter`. @@ -166,6 +160,7 @@ def __init__(self, n_dim, kernel, init='random', optim_mask=None, self.ztzG_approx = ztzG_approx # Optimizer parameters + self.fitted = False self.kernel = kernel self.solver = optim self.max_iter = max_iter @@ -234,7 +229,6 @@ def fit(self, events, end_time): Returns ------- - TODO: attributes self : object Fitted parameters. """ @@ -247,7 +241,8 @@ def fit(self, events, end_time): print('number of events is:', n_ground_events) n_ground_events = torch.tensor(n_ground_events) # Initialize Hawkes parameters - self.params_intens_ = init_hawkes_params_fadin( + # _params_intens: [baseline, alpha, kernel_params] + self._params_intens = init_hawkes_params_fadin( self, self.init, events, @@ -257,7 +252,7 @@ def fit(self, events, end_time): # Initialize optimizer self.opt = optimizer_fadin( - self.params_intens_, + self._params_intens, self.params_solver, solver=self.solver ) @@ -288,11 +283,11 @@ def fit(self, events, end_time): self.n_dim, self.n_dim ) - self.param_baseline_[0] = self.params_intens_[0].detach() - self.param_alpha_[0] = self.params_intens_[1].detach() + self.param_baseline_[0] = self._params_intens[0].detach() + self.param_alpha_[0] = self._params_intens[1].detach() for i in range(self.n_kernel_params): - self.param_kernel_[i, 0] = self.params_intens_[2 + i].detach() + self.param_kernel_[i, 0] = self._params_intens[2 + i].detach() #################################################### start = time.time() @@ -308,17 +303,17 @@ def fit(self, events, end_time): self.opt.step() # Save parameters - self.params_intens_[0].data = self.params_intens_[0].data.clip(0) * \ + self._params_intens[0].data = self._params_intens[0].data.clip(0) * \ self.baseline_mask - self.params_intens_[1].data = self.params_intens_[1].data.clip(0) * \ + self._params_intens[1].data = self._params_intens[1].data.clip(0) * \ self.alpha_mask - self.param_baseline_[i + 1] = self.params_intens_[0].detach() - self.param_alpha_[i + 1] = self.params_intens_[1].detach() + self.param_baseline_[i + 1] = self._params_intens[0].detach() + self.param_alpha_[i + 1] = self._params_intens[1].detach() for j in range(self.n_kernel_params): - self.params_intens_[2 + j].data = \ - self.params_intens_[2 + j].data.clip(0) + self._params_intens[2 + j].data = \ + self._params_intens[2 + j].data.clip(0) self.param_kernel_[j, i + 1] = \ - self.params_intens_[2 + j].detach() + self._params_intens[2 + j].detach() # Early stopping if i % 100 == 0: @@ -338,9 +333,27 @@ def fit(self, events, end_time): self.param_kernel_[j] = self.param_kernel_[j, i + 1] break print('iterations in ', time.time() - start) - + self.fitted = True return self + @property + def alpha_(self): + """Return the fitted alpha parameter of the Hawkes process.""" + assert self.fitted, "Solver must be fitted before accessing alpha_" + return self.param_alpha_[-1] + + @property + def baseline_(self): + """Return the fitted baseline parameter of the Hawkes process.""" + assert self.fitted, "Solver must be fitted before accessing baseline_" + return self.param_baseline_[-1] + + @property + def kernel_(self): + """Return the fitted kernel parameters of the Hawkes process.""" + assert self.fitted, "Solver must be fitted before accessing kernel_" + return [self.param_kernel_[j, -1] for j in range(self.n_kernel_params)] + class UNHaP(object): """Define the UNHaP framework for estimated mixture of Hawkes and @@ -444,32 +457,27 @@ class UNHaP(object): Attributes ---------- - params_intens_ : `list` of `tensor` - List of parameters of the Hawkes process at the end of the fit. - The list contains the following tensors in this order: - - `baseline`: `tensor`, shape (n_dim,): Baseline parameter. - - `baseline_noise`: `tensor`, shape (n_dim,): Baseline noise parameter. - - `alpha`: `tensor`, shape (n_dim, n_dim): Weight parameter. - - `kernel`: `list` of `tensor`, shape (n_dim, n_dim): - Kernel parameters. The size of the list varies depending - the number of parameters. The shape of each tensor is - `(n_dim, n_dim)`. + + baseline_ : `tensor`, shape (n_dim,) + Final baseline parameter of the Hawkes process after fitting. + baseline_noise_ : `tensor`, shape (n_dim,) + Final baseline noise parameter of the Hawkes process after fitting. + alpha_ : `tensor`, shape (n_dim, n_dim) + Final weight parameter of the Hawkes process kernel after fitting. + kernel_ : `list` of `tensor` + Final kernels parameters values after fitting. param_baseline_ : `tensor`, shape (max_iter, n_dim) Baseline parameter of the Hawkes process for each fit iteration. - param_baseline_noise_ : `tensor`, shape (max_iter, n_dim) Baseline parameter of the Hawkes process for each fit iteration. - param_alpha_ : `tensor`, shape (max_iter, n_dim, n_dim) Weight parameter of the Hawkes process for each fit iteration. - param_kernel_ : `list` of `tensor` list containing tensor array of kernels parameters for each fit iteration. The size of the list varies depending the number of parameters. The shape of each tensor is `(n_dim, n_dim)`. - param_rho_ : `tensor`, shape (n_dim, n_grid) Final latent variable rho of the mixture model. @@ -492,6 +500,7 @@ def __init__(self, n_dim, kernel, init='random', optim_mask=None, self.ztzG_approx = ztzG_approx # Set optimizer parameters + self.fitted = False self.kernel = kernel self.solver = optim self.max_iter = max_iter @@ -510,10 +519,6 @@ def __init__(self, n_dim, kernel, init='random', optim_mask=None, self.density_noise = density_noise self.stoc_classif = stoc_classif - # Set model parameters - - self.baseline = torch.rand(self.n_dim) - # Set optimization masks if optim_mask is None: optim_mask = { @@ -618,7 +623,8 @@ def fit(self, events, end_time): self.params_mixture = [self.rho] # Smart initialization of solver parameters - self.params_intens_ = init_hawkes_params_unhap( + # _params_intens: [baseline, baseline_noise, alpha, kernel_params] + self._params_intens = init_hawkes_params_unhap( self, self.init, marked_events, @@ -628,7 +634,7 @@ def fit(self, events, end_time): ) self.opt_intens, self.opt_mixture = optimizer_unhap( - [self.params_intens_, self.params_mixture], + [self._params_intens, self.params_mixture], self.params_solver, solver=self.solver ) @@ -665,11 +671,11 @@ def fit(self, events, end_time): self.n_dim ) - self.param_baseline_[0] = self.params_intens_[0].detach() - self.param_baseline_noise_[0] = self.params_intens_[1].detach() - self.param_alpha_[0] = self.params_intens_[2].detach() + self.param_baseline_[0] = self._params_intens[0].detach() + self.param_baseline_noise_[0] = self._params_intens[1].detach() + self.param_alpha_[0] = self._params_intens[2].detach() for i in range(self.n_kernel_params): - self.param_kernel_[i, 0] = self.params_intens_[3 + i].detach() + self.param_kernel_[i, 0] = self._params_intens[3 + i].detach() #################################################### start = time.time() @@ -681,7 +687,7 @@ def fit(self, events, end_time): # Update kernel and grad values kernel, grad_eta = self.kernel_model.kernel_and_grad( - self.params_intens_[3:3 + self.n_kernel_params], + self._params_intens[3:3 + self.n_kernel_params], discretization ) @@ -690,10 +696,10 @@ def fit(self, events, end_time): #################################################### # Update baseline, baseline noise and alpha - self.params_intens_[0].grad, self.params_intens_[1].grad, \ - self.params_intens_[2].grad = compute_base_gradients( + self._params_intens[0].grad, self._params_intens[1].grad, \ + self._params_intens[2].grad = compute_base_gradients( precomputations, - self.params_intens_, + self._params_intens, kernel, self.delta, end_time, @@ -704,11 +710,11 @@ def fit(self, events, end_time): # Update kernel for j in range(self.n_kernel_params): - self.params_intens_[3 + j].grad = \ + self._params_intens[3 + j].grad = \ get_grad_eta_mixture( precomputations, - self.params_intens_[0], - self.params_intens_[2], + self._params_intens[0], + self._params_intens[2], kernel, grad_eta[j], self.delta, @@ -726,7 +732,7 @@ def fit(self, events, end_time): marks_grid, kernel, marked_quantities[0], - self.params_intens_, + self._params_intens, self.delta, mask_void, n_ground_events, @@ -757,21 +763,47 @@ def fit(self, events, end_time): ) # Save and clip parameters - self.params_intens_[0].data = self.params_intens_[0].data.clip(0) * \ + self._params_intens[0].data = self._params_intens[0].data.clip(0) * \ self.baseline_mask - self.params_intens_[1].data = self.params_intens_[1].data.clip(0) - self.params_intens_[2].data = self.params_intens_[2].data.clip(0) * \ + self._params_intens[1].data = self._params_intens[1].data.clip(0) + self._params_intens[2].data = self._params_intens[2].data.clip(0) * \ self.alpha_mask - self.param_baseline_[i+1] = self.params_intens_[0].detach() - self.param_alpha_[i+1] = self.params_intens_[2].detach() - self.param_baseline_noise_[i+1] = self.params_intens_[1].detach() + self.param_baseline_[i+1] = self._params_intens[0].detach() + self.param_alpha_[i+1] = self._params_intens[2].detach() + self.param_baseline_noise_[i+1] = self._params_intens[1].detach() for j in range(self.n_kernel_params): - self.params_intens_[3+j].data = \ - self.params_intens_[3+j].data.clip(1e-3) - self.param_kernel_[j, i+1] = self.params_intens_[3+j].detach() + self._params_intens[3+j].data = \ + self._params_intens[3+j].data.clip(1e-3) + self.param_kernel_[j, i+1] = self._params_intens[3+j].detach() print('iterations in ', time.time() - start) - + # Save final parameters + self.fitted = True return self + + @property + def alpha_(self): + """Return the fitted alpha parameter of the Hawkes process.""" + assert self.fitted, "Solver must be fitted before accessing alpha_" + return self.param_alpha_[-1] + + @property + def baseline_(self): + """Return the fitted baseline parameter of the Hawkes process.""" + assert self.fitted, "Solver must be fitted before accessing baseline_" + return self.param_baseline_[-1] + + @property + def baseline_noise_(self): + """Return the fitted baseline noise parameter of the Hawkes process.""" + assert self.fitted, \ + "Solver must be fitted before accessing baseline_noise_" + return self.param_baseline_noise_[-1] + + @property + def kernel_(self): + """Return the fitted kernel parameters of the Hawkes process.""" + assert self.fitted, "Solver must be fitted before accessing kernel_" + return [self.param_kernel_[j, -1] for j in range(self.n_kernel_params)] From 95481e180ae7ce4bb5e753521233550db6dd6c73 Mon Sep 17 00:00:00 2001 From: vloison Date: Wed, 9 Jul 2025 18:59:59 +0200 Subject: [PATCH 19/41] update utils examples and tests with new solver attributes --- examples/plot_multivariate_fadin.py | 19 +++++++++---------- examples/plot_unhap.py | 24 ++++++++++++------------ examples/plot_univariate_fadin.py | 20 ++++++++------------ fadin/loss_and_gradient.py | 28 ++++++++++++++-------------- fadin/tests/test_mask.py | 8 ++++---- fadin/utils/vis.py | 14 +++++++------- 6 files changed, 54 insertions(+), 59 deletions(-) diff --git a/examples/plot_multivariate_fadin.py b/examples/plot_multivariate_fadin.py index 8e1f886..d7b0e03 100644 --- a/examples/plot_multivariate_fadin.py +++ b/examples/plot_multivariate_fadin.py @@ -49,7 +49,7 @@ events = simu_hawkes_cluster(T, baseline, alpha, kernel) ############################################################################### -# Here, we apply FaDIn. +# Here, we initiate FaDIn and fit it to the simulated data. solver = FaDIn( n_dim=n_dim, @@ -61,24 +61,23 @@ ) solver.fit(events, T) -# We average on the 10 last values of the optimization. +# We can now access the estimated parameters of the model. -estimated_baseline = solver.param_baseline[-10:].mean(0) -estimated_alpha = solver.param_alpha[-10:].mean(0) -param_kernel = [solver.param_kernel[0][-10:].mean(0)] +estimated_baseline = solver.param_baseline_[-10:].mean(0) +estimated_alpha = solver.param_alpha_[-10:].mean(0) +param_kernel = [solver.param_kernel_[0][-10:].mean(0)] -print('Estimated baseline is:', estimated_baseline) -print('Estimated alpha is:', estimated_alpha) +print('Estimated baseline is:', solver.baseline_) +print('Estimated alpha is:', solver.alpha_) print('Estimated parameters of the truncated Exponential kernel is:', - param_kernel[0]) + solver.kernel_) ############################################################################### # Here, we plot the values of the estimated kernels with FaDIn. kernel = DiscreteKernelFiniteSupport(dt, n_dim, kernel='truncated_exponential', kernel_length=kernel_length) -kernel_values = kernel.kernel_eval(param_kernel, - discretization) +kernel_values = kernel.kernel_eval(solver.kernel_, discretization) plt.subplots(figsize=(12, 8)) for i in range(n_dim): diff --git a/examples/plot_unhap.py b/examples/plot_unhap.py index c1a5c43..01eed23 100644 --- a/examples/plot_unhap.py +++ b/examples/plot_unhap.py @@ -63,20 +63,20 @@ # %% Print estimated parameters -print("Estimated baseline is: ", solver.param_baseline[-10:].mean().item()) -print("Estimated alpha is: ", solver.param_alpha[-10:].mean().item()) -print("Estimated kernel mean is: ", (solver.param_kernel[0][-10:].mean().item())) -print("Estimated kernel sd is: ", solver.param_kernel[1][-10:].mean().item()) -print("Estimated noise baseline is: ", solver.param_baseline_noise[-10:].mean().item()) +print("Estimated baseline is: ", solver.baseline_.item()) +print("Estimated alpha is: ", solver.alpha_.item()) +print("Estimated kernel mean is: ", solver.kernel_[0].item()) +print("Estimated kernel sd is: ", solver.kernel_[1].item()) +print("Estimated noise baseline is: ", solver.baseline_noise_.item()) # error on params -error_baseline = (solver.param_baseline[-10:].mean().item() - baseline.item()) ** 2 -error_baseline_noise = ( - solver.param_baseline_noise[-10:].mean().item() - baseline_noise.item() +error_bl = (solver.baseline_.item() - baseline.item()) ** 2 +error_bl_noise = ( + solver.baseline_noise_.item() - baseline_noise.item() ) ** 2 -error_alpha = (solver.param_alpha[-10:].mean().item() - alpha.item()) ** 2 -error_mu = (solver.param_kernel[0][-10:].mean().item() - 0.5) ** 2 -error_sigma = (solver.param_kernel[1][-10:].mean().item() - 0.1) ** 2 -sum_error = error_baseline + error_baseline_noise + error_alpha + error_mu + error_sigma +error_alpha = (solver.alpha_.item() - alpha.item()) ** 2 +error_mu = (solver.kernel_[0].item() - mu.item()) ** 2 +error_sigma = (solver.kernel_[1].item() - sigma.item()) ** 2 +sum_error = error_bl + error_bl_noise + error_alpha + error_mu + error_sigma error_params = np.sqrt(sum_error) print("L2 square error of the vector of parameters is:", error_params) diff --git a/examples/plot_univariate_fadin.py b/examples/plot_univariate_fadin.py index 9c47e61..45b1d17 100644 --- a/examples/plot_univariate_fadin.py +++ b/examples/plot_univariate_fadin.py @@ -37,8 +37,8 @@ # Here, we set the parameters of a Hawkes process with an Exponential(1) distribution. baseline = np.array([.4]) -alpha = np.array([[0.8]]) -beta = 2. +alpha = np.array([[0.9]]) +beta = 0.8 ############################################################################### # Here, we simulate the data. @@ -60,15 +60,12 @@ ) solver.fit(events, T) -# We average on the 10 last values of the optimization. - -estimated_baseline = solver.param_baseline[-10:].mean().item() -estimated_alpha = solver.param_alpha[-10:].mean().item() -param_kernel = [solver.param_kernel[0][-10:].mean().item()] +############################################################################### +# Here, we print the estimated parameters of the Hawkes process. -print('Estimated baseline is:', estimated_baseline) -print('Estimated alpha is:', estimated_alpha) -print('Estimated beta parameter of the exponential kernel is:', param_kernel[0]) +print('Estimated baseline is:', solver.baseline_.item()) +print('Estimated alpha is:', solver.alpha_.item()) +print('Estimated beta parameter of the exponential kernel is:', solver.kernel_) ############################################################################### @@ -76,8 +73,7 @@ kernel = DiscreteKernelFiniteSupport(dt, n_dim, kernel='truncated_exponential', kernel_length=kernel_length) -kernel_values = kernel.kernel_eval([torch.Tensor([param_kernel])], - discretization) +kernel_values = kernel.kernel_eval(solver.kernel_, discretization) plt.plot(discretization[1:], kernel_values.squeeze()[1:]/kernel_length, label='FaDIn\' estimated kernel') diff --git a/fadin/loss_and_gradient.py b/fadin/loss_and_gradient.py index 11b853f..1364b6f 100644 --- a/fadin/loss_and_gradient.py +++ b/fadin/loss_and_gradient.py @@ -27,11 +27,11 @@ def compute_gradient_fadin(solver, events_grid, discretization, """ # Compute kernel and gradient kernel = solver.kernel_model.kernel_eval( - solver.params_intens[2:], + solver._params_intens[2:], discretization ) grad_theta = solver.kernel_model.grad_eval( - solver.params_intens[2:], + solver._params_intens[2:], discretization ) @@ -39,42 +39,42 @@ def compute_gradient_fadin(solver, events_grid, discretization, solver.zG, solver.zN, solver.ztzG, - solver.params_intens[0], - solver.params_intens[1], + solver._params_intens[0], + solver._params_intens[1], kernel, n_events, solver.delta, end_time ).detach() # Update baseline gradient - solver.params_intens[0].grad = get_grad_baseline( + solver._params_intens[0].grad = get_grad_baseline( solver.zG, - solver.params_intens[0], - solver.params_intens[1], + solver._params_intens[0], + solver._params_intens[1], kernel, solver.delta, n_events, end_time ) # Update alpha gradient - solver.params_intens[1].grad = get_grad_alpha( + solver._params_intens[1].grad = get_grad_alpha( solver.zG, solver.zN, solver.ztzG, - solver.params_intens[0], - solver.params_intens[1], + solver._params_intens[0], + solver._params_intens[1], kernel, solver.delta, n_events ) # Update kernel gradient for j in range(solver.n_kernel_params): - solver.params_intens[2 + j].grad = \ + solver._params_intens[2 + j].grad = \ get_grad_eta( solver.zG, solver.zN, solver.ztzG, - solver.params_intens[0], - solver.params_intens[1], + solver._params_intens[0], + solver._params_intens[1], kernel, grad_theta[j], solver.delta, @@ -570,7 +570,7 @@ def get_grad_eta_mixture(precomputations, baseline, alpha, kernel, grad_theta_ = torch.zeros(n_dim, n_dim) for m in range(n_dim): - cst = 2 * delta * (baseline[m] * square_int_hawkes[m]) # + baseline_noise[m]) + cst = 2 * delta * (baseline[m] * square_int_hawkes[m]) for n in range(n_dim): grad_theta_[m, n] = cst * alpha[m, n] * (grad_kernel[m, n] @ phi_tilde[n]) diff --git a/fadin/tests/test_mask.py b/fadin/tests/test_mask.py index 86a0346..8399870 100644 --- a/fadin/tests/test_mask.py +++ b/fadin/tests/test_mask.py @@ -23,9 +23,9 @@ def maskedsolver(events, T, kernel, max_iter=1000, ztzG_approx=False, random_state=random_state ) solver.fit(events, T) - estimated_baseline = solver.params_intens[0] - estimated_alpha = solver.params_intens[1] - param_kernel = solver.params_intens[2:] + estimated_baseline = solver.baseline_ + estimated_alpha = solver.alpha_ + param_kernel = solver.kernel_ if kernel == 'raised_cosine': # multiply alpha by 2* sigma estimated_alpha = 2 * estimated_alpha * param_kernel[1] @@ -129,4 +129,4 @@ def test_rc_mask(): ) assert torch.allclose(rc_bl, torch.Tensor([0., 0.])) assert torch.allclose(rc_alpha * torch.Tensor([[1., 1.], [0., 1.]]), - torch.zeros(2, 2)) + torch.zeros(2, 2)) diff --git a/fadin/utils/vis.py b/fadin/utils/vis.py index 07627bf..16b1004 100644 --- a/fadin/utils/vis.py +++ b/fadin/utils/vis.py @@ -47,7 +47,7 @@ def plot(solver, plotfig=False, bl_noise=False, title=None, ch_names=None, kernel=solver.kernel, kernel_length=solver.kernel_length) - kappa_values = kernel.kernel_eval(solver.params_intens[-2:], + kappa_values = kernel.kernel_eval(solver.kernel_, discretization).detach() # Plot if ch_names is None: @@ -62,9 +62,9 @@ def plot(solver, plotfig=False, bl_noise=False, title=None, ch_names=None, for j in range(solver.n_dim): # Plot baseline label = (rf'$\mu_{{{ch_names[i]}}}$=' + - f'{round(solver.baseline[i].item(), 2)}') + f'{round(solver.baseline_[i].item(), 2)}') axs[i, j].hlines( - y=solver.baseline[i].item(), + y=solver.baseline_[i].item(), xmin=0, xmax=solver.kernel_length, label=label, @@ -73,10 +73,10 @@ def plot(solver, plotfig=False, bl_noise=False, title=None, ch_names=None, ) if bl_noise: # Plot noise baseline - mutilde = round(solver.baseline_noise[i].item(), 2) + mutilde = round(solver.baseline_noise_[i].item(), 2) label = rf'$\tilde{{\mu}}_{{{ch_names[i]}}}$={mutilde}' axs[i, j].hlines( - y=solver.baseline_noise[i].item(), + y=solver.baseline_noise_[i].item(), xmin=0, xmax=solver.kernel_length, label=label, @@ -84,7 +84,7 @@ def plot(solver, plotfig=False, bl_noise=False, title=None, ch_names=None, linewidth=4 ) # Plot kernel (i, j) - phi_values = solver.alpha[i, j].item() * kappa_values[i, j, 1:] + phi_values = solver.alpha_[i, j].item() * kappa_values[i, j, 1:] axs[i, j].plot( discretization[1:], phi_values, @@ -93,7 +93,7 @@ def plot(solver, plotfig=False, bl_noise=False, title=None, ch_names=None, ) if solver.kernel == 'truncated_gaussian': # Plot mean of gaussian kernel - mean = round(solver.params_intens[-2][i, j].item(), 2) + mean = round(solver.kernel_[i, j].item(), 2) axs[i, j].vlines( x=mean, ymin=0, From ce378acc6cc4f1501e6055fbdeb864f7b4d3475a Mon Sep 17 00:00:00 2001 From: vloison Date: Wed, 9 Jul 2025 19:00:08 +0200 Subject: [PATCH 20/41] add attributes test --- fadin/tests/test_attributes.py | 60 ++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 fadin/tests/test_attributes.py diff --git a/fadin/tests/test_attributes.py b/fadin/tests/test_attributes.py new file mode 100644 index 0000000..9252ce2 --- /dev/null +++ b/fadin/tests/test_attributes.py @@ -0,0 +1,60 @@ +import numpy as np +from fadin.solver import FaDIn, UNHaP +from fadin.utils.utils_simu import simu_hawkes_cluster, simulate_marked_data + + +def test_fadin_attr(): + events = simu_hawkes_cluster( + end_time=1000, + baseline=np.array([0.4]), + alpha=np.array([[0.9]]), + kernel='expon', + params_kernel={'scale': 1/4.}, + random_state=0 + ) + solver = FaDIn( + n_dim=1, + kernel="truncated_exponential", + kernel_length=1, + delta=0.01, + optim="RMSprop", + params_optim={'lr': 1e-3}, + max_iter=20 + ) + solver.fit(events, 1000) + assert hasattr(solver, 'baseline_'), \ + "FaDIn should have baseline_ attribute" + assert hasattr(solver, 'alpha_'), "FaDIn should have alpha_ attribute" + assert hasattr(solver, 'kernel_'), "FaDIn should have kernel_ attribute" + + +def test_unhap_attr(): + baseline = np.array([0.3]) + baseline_noise = np.array([0.05]) + alpha = np.array([[1.45]]) + mu = np.array([[0.4]]) + sigma = np.array([[0.1]]) + + end_time = 1000 + seed = 0 + max_iter = 20 + batch_rho = 5 + + ev, noisy_marks, true_rho = simulate_marked_data( + baseline, baseline_noise.item(), alpha, end_time, mu, sigma, seed=seed + ) + + solver = UNHaP( + n_dim=1, + kernel="truncated_gaussian", + init="random", + max_iter=max_iter, + batch_rho=batch_rho + ) + solver.fit(ev, end_time) + assert hasattr(solver, 'baseline_'), \ + "UNHaP should have baseline_ attribute" + assert hasattr(solver, 'baseline_noise_'), \ + "UNHaP should have baseline_noise_ attribute" + assert hasattr(solver, 'alpha_'), "UNHaP should have alpha_ attribute" + assert hasattr(solver, 'kernel_'), "UNHaP should have kernel_ attribute" From 24aab0e1c462993aab7ee6c76f8b0664828f4ff7 Mon Sep 17 00:00:00 2001 From: vloison Date: Wed, 9 Jul 2025 19:03:15 +0200 Subject: [PATCH 21/41] minor changes --- examples/plot_multivariate_fadin.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/examples/plot_multivariate_fadin.py b/examples/plot_multivariate_fadin.py index d7b0e03..e01b745 100644 --- a/examples/plot_multivariate_fadin.py +++ b/examples/plot_multivariate_fadin.py @@ -50,7 +50,7 @@ ############################################################################### # Here, we initiate FaDIn and fit it to the simulated data. - +print("Fitting FaDIn solver...") solver = FaDIn( n_dim=n_dim, kernel="truncated_exponential", @@ -60,13 +60,9 @@ max_iter=10000 ) solver.fit(events, T) - +print("FaDIn solver fitted.") # We can now access the estimated parameters of the model. -estimated_baseline = solver.param_baseline_[-10:].mean(0) -estimated_alpha = solver.param_alpha_[-10:].mean(0) -param_kernel = [solver.param_kernel_[0][-10:].mean(0)] - print('Estimated baseline is:', solver.baseline_) print('Estimated alpha is:', solver.alpha_) print('Estimated parameters of the truncated Exponential kernel is:', From ccd659e0e4271df48d52fa5fc19d5385697c0a3d Mon Sep 17 00:00:00 2001 From: vloison Date: Wed, 9 Jul 2025 19:10:12 +0200 Subject: [PATCH 22/41] fix plot --- fadin/utils/vis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fadin/utils/vis.py b/fadin/utils/vis.py index 16b1004..0faee76 100644 --- a/fadin/utils/vis.py +++ b/fadin/utils/vis.py @@ -93,7 +93,7 @@ def plot(solver, plotfig=False, bl_noise=False, title=None, ch_names=None, ) if solver.kernel == 'truncated_gaussian': # Plot mean of gaussian kernel - mean = round(solver.kernel_[i, j].item(), 2) + mean = round(solver.kernel_[0].item(), 2) axs[i, j].vlines( x=mean, ymin=0, From 26fecd9e05421f55071114f4a707e7d0acd2bef7 Mon Sep 17 00:00:00 2001 From: vloison Date: Wed, 9 Jul 2025 21:20:14 +0200 Subject: [PATCH 23/41] rename UNHaP rho_ attribute --- fadin/solver.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/fadin/solver.py b/fadin/solver.py index 14e8ed9..84b879c 100644 --- a/fadin/solver.py +++ b/fadin/solver.py @@ -466,6 +466,8 @@ class UNHaP(object): Final weight parameter of the Hawkes process kernel after fitting. kernel_ : `list` of `tensor` Final kernels parameters values after fitting. + rho_ : `tensor`, shape (n_dim, n_grid) + Final latent variable rho of the mixture model. param_baseline_ : `tensor`, shape (max_iter, n_dim) Baseline parameter of the Hawkes process for each fit iteration. @@ -478,8 +480,7 @@ class UNHaP(object): iteration. The size of the list varies depending the number of parameters. The shape of each tensor is `(n_dim, n_dim)`. - param_rho_ : `tensor`, shape (n_dim, n_grid) - Final latent variable rho of the mixture model. + v_loss_ : `tensor`, shape (n_iter) loss accross iterations. @@ -588,7 +589,7 @@ def fit(self, events, end_time): self : object Fitted parameters. """ - n_grid = int(1 / self.delta) * end_time + 1 + n_grid = int(1 / self.delta * end_time) + 1 discretization = torch.linspace(0, self.kernel_length, self.L) events_grid, marks_grid, marked_events, _ = \ @@ -639,7 +640,7 @@ def fit(self, events, end_time): solver=self.solver ) - self.param_rho_ = self.params_mixture[0].detach() + self.rho_ = self.params_mixture[0].detach() #################################################### # Precomputations @@ -650,7 +651,7 @@ def fit(self, events, end_time): precomputations = compute_constants_unhap(z_tilde, marks_grid, events_grid, - self.param_rho_, + self.rho_, marked_quantities[1], self.L) print('precomput:', time.time() - start) @@ -703,7 +704,7 @@ def fit(self, events, end_time): kernel, self.delta, end_time, - self.param_rho_, + self.rho_, marked_quantities, n_ground_events ) @@ -747,16 +748,16 @@ def fit(self, events, end_time): z_tilde = marks_grid * torch.round(self.params_mixture[0].data) precomputations = compute_constants_unhap( z_tilde, marks_grid, events_grid, - self.param_rho_, marked_quantities[1], self.L) + self.rho_, marked_quantities[1], self.L) if not self.stoc_classif: # vanilla UNHaP - self.param_rho_ = torch.round( + self.rho_ = torch.round( self.params_mixture[0].detach() ) else: # StocUNHaP random_rho = torch.rand(self.n_dim, n_grid) - self.param_rho_ = torch.where( + self.rho_ = torch.where( random_rho < self.params_mixture[0].data, 1., 0. From 4157a85188e1d0aa0456ebe2e5ad7d5b134616c8 Mon Sep 17 00:00:00 2001 From: vloison Date: Wed, 9 Jul 2025 21:21:20 +0200 Subject: [PATCH 24/41] update ecg experiment with attributes --- experiments/unhap/ecg_gait/utils_ecg.py | 27 +++++++++++-------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/experiments/unhap/ecg_gait/utils_ecg.py b/experiments/unhap/ecg_gait/utils_ecg.py index 8558c97..29b6c4b 100644 --- a/experiments/unhap/ecg_gait/utils_ecg.py +++ b/experiments/unhap/ecg_gait/utils_ecg.py @@ -62,7 +62,7 @@ def plot_signal_rho(solver, begin_s, end_s, figtitle="", save_fig=None): - rho = solver.param_rho[0] + rho = solver.rho_[0] print("rho", rho) print("shape rho", rho.shape) begin_rho = int(begin_s / solver.delta) @@ -705,31 +705,28 @@ def fadin_fit( solver.fit(events, T) # Print estimated solver parameters - estimated_baseline = solver.baseline - estimated_alpha = solver.alpha - param_kernel = [ - solver.param_kernel[0][-1].item(), - solver.param_kernel[1][-1].item(), - ] # [mean, std] + estimated_baseline = solver.baseline_ + estimated_alpha = solver.alpha_ + param_kernel = solver.kernel_ if param_kernel[0] == 0: abs_error = np.nan else: - abs_error = np.abs(60 / param_kernel[0] - np.mean(cleandf_hr)) + abs_error = np.abs(60 / param_kernel[0].item() - np.mean(cleandf_hr)) rel_abs_error = abs_error / np.mean(cleandf_hr) if return_ll: # Compute solver intensity function on events if model in [UNHaP]: - intens_bl = estimated_baseline + solver.baseline_noise + intens_bl = estimated_baseline + solver.baseline_noise_ else: intens_bl = estimated_baseline _, events_grid = projected_grid_marked(events, delta, L * T + 1) intens = evaluate_intensity( [intens_bl], - solver.alpha, - solver.param_kernel[0][-1], - solver.param_kernel[1][-1], + solver.alpha_, + solver.kernel_[0], + solver.kernel_[1], delta, events_grid, ) @@ -745,7 +742,7 @@ def fadin_fit( if model in [UNHaP]: print( "Estimated noise baseline is:", - torch.round(solver.baseline_noise, decimals=4).item(), + torch.round(solver.baseline_noise_, decimals=4).item(), ) print( "Estimated alpha is:", @@ -753,11 +750,11 @@ def fadin_fit( ) print( "Estimated mean of gaussian kernel:", - np.round(param_kernel[0], decimals=4) + np.round(param_kernel[0].item(), decimals=4) ) print( "Estimated std of gaussian kernel:", - np.round(param_kernel[1], decimals=4) + np.round(param_kernel[1].item(), decimals=4) ) # Print monitor QRS interval for comparison print( From 8f1cc0f86e36249f39fe5835b2065e9f3a04a908 Mon Sep 17 00:00:00 2001 From: vloison Date: Wed, 9 Jul 2025 21:21:42 +0200 Subject: [PATCH 25/41] add gait experiments --- experiments/unhap/ecg_gait/run_gait.py | 398 +++++++++++++++++++++++ experiments/unhap/ecg_gait/utils_gait.py | 170 ++++++++++ pyproject.toml | 3 + 3 files changed, 571 insertions(+) create mode 100644 experiments/unhap/ecg_gait/run_gait.py create mode 100644 experiments/unhap/ecg_gait/utils_gait.py diff --git a/experiments/unhap/ecg_gait/run_gait.py b/experiments/unhap/ecg_gait/run_gait.py new file mode 100644 index 0000000..296a6a8 --- /dev/null +++ b/experiments/unhap/ecg_gait/run_gait.py @@ -0,0 +1,398 @@ +""" +Adapted from +https://alphacsc.github.io/stable/auto_examples/dicodile/plot_gait.html#sphx-glr-auto-examples-dicodile-plot-gait-py +Code to benchmark on gait data: +- CDL + UNHaP +- CDL + FaDIn +- Neurokit +- pyHRV +""" +import time +from pathlib import Path + +import torch +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt + +import neurokit2 as nk +import pyhrv.time_domain as td +from biosppy.signals.ecg import ecg +from dicodile.data.gait import get_gait_data + +try: + from alphacsc import BatchCDL + from alphacsc.utils.convolution import construct_X_multi +except ImportError: + raise ImportError( + 'Please install the package alphacsc to run this example.\n' + 'pip install alphacsc' + ) + +from fadin.solver import FaDIn, UNHaP +from fadin.utils.vis import plot + +from fadin.utils.utils import projected_grid_marked +from utils_ecg import process_cdl_actis +from utils_gait import SLOTS +from utils_gait import save_errors +############################################################################## +# PARAMETERS +############################################################################## + +RUN_BASELINES = True + +# DATA PARAMETERS +SFREQ = 100 + +# CDL PARAMETERS +N_ATOMS = 1 +# Set individual atom (patch) size. +N_TIMES_ATOM = 150 # samples +N_ITERATIONS = 20 +REG = 0.5 + +# UNHAP AND FADIN PARAMETERS +MODEL_LIST = [FaDIn, UNHaP, 'StocUNHaP'] +KERNEL_LENGTH = 2 +DELTA = 0.01 +LR = 1e-3 +L = int(1 / DELTA) +KERNEL_TYPE = 'truncated_gaussian' +N_ITERATIONS_FADIN = 20_000 +RANDOM_STATE = 40 # Only used if 'moment_matching' in INIT +BATCH_RHO = 1000 +INIT = 'moment_matching_mean' + +# PLOT PARAMETERS +PLOTFIG = True # Whether plots are drawn. +SHOWFIG = False # Whether plots are shown at each iteration. +SAVEFIG = True # Whether plots are saved. Only read if PLOTFIG is True. +RESPATH = Path('results') / 'Gait' +RESPATH.mkdir(parents=True, exist_ok=True) + +# DataFrames to save errors +ABS_ERROR = pd.DataFrame( + data=100., + index=SLOTS.keys(), + columns=[str(m) for m in MODEL_LIST] + ['pyHRV'] + ['Neurokit'] +) +REL_ABS_ERROR = pd.DataFrame( + data=1., + index=SLOTS.keys(), + columns=[str(m) for m in MODEL_LIST] + ['pyHRV'] + ['Neurokit'] +) + +for KEY in SLOTS.keys(): + + SUBJECT = SLOTS[KEY][0] + TRIAL = SLOTS[KEY][1] + print(f'\nsubject {SUBJECT}, trial {TRIAL}') + + # Save results regularly + if TRIAL == '1': + print('Saving DataFrames') + ABS_ERROR.to_csv(RESPATH / f'abs_error{MODEL_LIST}') + REL_ABS_ERROR.to_csv(RESPATH / f'rel_abs_error{MODEL_LIST}') + + if SLOTS[KEY] in [['4', '5']]: + continue + # LOAD DATA + trial = get_gait_data(subject=SUBJECT, trial=TRIAL) # dictionary + data = trial['data'] + + ########################################################################## + # GROUND TRUTH + ########################################################################## + step_start = np.array(trial['RightFootActivity'])[:, 0] + step_delta = np.diff(step_start) + gt_mean = np.mean(step_delta) / SFREQ + print('GT inter-step interval', round(gt_mean, 3), 'seconds') + # Data keys: + # First letter: Left or Right foot. + # Second letter : A = acceleration, R = angular velocity. + # Third letter : axis (X, Y, Z or V). + # V is aligned with gravity, others are in the capteur referential. + X = trial['data']['RAV'].to_numpy() + + # Plot histogram of ground truth step delta + if PLOTFIG: + fig_gt, ax_gt = plt.subplots(layout='tight') + fig_gt.suptitle( + f's{SUBJECT} t{TRIAL} ground truth of inter-step interval' + ) + ax_gt.hist(step_delta / SFREQ) + ax_gt.set_xlim([0, KERNEL_LENGTH]) + ax_gt.set_xlabel('Time (s)') + ax_gt.set_ylabel('Histogram') + if SAVEFIG: + fig_gt.savefig(RESPATH / f's{SUBJECT}t{TRIAL}_GT.svg') + + ########################################################################## + # ECG METHODS + ########################################################################## + if RUN_BASELINES: + # PYHRV + # Get R-peaks series using biosppy + t, _, rpeaks = ecg(X, sampling_rate=SFREQ, show=False)[:3] + # Compute inter-step interval using R-peak series + nn_stats = td.nni_parameters(rpeaks=t[rpeaks]) + mean_interval_pyhrv = nn_stats['nni_mean'] / 1000 # seconds + + save_errors( + mean_interval_pyhrv, + gt_mean, + KEY, + ABS_ERROR, + REL_ABS_ERROR, + method='pyHRV' + ) + + # Neurokit + nk_start = time.time() + signals, info = nk.ecg_process( + X, + sampling_rate=SFREQ, + method='pantompkins1985' + ) + mean_interval_nk = (60. / signals['ECG_Rate']).mean() # seconds + save_errors( + mean_interval_nk, + gt_mean, + KEY, + ABS_ERROR, + REL_ABS_ERROR, + method='Neurokit' + ) + nk_time = time.time() - nk_start + print('Neurokit runtime:', round(nk_time, 3), 'seconds') + + ########################################################################## + # CDL + UNHaP + ########################################################################## + # Reshape X to (n_trials, n_channels, n_times) for alphaCSC + X = X.reshape(1, 1, *X.shape) + # DICTIONARY LEARNING + # Inialize CDL solver + cdl = BatchCDL( + # Shape of the dictionary + N_ATOMS, + N_TIMES_ATOM, + rank1=False, + n_iter=N_ITERATIONS, + n_jobs=4, + window=True, + sort_atoms=True, + random_state=20, + lmbd_max='fixed', + reg=REG, + verbose=1 + ) + # Fit cdl solver to data + res = cdl.fit(X) + D_hat = res._D_hat + + # RECONSTRUCTED SIGNAL + z_hat = res._z_hat + X_hat = construct_X_multi(z_hat, D_hat) + + if PLOTFIG: + # Plot learned atom + fig_atom, ax_atom = plt.subplots(layout='tight') + ax_atom.plot(np.arange(D_hat.shape[2]) / SFREQ, D_hat.squeeze()) + ax_atom.set_xlabel('Time (s)') + ax_atom.set_title('Temporal atom of CDL') + if SAVEFIG: + fig_atom.savefig( + RESPATH / f's{SUBJECT}t{TRIAL}_temp_atom.svg' + ) + + # POINT PROCESSES + # FaDIn + if FaDIn in MODEL_LIST: + # Initialize solver + solver_fadin = FaDIn( + n_dim=1, + kernel=KERNEL_TYPE, + kernel_length=KERNEL_LENGTH, + params_optim={'lr': LR}, + delta=DELTA, + max_iter=N_ITERATIONS_FADIN, + random_state=RANDOM_STATE, + ) + # Convert CDL events (# z_hat is normalized in process_cdl_actis) + events_nm = process_cdl_actis(z_hat, freq=SFREQ, reg=REG, marked=False) + T = int(len(z_hat.squeeze()) / SFREQ) + 1 + # fit solver to events + solver_fadin.fit(events_nm, T) + # Compute and save errors + mean_interval_fadin = solver_fadin.kernel_[0].item() + save_errors( + mean_interval_fadin, + gt_mean, + KEY, + ABS_ERROR, + REL_ABS_ERROR, + str(FaDIn) + ) + # Plot learned kernel + if PLOTFIG: + plot( + solver_fadin, + plotfig=True, + bl_noise=False, + title=f'FaDIn subject {SUBJECT}, trial {TRIAL}', + savefig=RESPATH / f's{SUBJECT}t{TRIAL}_fadin_kernel.svg' + ) + del solver_fadin + # UNHaP + if UNHaP in MODEL_LIST or 'StocUNHaP' in MODEL_LIST: + # Initialize solver + solver_unhap = UNHaP( + n_dim=1, + kernel=KERNEL_TYPE, + kernel_length=KERNEL_LENGTH, + params_optim={'lr': LR}, + delta=DELTA, + max_iter=N_ITERATIONS_FADIN, + random_state=RANDOM_STATE, + batch_rho=BATCH_RHO, + init=INIT + ) + # Convert CDL events (# z_hat is normalized in process_cdl_actis) + events = process_cdl_actis(z_hat, freq=SFREQ, reg=REG, marked=True) + T = len(z_hat.squeeze()) / SFREQ + # Fit solver to events + solver_unhap.fit(events, T) + # Compute and save errors + mean_interval_unhap = solver_unhap.kernel_[0].item() + save_errors( + mean_interval_unhap, + gt_mean, + KEY, + ABS_ERROR, + REL_ABS_ERROR, + str(UNHaP) + ) + if PLOTFIG: + # Plot learned kernel + solver_unhap.kernel_length = KERNEL_LENGTH + fig_kernel, ax_kernel = plot( + solver_unhap, + plotfig=True, + bl_noise=True, + title=f'UNHaP subject {SUBJECT}, trial {TRIAL}', + savefig=RESPATH / f's{SUBJECT}t{TRIAL}_kernel.svg' + ) + # TODO: Twin x-axis for independent y-axes + # axes = [ax_kernel, ax_kernel.twinx()] + # axes[1].hist(step_delta / SFREQ) + + # Compute PP intensity function and rho on data + n_grid = int(1 / solver_unhap.delta * T) + 1 + events_grid, events_grid_wm = projected_grid_marked( + events, solver_unhap.delta, n_grid + ) + discretization = torch.linspace( + 0, solver_unhap.kernel_length, solver_unhap.L + ) + intens = solver_unhap.kernel_model.intensity_eval( + solver_unhap.baseline_, + solver_unhap.alpha_, + solver_unhap.kernel_, + events_grid, + discretization + ).detach().numpy().squeeze() + rho = solver_unhap.rho_.numpy().squeeze() + + # Plot signal, CDL activations, Hawkes intensity + fig_hat, ax_hat = plt.subplots( + 3, figsize=(12, 8), sharex='all', layout='tight' + ) + fig_hat.suptitle(f'subject {SUBJECT}, trial {TRIAL}') + # Original and constructed signal + ax_hat[0].plot(X.squeeze(), label='ORIGINAL') + ax_hat[0].plot(X_hat.squeeze(), label='RECONSTRUCTED') + ax_hat[0].set_xlabel('time (x10ms)') + ax_hat[0].legend(loc='best') + ax_hat[0].set_title('Signal') + # Activations + ax_hat[1].stem( + z_hat[0][0], linefmt='red', label='noise activations' + ) + ax_hat[1].stem( + z_hat.squeeze() * rho[:len(z_hat.squeeze())], + linefmt='green', + label='PP activations' + ) + ax_hat[1].legend(loc='best') + ax_hat[1].set_title('CDL activations') + # Intensity + ax_hat[2].plot( + DELTA * SFREQ * np.arange(len(intens)), + rho * intens, label='Intensity function' + ) + ax_hat[2].plot( + DELTA * SFREQ * np.arange(len(rho)), + rho, label='rho' + ) + ax_hat[2].set_title('Intensity function') + ax_hat[2].legend(loc='best') + if SAVEFIG: + fig_hat.savefig( + RESPATH / f's{SUBJECT}t{TRIAL}_actis_intens.svg' + ) + del solver_unhap, cdl, res + + if 'StocUNHaP' in MODEL_LIST: + # Initialize StocUNHaP solver + solver_stoc_unhap = UNHaP( + n_dim=1, + kernel=KERNEL_TYPE, + kernel_length=KERNEL_LENGTH, + params_optim={'lr': LR}, + delta=DELTA, + max_iter=N_ITERATIONS_FADIN, + random_state=RANDOM_STATE, + init=INIT, + batch_rho=BATCH_RHO, + stoc_classif=True + ) + # Convert CDL events (# z_hat is normalized in process_cdl_actis) + events = process_cdl_actis(z_hat, freq=SFREQ, reg=REG, marked=True) + T = len(z_hat.squeeze()) / SFREQ + # Fit solver to events + solver_stoc_unhap.fit(events, T) + # Compute and save errors + mean_interval_stoc_unhap = solver_stoc_unhap.kernel_[0].item() + save_errors( + mean_interval_stoc_unhap, + gt_mean, + KEY, + ABS_ERROR, + REL_ABS_ERROR, + 'StocUNHaP' + ) + if PLOTFIG: + # Plot learned kernel + fig_kernel, ax_kernel = plot( + solver_stoc_unhap, + plotfig=True, + bl_noise=True, + title=f'StocUNHaP subject {SUBJECT}, trial {TRIAL}', + savefig=RESPATH / f's{SUBJECT}t{TRIAL}_stoc_kernel.svg' + ) + if SHOWFIG: + plt.show(block=True) + +# Compute error statistics +mae = ABS_ERROR.mean(axis=0) +mrae = REL_ABS_ERROR.mean(axis=0) +std_ae = ABS_ERROR.std(axis=0) +std_rae = REL_ABS_ERROR.std(axis=0) +print('Mean Absolute Error', mae) +print('std Absolute Error', std_ae) + +# Save metrics +ABS_ERROR.to_csv(RESPATH / f'abs_error{MODEL_LIST}.csv') +REL_ABS_ERROR.to_csv(RESPATH / f'rel_abs_error{MODEL_LIST}.csv') diff --git a/experiments/unhap/ecg_gait/utils_gait.py b/experiments/unhap/ecg_gait/utils_gait.py new file mode 100644 index 0000000..873b87c --- /dev/null +++ b/experiments/unhap/ecg_gait/utils_gait.py @@ -0,0 +1,170 @@ + +import json +from zipfile import ZipFile + +import numpy as np +import pandas as pd +from download import download + + +SLOTS = { + 1: [1, 1], + 2: [2, 1], + 3: [3, 1], + 4: [4, 1], + 5: [5, 1], + 6: [6, 1], + 7: [7, 1], + 8: [8, 1], + 9: [9, 1], + 10: [10, 1], + 11: [11, 1], + 12: [12, 1], + 13: [13, 1], + 14: [14, 1], + 15: [15, 1], + 16: [16, 1], + 17: [17, 1], + 18: [18, 1], + 19: [19, 1], + 20: [20, 1], + 21: [21, 1], + 22: [22, 1], + 23: [23, 1], + 24: [24, 1], + 25: [25, 1], + 26: [26, 1], + 27: [27, 1], + 28: [28, 1], + 29: [29, 1], + 30: [30, 1], + 40: [40, 1], + 41: [41, 1], + 42: [42, 1], + 43: [43, 1], + 44: [44, 1], + 45: [45, 1], + 46: [46, 1], + 47: [47, 1], + 48: [48, 1], + 49: [49, 1], +} +# DATA_HOME = Path("gait_data") +# GAIT_RECORD_ID_LIST_FNAME = DATA_HOME / "gait_record_id_list.json" +# GAIT_PARTICIPANTS_FNAME = DATA_HOME / "gait_participants.tsv" +# SLOTS_FNAME = DATA_HOME / "gait_slots.json" + + +# def download_gait(verbose=True): +# gait_dir = DATA_HOME / "gait" +# gait_dir.mkdir(parents=True, exist_ok=True) +# gait_zip = download( +# "http://dev.ipol.im/~truong/GaitData.zip", +# gait_dir / "GaitData.zip", +# replace=False, +# verbose=verbose +# ) + +# return gait_zip + + +def get_gait_data(subject=1, trial=1, only_meta=False, verbose=True): + """ + Retrieve gait data from this `dataset`_. + + Parameters + ---------- + subject: int, defaults to 1 + Subject identifier. + Valid subject-trial pairs can be found in this `list`_. + trial: int, defaults to 1 + Trial number. + Valid subject-trial pairs can be found in this `list`_. + only_meta: bool, default to False + If True, only returns the subject metadata + verbose : bool, default to True + Whether to print download status to the screen. + + Returns + ------- + dictDATA_HOME = Path("gait_data") +# GAIT_RECORD_ID_LIST_FNAME = DATA_HOME / "gait_record_id_list.json" +# GAIT_PARTICIPANTS_FNAME = DATA_HOME / "gait_participants.tsv" +# SLOTS_FNAME = DATA_HOME / "gait_slots.json" + + +# def download_gait(verbose=True): +# gait_dir = DATA_HOME / "gait" +# gait_dir.mkdir(parents=True, exist_ok=True) +# gait_zip = download( +# "http://dev.ipol.im/~truong/GaitData.zip", +# gait_dir / "GaitData.zip", +# replace=False, +# verbose=verbose +# ) + +# return gait_zip + + A dictionary containing metadata and data relative + to a trial. The 'data' attribute contains time + series for the trial, as a Pandas dataframe. + + + .. _dataset: https://github.com/deepcharles/gait-data + .. _list: + https://github.com/deepcharles/gait-data/blob/master/code_list.json + """ + # coerce subject and trial + subject = int(subject) + trial = int(trial) + + gait_zip = download_gait(verbose=verbose) + + with ZipFile(gait_zip) as zf: + with zf.open(f"GaitData/{subject}-{trial}.json") as meta_file, \ + zf.open(f"GaitData/{subject}-{trial}.csv") as data_file: + meta = json.load(meta_file) + if not only_meta: + data = pd.read_csv(data_file, sep=',', header=0) + meta['data'] = data + return meta + + +def save_errors(est, gt, key, ae_df, rae_df, method='', verbose=True, + unit='seconds'): + """Computes and saves Absolute Error and Relative Error between + est and gt into DataFrames. + Parameters: + ----------- + est: float + Estimated value which will be conpared to the ground truth value. + gt: float + Ground truth value. + key: Any + Key of DataFrame line in which the errors will be saved. + ae_df: pandas.DataFrame + DataFrame where the Absolute Error will be saved. Must be created + before calling this function. + rae_df: pandas.DataFrame + DataFrame where the Relative Absolute Error will be saved. + Must be created before calling this function. + method: `str` + Name of method used to compute `est`. Must be a column of dataframes + ae_df and rae_df. + verbose: `bool` + If True, the Absolute Error is printed.DATA_HOME = Path("gait_data") + + Returns: + -------- + ae: float + Absolute Error. + rae: float + Relative Absolute Error. + """ + ae = np.abs(est - gt) + rae = ae / gt + ae_df.loc[key, method] = ae + rae_df.loc[key, method] = rae + if verbose: + print(f'{method} Absolute Error:', round(ae, 3), unit) + return ae, rae diff --git a/pyproject.toml b/pyproject.toml index e6f5458..c2124cc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,6 +41,9 @@ test = [ "pytest>=7.2", ] +experiments = [ + "dicodile>=0.3", +] doc = [ "furo", "matplotlib", From eb67efc8c0352682b016e0d326f74fdc49e4102f Mon Sep 17 00:00:00 2001 From: vloison Date: Wed, 9 Jul 2025 21:40:57 +0200 Subject: [PATCH 26/41] remove unused function in gait experiment --- experiments/unhap/ecg_gait/utils_gait.py | 63 ------------------------ 1 file changed, 63 deletions(-) diff --git a/experiments/unhap/ecg_gait/utils_gait.py b/experiments/unhap/ecg_gait/utils_gait.py index 873b87c..f35c75c 100644 --- a/experiments/unhap/ecg_gait/utils_gait.py +++ b/experiments/unhap/ecg_gait/utils_gait.py @@ -4,7 +4,6 @@ import numpy as np import pandas as pd -from download import download SLOTS = { @@ -68,68 +67,6 @@ # return gait_zip -def get_gait_data(subject=1, trial=1, only_meta=False, verbose=True): - """ - Retrieve gait data from this `dataset`_. - - Parameters - ---------- - subject: int, defaults to 1 - Subject identifier. - Valid subject-trial pairs can be found in this `list`_. - trial: int, defaults to 1 - Trial number. - Valid subject-trial pairs can be found in this `list`_. - only_meta: bool, default to False - If True, only returns the subject metadata - verbose : bool, default to True - Whether to print download status to the screen. - - Returns - ------- - dictDATA_HOME = Path("gait_data") -# GAIT_RECORD_ID_LIST_FNAME = DATA_HOME / "gait_record_id_list.json" -# GAIT_PARTICIPANTS_FNAME = DATA_HOME / "gait_participants.tsv" -# SLOTS_FNAME = DATA_HOME / "gait_slots.json" - - -# def download_gait(verbose=True): -# gait_dir = DATA_HOME / "gait" -# gait_dir.mkdir(parents=True, exist_ok=True) -# gait_zip = download( -# "http://dev.ipol.im/~truong/GaitData.zip", -# gait_dir / "GaitData.zip", -# replace=False, -# verbose=verbose -# ) - -# return gait_zip - - A dictionary containing metadata and data relative - to a trial. The 'data' attribute contains time - series for the trial, as a Pandas dataframe. - - - .. _dataset: https://github.com/deepcharles/gait-data - .. _list: - https://github.com/deepcharles/gait-data/blob/master/code_list.json - """ - # coerce subject and trial - subject = int(subject) - trial = int(trial) - - gait_zip = download_gait(verbose=verbose) - - with ZipFile(gait_zip) as zf: - with zf.open(f"GaitData/{subject}-{trial}.json") as meta_file, \ - zf.open(f"GaitData/{subject}-{trial}.csv") as data_file: - meta = json.load(meta_file) - if not only_meta: - data = pd.read_csv(data_file, sep=',', header=0) - meta['data'] = data - return meta - - def save_errors(est, gt, key, ae_df, rae_df, method='', verbose=True, unit='seconds'): """Computes and saves Absolute Error and Relative Error between From a93c4d9d0a4aefb5aa36fe7d6b8005b303b4efa5 Mon Sep 17 00:00:00 2001 From: vloison Date: Wed, 9 Jul 2025 21:42:14 +0200 Subject: [PATCH 27/41] minor changes --- experiments/unhap/ecg_gait/utils_gait.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/experiments/unhap/ecg_gait/utils_gait.py b/experiments/unhap/ecg_gait/utils_gait.py index f35c75c..02e1071 100644 --- a/experiments/unhap/ecg_gait/utils_gait.py +++ b/experiments/unhap/ecg_gait/utils_gait.py @@ -1,9 +1,6 @@ -import json -from zipfile import ZipFile import numpy as np -import pandas as pd SLOTS = { From a1b05ac6534f666c69a92504641d8e46deb65085 Mon Sep 17 00:00:00 2001 From: vloison Date: Fri, 11 Jul 2025 19:06:40 +0200 Subject: [PATCH 28/41] add unhap experiments on simulated data --- .../unhap/simulated_data/plot_infer_noise.py | 77 + .../unhap/simulated_data/plot_infer_rho.py | 89 + .../simulated_data/run_benchmark_marked.py | 465 +++++ .../unhap/simulated_data/run_benchmark_tpp.py | 305 ++++ .../unhap/simulated_data/run_infer_noise.py | 1570 +++++++++++++++++ .../unhap/simulated_data/run_rho_error.py | 163 ++ 6 files changed, 2669 insertions(+) create mode 100644 experiments/unhap/simulated_data/plot_infer_noise.py create mode 100644 experiments/unhap/simulated_data/plot_infer_rho.py create mode 100644 experiments/unhap/simulated_data/run_benchmark_marked.py create mode 100644 experiments/unhap/simulated_data/run_benchmark_tpp.py create mode 100644 experiments/unhap/simulated_data/run_infer_noise.py create mode 100644 experiments/unhap/simulated_data/run_rho_error.py diff --git a/experiments/unhap/simulated_data/plot_infer_noise.py b/experiments/unhap/simulated_data/plot_infer_noise.py new file mode 100644 index 0000000..5ef1208 --- /dev/null +++ b/experiments/unhap/simulated_data/plot_infer_noise.py @@ -0,0 +1,77 @@ +# %% +import numpy as np +import pandas as pd +import matplotlib +import matplotlib.pyplot as plt +from matplotlib.lines import Line2D + + +FONTSIZE = 14 +plt.rcParams["figure.figsize"] = (5, 3.2) +plt.rcParams["axes.grid"] = False +plt.rcParams["axes.grid.axis"] = "y" +plt.rcParams["grid.linestyle"] = "--" +plt.rcParams['xtick.labelsize'] = FONTSIZE +plt.rcParams['ytick.labelsize'] = FONTSIZE +plt.rcParams['font.size'] = FONTSIZE +plt.rcParams['mathtext.fontset'] = 'cm' +plt.rc('legend', fontsize=FONTSIZE - 1) +matplotlib.rcParams['pdf.fonttype'] = 42 +matplotlib.rcParams['ps.fonttype'] = 42 + + +colors = [matplotlib.cm.viridis(x) for x in np.linspace(0, 1, 5)][1:] +setting = 'high' + +df = pd.read_csv(f'results/error_denoising_infer_unbalanced_{setting}.csv') + +T = df["end_time"].unique() +T.sort() +ls = [':', '-'] +models = ['jointfadin', 'mixture'] + +lw = 4 +# %% +for j, model in enumerate(models): + for i, t in enumerate(T): + this_df = df.query("end_time == @t") + curve = this_df.groupby("noise")[f"err_norm2_{model}"].quantile( + [0.25, 0.5, 0.75]).unstack() + plt.semilogy( + curve.index, curve[0.5], lw=lw, c=colors[2-i], + markersize=10, markevery=3, linestyle=ls[j] + ) + +plt.xlim(0.1, 1.5) +plt.xticks(fontsize=10) +plt.yticks(fontsize=10) +plt.ylabel('Estimation Error', labelpad=0) +plt.xlabel('Noise level', labelpad=0, size=12) +plt.savefig(f'plots/error_denoising_infer_unbalanced_{setting}.pdf') +plt.show() +# %% +fig1, ax1 = plt.subplots(1, 1, figsize=(8, 0.5)) +custom_lines_model = [ + Line2D([], [], color='k', lw=3, linestyle=ls[i]) for i in range(2) +] +custom_lines_T = [ + Line2D([], [], color=colors[2-j], lw=3,) for j in range(3) +] + +leg = ax1.legend( + custom_lines_model, + [f"{models[i]}" for i in range(2)], + title="Models", loc="lower center", + bbox_to_anchor=(-0.05, -0.3, 0.65, 0.01), ncol=3, prop={'size': 9} + ) +leg2 = ax1.legend( + custom_lines_T, + [f"{T[i]}" for i in range(3)], + title=r"$T$", loc="lower center", + bbox_to_anchor=(0.1, -0.3, 1.2, 0.01), ncol=3, prop={'size': 9} + ) +ax1.add_artist(leg) +ax1.add_artist(leg2) +ax1.axis("off") +plt.savefig('plots/legend_infer_noise.pdf') +# %% diff --git a/experiments/unhap/simulated_data/plot_infer_rho.py b/experiments/unhap/simulated_data/plot_infer_rho.py new file mode 100644 index 0000000..305ccf4 --- /dev/null +++ b/experiments/unhap/simulated_data/plot_infer_rho.py @@ -0,0 +1,89 @@ +# %% +import numpy as np +import pandas as pd +import matplotlib +import matplotlib.pyplot as plt + + +FONTSIZE = 14 +plt.rcParams["figure.figsize"] = (5, 3.2) +plt.rcParams["axes.grid"] = False +plt.rcParams["axes.grid.axis"] = "y" +plt.rcParams["grid.linestyle"] = "--" +plt.rcParams['xtick.labelsize'] = FONTSIZE +plt.rcParams['ytick.labelsize'] = FONTSIZE +plt.rcParams['font.size'] = FONTSIZE +plt.rcParams['mathtext.fontset'] = 'cm' +plt.rc('legend', fontsize=FONTSIZE - 1) +matplotlib.rcParams['pdf.fonttype'] = 42 +matplotlib.rcParams['ps.fonttype'] = 42 + +palette = [matplotlib.cm.viridis_r(x) for x in np.linspace(0, 1, 5)][1:] + +df = pd.read_csv('results/error_rho_infer_5000.csv') + +T = df["end_time"].unique() +T.sort() + +noise = .5 +_ = plt.figure() +for i, t in enumerate(T): + df = df[df['noise'] == noise] + this_df = df.query("end_time == @t") + curve = this_df.groupby("alpha")["rec_score"].quantile( + [0.35, 0.5, 0.65]).unstack() + plt.plot( + curve.index, curve[0.5], lw=4, c=palette[i], + markersize=10, markevery=3, label=f'$T$ = {T[i]}' + ) + plt.fill_between( + curve.index, curve[0.35], curve[0.65], alpha=0.2, + color=palette[i] + ) +plt.xlim(0, 1.47) +plt.xticks(fontsize=7) +plt.yticks(fontsize=10) +plt.ylabel(r'rc score', labelpad=0) +plt.xlabel(r'$\alpha$', labelpad=0, size=15) +plt.legend() +plt.savefig(f'plots/infer_rho_5000_{noise}.pdf') +plt.show() +# %% + + +df = pd.read_csv('results/error_rho_infer_5000.csv') + +mk = ['x', 'o', '^'] +fig, ax = plt.subplots() +noise = .1 +# Assuming T is defined somewhere in your code +for i, t in enumerate(T): + df = df[df['noise'] == noise] + this_df = df.query("end_time == @t") + curve1 = this_df.groupby("alpha")["pr_score"].quantile( + [0.25, 0.5, 0.75]).unstack() + curve2 = this_df.groupby("alpha")["rec_score"].quantile( + [0.25, 0.5, 0.75]).unstack() + + pr = curve1[0.5].values[::2] + rec = curve2[0.5].values[::2] + palette = np.linspace(0, 1.5, pr.shape[0]) + + # Scatter plot with specified marker and color + im = ax.scatter( + pr, + rec, + c=palette, + marker=mk[i], + label=f'T={t}', + cmap='plasma', + s=80 + ) + +ax.set_xlabel('Precision', size=15, labelpad=0) +ax.set_ylabel('Recall', size=15, labelpad=0) +ax.tick_params(axis='both', which='major', labelsize=7) +ax.tick_params(axis='both', which='minor', labelsize=7) + +plt.savefig('plots/scatter_rho_uniform.pdf') +plt.show() diff --git a/experiments/unhap/simulated_data/run_benchmark_marked.py b/experiments/unhap/simulated_data/run_benchmark_marked.py new file mode 100644 index 0000000..88b30b6 --- /dev/null +++ b/experiments/unhap/simulated_data/run_benchmark_marked.py @@ -0,0 +1,465 @@ +"""Warning: the pyhawkes and ttpp packages are difficult to install, install +and run at your own risk :) !""" + +# %% Import libraries +import itertools +import time +import numpy as np + +from joblib import Parallel, delayed +import pandas as pd +import torch +from fadin.kernels import DiscreteKernelFiniteSupport +from fadin.solver import FaDIn + + +from pyhawkes.models import ( + DiscreteTimeNetworkHawkesModelGammaMixture, + DiscreteTimeStandardHawkesModel, +) +from pyhawkes.internals.network import StochasticBlockModel + +from ttpp.data import load_dataset +from ttpp.training import train_tripp, train_rnn + +from fadin.utils.utils import projected_grid, projected_grid_marked +from fadin.utils.utils_simu import simu_marked_hawkes_cluster, \ + simu_multi_poisson +from fadin.solver import UNHaP +from fadin.loss_and_gradient import discrete_ll_loss_conv + + +def identity(x, **param): + return x + + +def linear_zero_one(x, **params): + temp = 2 * x + mask = x > 1 + temp[mask] = 0.0 + return temp + + +def reverse_linear_zero_one(x, **params): + temp = 2 - 2 * x + mask = x > 1 + temp[mask] = 0.0 + return temp + + +def truncated_gaussian(x, **params): + rc = DiscreteKernelFiniteSupport( + delta=0.01, n_dim=1, kernel="truncated_gaussian" + ) + mu = params["mu"] + sigma = params["sigma"] + kernel_values = rc.kernel_eval( + [torch.Tensor(mu), torch.Tensor(sigma)], torch.tensor(x) + ) + + return kernel_values.double().numpy() + + +def evaluate_intensity(baseline, alpha, mean, sigma, delta, events_grid): + L = int(1 / delta) + TG = DiscreteKernelFiniteSupport( + delta, n_dim=1, kernel="truncated_gaussian", lower=0, upper=1 + ) + + intens = TG.intensity_eval( + torch.tensor(baseline), + torch.tensor(alpha), + [torch.Tensor(mean), torch.Tensor(sigma)], + events_grid, + torch.linspace(0, 1, L), + ) + return intens + + +def simulate_data(baseline, baseline_noise, alpha, delta, end_time, seed=0): + + n_dim = len(baseline) + + time_kernel = truncated_gaussian + params_time_kernel = dict(mu=mu, sigma=sigma) + + marks_kernel = identity + marks_density = linear_zero_one + time_kernel = truncated_gaussian + + params_marks_density = dict(scale=1) + params_marks_kernel = dict(slope=1.2) + params_time_kernel = dict(mu=mu, sigma=sigma) + + marked_events, _ = simu_marked_hawkes_cluster( + end_time, + baseline, + alpha, + time_kernel, + marks_kernel, + marks_density, + params_time_kernel=params_time_kernel, + params_marks_kernel=params_marks_kernel, + params_marks_density=params_marks_density, + time_kernel_length=None, + marks_kernel_length=None, + upper_bound=None, + random_state=None, + ) + + noisy_events_ = simu_multi_poisson(end_time, baseline_noise) + + random_marks = [ + np.random.rand(noisy_events_[i].shape[0]) / 5.0 for i in range(n_dim) + ] + + noisy_events = [ + np.concatenate( + (noisy_events_[i].reshape(-1, 1), random_marks[i].reshape(-1, 1)), + axis=1 + ) + for i in range(n_dim) + ] + # marked Hawkes concatenated with marked Poisson + mev_cat = [ + np.concatenate((noisy_events[i], marked_events[i]), axis=0) + for i in range(n_dim) + ] + + # marked events concatenated and sorted + mev_cat_sorted = [mev_cat[i][mev_cat[i][:, 0].argsort(0)] + for i in range(n_dim)] + # events concatenated and sorted + ev_cat_sorted = [mev_cat_sorted[i][:, 0] for i in range(n_dim)] + + L = int(1 / delta) + events_grid = projected_grid(ev_cat_sorted, delta, L * end_time + 1) + _, events_grid_clean = projected_grid_marked( + marked_events, delta, L * end_time + 1 + ) + intens = evaluate_intensity(baseline, alpha, mu, sigma, delta, events_grid) + + return ( + ev_cat_sorted, + mev_cat_sorted, + intens, + events_grid, + marked_events, + events_grid_clean, + ) + + +def run_vb(S, events_grid_test, size_grid, dt, seed): + np.random.seed(seed) + init_len = size_grid + start = time.time() + init_model = DiscreteTimeStandardHawkesModel( + K=1, dt=dt, dt_max=1, B=5, alpha=1.0, beta=1.0 + ) + init_model.add_data(S[:init_len, :]) + + init_model.initialize_to_background_rate() + init_model.fit_with_bfgs() + + ########################################################### + # Create a test weak spike-and-slab model + ########################################################### + + test_network = StochasticBlockModel(K=1, C=1) + test_model = DiscreteTimeNetworkHawkesModelGammaMixture( + K=1, dt=dt, dt_max=1, B=5, network=test_network + ) + test_model.add_data(S) + + # Initialize with the standard model parameters + if init_model is not None: + test_model.initialize_with_standard_model(init_model) + + ########################################################### + # Fit the test model with variational Bayesian inference + ########################################################### + # VB coordinate descent + + N_iters = 1000 + vlbs = [] + samples = [] + for itr in range(N_iters): + vlbs.append(test_model.meanfield_coordinate_descent_step()) + + # Resample from variational distribution and plot + test_model.resample_from_mf() + samples.append(test_model.copy_sample()) + + results = dict( + intens_vb=test_model.compute_rate(), + ll=test_model.heldout_log_likelihood( + events_grid_test.T.numpy().astype(np.int64) + ), + time=time.time() - start, + ) + + return results + + +def run_experiment(baseline, baseline_noise, alpha, end_time, delta, seed): + ev_cat_sorted, mev_cat_sorted, true_intens, events_grid, \ + marked_events, events_grid_clean = simulate_data( + baseline, + baseline_noise, + alpha, + delta, + end_time=end_time, + seed=seed + ) + ev_cat_sorted_test, mev_cat_sorted_test, true_intens_test, _, \ + marked_events_test, events_grid_clean_test = \ + simulate_data( + baseline, + baseline_noise, + alpha, + delta, + end_time=end_time, + seed=seed + ) + print(events_grid.sum()) + print(events_grid_clean_test.sum()) + true_LL = discrete_ll_loss_conv(true_intens_test, events_grid_clean_test, + delta, end_time) + + ########################################################## + max_iter = 10_000 + start = time.time() + solver = UNHaP( + n_dim=1, + kernel="truncated_gaussian", + kernel_length=1.0, + delta=delta, + optim="RMSprop", + params_optim={"lr": 1e-3}, + max_iter=max_iter, + batch_rho=100, + init='moment_matching_mean' + ) + + solver.fit(mev_cat_sorted, end_time) + comp_time_mix = time.time() - start + baseline_mix = solver.param_baseline[-10:].mean().item() + alpha_mix = solver.param_alpha[-10:].mean().item() + mu_mix = solver.param_kernel[0][-10:].mean().item() + sigma_mix = solver.param_kernel[1][-10:].mean().item() + print(baseline_mix) + print(alpha_mix) + print(mu_mix) + print(sigma_mix) + mixture_intens = evaluate_intensity( + [baseline_mix], + [[alpha_mix]], + [[mu_mix]], + [[sigma_mix]], + delta, + events_grid_clean_test, + ) + + err_mixture = np.absolute( + mixture_intens.numpy() - true_intens.numpy() + ).mean() + + ll_mixture = discrete_ll_loss_conv( + mixture_intens, events_grid_clean_test, delta, end_time + ) + + results = dict( + err_mixture=err_mixture, + comp_time_mixture=comp_time_mix, + ll_mixture=ll_mixture.item(), + ) + + ########################################################## + start = time.time() + # StocUNHaP + solver = UNHaP( + n_dim=1, + kernel="truncated_gaussian", + kernel_length=1.0, + delta=delta, + optim="RMSprop", + params_optim={"lr": 1e-3}, + max_iter=max_iter, + batch_rho=100, + init='moment_matching_mean', + stoc_classif=True + ) + solver.fit(mev_cat_sorted, end_time) + comp_time_stocunhap = time.time() - start + baseline_stocunhap = solver.param_baseline[-10:].mean().item() + alpha_stocunhap = solver.param_alpha[-10:].mean().item() + mu_stocunhap = solver.param_kernel[0][-10:].mean().item() + sigma_stocunhap = solver.param_kernel[1][-10:].mean().item() + print(baseline_stocunhap) + print(alpha_stocunhap) + print(mu_stocunhap) + print(sigma_stocunhap) + stocunhap_intens = evaluate_intensity( + [baseline_stocunhap], + [[alpha_stocunhap]], + [[mu_stocunhap]], + [[sigma_stocunhap]], + delta, + events_grid_clean_test, + ) + + err_stocunhap = np.absolute( + stocunhap_intens.numpy() - stocunhap_intens.numpy() + ).mean() + + ll_stocunhap = discrete_ll_loss_conv( + stocunhap_intens, events_grid_clean_test, delta, end_time + ) + results["err_stocunhap"] = err_stocunhap + results["comp_time_stocunhap"] = comp_time_stocunhap + results["ll_stocunhap"] = ll_stocunhap.item() + + results["end_time"] = end_time + results["baseline_noise"] = baseline_noise.item() + + ########################################################## + start = time.time() + solver = FaDIn( + n_dim=1, + kernel="truncated_gaussian", + kernel_length=1.0, + delta=delta, + optim="RMSprop", + params_optim={"lr": 1e-3}, + max_iter=max_iter + ) + # Change the mark to one to use FaDIn + mev_cat_sorted_ = [mev_cat_sorted[i].copy() for i in range(1)] + for i in range(1): + mev_cat_sorted_[i][:, 1] = 1.0 + solver.fit(mev_cat_sorted_, end_time) + comp_time_fadin = time.time() - start + baseline_fadin = solver.param_baseline[-10:].mean().item() + alpha_fadin = solver.param_alpha[-10:].mean().item() + mu_fadin = solver.param_kernel[0][-10:].mean().item() + sigma_fadin = solver.param_kernel[1][-10:].mean().item() + print(baseline_fadin) + print(alpha_fadin) + print(mu_fadin) + print(sigma_fadin) + fadin_intens = evaluate_intensity( + [baseline_fadin], + [[alpha_fadin]], + [[mu_fadin]], + [[sigma_fadin]], + delta, + events_grid_clean_test, + ) + + err_fadin = np.absolute(fadin_intens.numpy() - true_intens.numpy()).mean() + + ll_fadin = discrete_ll_loss_conv( + fadin_intens, events_grid_clean_test, delta, end_time + ) + results["err_fadin"] = err_fadin + results["comp_time_fadin"] = comp_time_fadin + results["ll_fadin"] = ll_fadin.item() + results["end_time"] = end_time + results["baseline_noise"] = baseline_noise.item() + + ########################################################## + S = events_grid.T.numpy().astype(np.int64) + size_grid = events_grid.shape[1] + dic_res = run_vb(S, events_grid_clean_test, size_grid, delta, seed) + + results["err_vb"] = np.absolute( + true_intens.numpy() - dic_res["intens_vb"].T + ).mean() + results["ll_vb"] = -dic_res["ll"] / end_time + results["comp_time_vb"] = dic_res["time"] + ########################################################## + # TRIPP + ########################################################## + events_format = [] + mean = 0 + n_dim = 1 + for i in range(n_dim): + events_format.append({"arrival_times": ev_cat_sorted[i]}) + mean += ev_cat_sorted[i].shape[0] + mean /= n_dim + data = {"sequences": [], "t_max": end_time, "mean_number_items": mean} + data["sequences"] = events_format + + dset = load_dataset("None", data) + for i in range(n_dim): + events_format.append({"arrival_times": ev_cat_sorted_test[i]}) + mean += ev_cat_sorted_test[i].shape[0] + mean /= n_dim + data = {"sequences": [], "t_max": end_time, "mean_number_items": mean} + data["sequences"] = events_format + + dset_test = load_dataset("None", data) + d_train, _, _ = dset.train_val_test_split( + train_size=1, val_size=0, test_size=0 + ) + start = time.time() + tripp = train_tripp( + d_train, d_train, + T=end_time, + n_knots=20, + block_size=16, + n_blocks=4, + n_iter=500 + ) + comp_time_tripp = time.time() - start + start = time.time() + rnn = train_rnn( + d_train, d_train, T=end_time, n_knots=20, hidden_size=32, n_iter=1000 + ) + comp_time_rnn = time.time() - start + d_test = torch.utils.data.DataLoader( + dset_test, batch_size=1, shuffle=False + ) + + for x, mask in d_test: + nll_loss_tripp = -(tripp.log_prob(x, mask) / end_time) + + for x, mask in d_test: + nll_loss_rnn = -(rnn.log_prob(x, mask) / end_time) + + results["ll_tripp"] = nll_loss_tripp.detach().item() + results["ll_rnn"] = nll_loss_rnn.detach().item() + + results["comp_time_tripp"] = comp_time_tripp + results["comp_time_rnn"] = comp_time_rnn + + results["true_ll"] = true_LL.item() + return results + + +baseline = np.array([0.1]) +mu = np.array([[0.5]]) +sigma = np.array([[0.1]]) +alpha = np.array([[1.0]]) +delta = 0.01 +end_time_list = [100, 500, 1_000] +baseline_noise_list = [np.array([1.0])] +seeds = np.arange(10) + + +n_jobs = 70 +all_results = Parallel(n_jobs=n_jobs, verbose=10)( + delayed(run_experiment)(baseline, baseline_noise, + alpha, end_time, + delta, seed=seed) + for end_time, baseline_noise, seed in itertools.product( + end_time_list, baseline_noise_list, seeds + ) +) + +# save results +df = pd.DataFrame(all_results) +df.to_csv( + f"results/benchmark_marked_bl_noise{baseline_noise_list}.csv", + index=False +) diff --git a/experiments/unhap/simulated_data/run_benchmark_tpp.py b/experiments/unhap/simulated_data/run_benchmark_tpp.py new file mode 100644 index 0000000..fc069a1 --- /dev/null +++ b/experiments/unhap/simulated_data/run_benchmark_tpp.py @@ -0,0 +1,305 @@ +# %% Import libraries +import itertools +import time +import numpy as np +import pandas as pd +from joblib import Parallel, delayed +import torch + +from fadin.kernels import DiscreteKernelFiniteSupport +from fadin.loss_and_gradient import discrete_ll_loss_conv +from fadin.utils.utils import optimizer_fadin, projected_grid_marked +from fadin.utils.utils import projected_grid +from fadin.utils.utils_simu import simu_marked_hawkes_cluster +from fadin.utils.utils_simu import simu_multi_poisson, custom_density +from fadin.utils.functions import identity, linear_zero_one +from fadin.utils.functions import reverse_linear_zero_one + + +class TPPSelect(): + def __init__(self, n_dim, delta=0.01, k=0.3, n_candidates=10, reg=0): + self.n_dim = n_dim + self.k = k + self.n_candidates = n_candidates + self.reg = reg + self.delta = delta + + self.kernel = DiscreteKernelFiniteSupport( + delta, n_dim=self.n_dim, kernel="truncated_gaussian", + lower=0, kernel_length=1.0 + ) + self.ker_grid = torch.arange(0, 1.00001, delta) + self.L = len(self.ker_grid) + + def intensity(self, events_grid): + """Compute the intensity of the Hawkes process on a grid. + + Parameters: + ----------- + events_grid : torch.Tensor, shape (n_dim, n_grid) + History of the events on a grid as a sparse vector. + + """ + n_dim, n_grid = events_grid.shape + + alpha = self.params[n_dim:n_dim*(n_dim+1)] + kernel_params = self.params[n_dim*(n_dim+1):].reshape(2, 1, 1) + + # Evaluate of the kernel on the grid, shape (n_dim, n_dim, ker_grid) + kernel_val = self.kernel.kernel_eval(kernel_params, self.ker_grid) + kernel_val = kernel_val * alpha.reshape(n_dim, n_dim, 1) + + return torch.conv_transpose1d( + events_grid, kernel_val.transpose(0, 1).float() + )[:, :-self.L + 1] + + def loss(self, events_grid, rho): + n_dim, n_grid = events_grid.shape + + intensity = self.intensity(events_grid) + n_events = events_grid.sum(1) + n_clustered = rho.sum(1) + + baseline = self.params[:n_dim] + integral = torch.norm( + intensity + baseline[:, None], p='fro' + )**2 * self.delta + contrib_cluster = torch.dot( + (events_grid * rho).ravel(), intensity.ravel() + ) + contrib_exo = baseline * (n_events - n_clustered) + return integral - 2 * (contrib_cluster + contrib_exo) + + def inner_min(self, events_grid, rho, max_iter=20): + self.opt = optimizer_fadin( + [self.params], {'lr': 1e-3}, solver='RMSprop' + ) + for i in range(20): + self.opt.zero_grad() + loss = self.loss(events_grid, rho) + loss.backward() + self.opt.step() + return self.loss(events_grid, rho) + + def fit(self, events, end_time): + + n_events = len(events[0]) + k = self.k if isinstance(self.k, int) else int(self.k * n_events) + print(f"Looking for {n_events - k} structured events.") + + n_grid = int(1 / self.delta * end_time) + 1 + events_grid, events_grid_wm = projected_grid_marked( + events, self.delta, n_grid + ) + + # Smart initialization of solver parameters + n_dim = self.n_dim + baseline = n_events / (end_time * (n_dim + 1)) + alpha = 1 / (n_dim + 1) + kernel_params = torch.ones(2) * 0.5 + self.params = torch.tensor([ + baseline, alpha, *kernel_params + ]).requires_grad_(True) + + rho = (events_grid > 0).float() + self.inner_min(events_grid_wm, rho, max_iter=200) + for i in range(k): + print(f"Fitting...{(i+1) / k:6.1%}\r", end='', flush=True) + candidate_events = torch.multinomial( + rho[0], self.n_candidates, replacement=False + ) + F, idx_select = torch.inf, -1 + for idx in candidate_events: + rho[0, idx] = 0 + LL = self.inner_min(events_grid_wm, rho) + rho[0, idx] = 1 + if LL < F: + F, idx_select = LL, idx + rho[0, idx_select] = 0 + + self.param_baseline = self.params[:1].detach() + self.param_alpha = self.params[1:2].detach() + self.param_kernel = self.params[2:].detach().reshape(2, 1) + self.param_rho = rho + print("Fitting... done") + + +def truncated_gaussian(x, **params): + rc = DiscreteKernelFiniteSupport(delta=0.01, n_dim=1, + kernel='truncated_gaussian') + mu = params['mu'] + sigma = params['sigma'] + kernel_values = rc.kernel_eval( + [torch.Tensor(mu), torch.Tensor(sigma)], torch.tensor(x)) + + return kernel_values.double().numpy() + + +def evaluate_intensity(baseline, alpha, mean, sigma, delta, events_grid): + L = int(1 / delta) + TG = DiscreteKernelFiniteSupport( + delta, + n_dim=1, + kernel='truncated_gaussian', + lower=0 + ) + + intens = TG.intensity_eval(torch.tensor(baseline), + torch.tensor([[0.97]]), + [torch.Tensor(mean), + torch.Tensor(sigma)], + events_grid, torch.linspace(0, 1, L)) + return intens + + +def simulate_data(baseline, baseline_noise, alpha, delta, end_time, seed=0): + + n_dim = len(baseline) + + time_kernel = truncated_gaussian + params_time_kernel = dict(mu=mu, sigma=sigma) + + marks_kernel = identity + marks_density = linear_zero_one + time_kernel = truncated_gaussian + + params_marks_density = dict(scale=1) + params_marks_kernel = dict(slope=1.2) + params_time_kernel = dict(mu=mu, sigma=sigma) + + marked_events, _ = simu_marked_hawkes_cluster( + end_time, baseline, alpha, time_kernel, marks_kernel, marks_density, + params_time_kernel=params_time_kernel, + params_marks_kernel=params_marks_kernel, + params_marks_density=params_marks_density, time_kernel_length=None, + marks_kernel_length=None, random_state=seed) + + noisy_events_ = simu_multi_poisson( + end_time, baseline_noise, random_state=seed + ) + np.random.seed(seed) + random_marks = [ + np.random.rand(noisy_events_[i].shape[0]) for i in range(n_dim)] + random_marks = [ + custom_density( + reverse_linear_zero_one, + dict(), + size=noisy_events_[i].shape[0], + kernel_length=1. + ) for i in range(n_dim) + ] + + noisy_events = [ + np.concatenate( + (noisy_events_[i].reshape(-1, 1), random_marks[i].reshape(-1, 1)), + axis=1 + ) for i in range(n_dim) + ] + # marked Hawkes concatenated with marked Poisson + mev_cat = [np.concatenate( + (noisy_events[i], marked_events[i]), axis=0) for i in range(n_dim)] + + # marked events concatenated and sorted + mev_cat_sorted = [ + mev_cat[i][mev_cat[i][:, 0].argsort(0)] for i in range(n_dim) + ] + # events concatenated and sorted + ev_cat_sorted = [mev_cat_sorted[i][:, 0] for i in range(n_dim)] + + L = int(1 / delta) + events_grid = projected_grid(ev_cat_sorted, delta, L * end_time + 1) + _, events_grid_clean = projected_grid_marked( + marked_events, delta, L * end_time + 1 + ) + intens = evaluate_intensity(baseline, alpha, mu, sigma, delta, events_grid) + + return ev_cat_sorted, mev_cat_sorted, intens, events_grid, \ + marked_events, events_grid_clean + + +def run_experiment(baseline, baseline_noise, alpha, end_time, delta, mu, sigma, + seed): + + ev_cat_sorted, mev_cat_sorted, true_intens, events_grid, \ + marked_events, events_grid_clean = simulate_data( + baseline, + baseline_noise, + alpha, + delta=delta, + end_time=end_time, + seed=seed + ) + ev_cat_sorted_test, mev_cat_sorted_test, true_intens_test, _, \ + marked_events_test, events_grid_clean_test = \ + simulate_data( + baseline, + baseline_noise, + alpha, + delta=delta, + end_time=end_time, + seed=seed + ) + + true_LL = discrete_ll_loss_conv(true_intens_test, events_grid_clean_test, + delta, end_time) + + start = time.time() + + solver = TPPSelect(n_dim=1, k=0.2, delta=delta) + + solver.fit(mev_cat_sorted, end_time) + comp_time_mfadin = time.time() - start + baseline_mfadin = solver.param_baseline[-10:].mean().item() + alpha_mfadin = solver.param_alpha[-10:].mean().item() + mu_mfadin = solver.param_kernel[0][-10:].mean().item() + sigma_mfadin = solver.param_kernel[1][-10:].mean().item() + + mfadin_intens = evaluate_intensity( + [baseline_mfadin], [[alpha_mfadin]], [[mu_mfadin]], + [[sigma_mfadin]], delta, events_grid_clean_test) + err_mfadin = np.absolute( + mfadin_intens.numpy() - true_intens.numpy() + ).mean() + + ll_mfadin = discrete_ll_loss_conv(mfadin_intens.clip(1e-5), + events_grid_clean_test, + delta, end_time) + + results = dict(err_mfadin=err_mfadin, + comp_time_mfadin=comp_time_mfadin, + ll_mfadin=ll_mfadin.item(), + true_ll=true_LL.item(), + end_time=end_time, + baseline_noise=baseline_noise.item(), + seed=seed + ) + + return results + + +baseline = np.array([.1]) +mu = np.array([[0.5]]) +sigma = np.array([[0.1]]) +alpha = np.array([[1.0]]) +delta = 0.01 +end_time_list = [100, 500, 1_000] +baseline_noise = [np.array([1.])] +seeds = np.arange(10) + + +n_jobs = 10 # run with 70 to go faster +all_results = Parallel(n_jobs=n_jobs, verbose=10)( + delayed(run_experiment)(baseline, baseline_noise, + alpha, end_time, + delta, mu, sigma, seed=seed) + for end_time, baseline_noise, seed in itertools.product( + end_time_list, baseline_noise, seeds + ) +) + +# save results +df = pd.DataFrame(all_results) +df.to_csv('results/benchmark_rebutall_tpp_select.csv', index=False) + +print('\n Mean \n:', df[df['baseline_noise'] == 1.].groupby('end_time').mean()) +print('\n std \n:', df[df['baseline_noise'] == 1.].groupby('end_time').std()) diff --git a/experiments/unhap/simulated_data/run_infer_noise.py b/experiments/unhap/simulated_data/run_infer_noise.py new file mode 100644 index 0000000..63db47d --- /dev/null +++ b/experiments/unhap/simulated_data/run_infer_noise.py @@ -0,0 +1,1570 @@ +# %% import stuff +import numpy as np +import itertools +from joblib import Parallel, delayed +import pandas as pd +import torch +import time + +from fadin.utils.utils_simu import simu_marked_hawkes_cluster +from fadin.utils.utils_simu import simu_multi_poisson +from fadin.solver import UNHaP +from fadin.kernels import DiscreteKernelFiniteSupport +from fadin.utils.utils import projected_grid_marked, optimizer_fadin +from fadin.init import momentmatching_fadin, momentmatching_unhap +from fadin.utils.compute_constants import get_zG, get_ztzG_approx, get_ztzG +from fadin.loss_and_gradient import discrete_l2_loss_precomputation +from fadin.loss_and_gradient import discrete_ll_loss_conv +from fadin.loss_and_gradient import discrete_l2_loss_conv + + +# %% Define JointFaDIn and JointFaDInDenoising solvers + +def get_zN_joint(events_grid, events_ground_grid, L): + """ + events_grid.shape = n_dim, n_grid + zN.shape = n_dim, n_dim, L + zLN.shape = n_dim, n_dim + """ + n_dim, _ = events_grid.shape + + zN = np.zeros(shape=(n_dim, n_dim, L)) + for i in range(n_dim): + ei = events_grid[i] # Count drived marks + for j in range(n_dim): + ej = events_grid[j] + zN[i, j, 0] = ej @ ei + for tau in range(1, L): + zN[i, j, tau] = ei[tau:] @ ej[:-tau] + + return 2 * zN # to get a density for the distribution of the mark + + +def get_grad_baseline_joint(zG, baseline, alpha, kernel, + delta, end_time, sum_marks, + square_int_marks, n_ground_events): + """Return the gradient of the discrete l2 loss w.r.t. the baseline. + + .. math:: + N_T\\frac{\\partial\\mathcal{L}_G}{\\partial \\mu_{m}} = + 2 T \\mu_m - 2 N_T^m + 2 \\Delta\\sum_{j=1}^{p} \\sum_{\\tau=1}^{L} + \\phi_{mj}^\\Delta[\\tau]\\Phi_{j}(\\tau; G) + + Parameters + ---------- + zG : tensor, shape (n_dim, L) + + baseline : tensor, shape (n_dim,) + Baseline parameter of the intensity of the Hawkes process. + + alpha : tensor, shape (n_dim, n_dim) + Alpha parameter of the intensity of the Hawkes process. + + kernel : tensor, shape (n_dim, n_dim, L) + Kernel values on the discretization. + + delta : float + Step size of the discretization grid. + + sum_marks : tensor, shape (n_dim,) + Sum of marks for each dimension. + + end_time : float + The end time of the Hawkes process. + + n_ground_events : tensor, shape (n_dim,) + Number of events for each dimension. + + Returns + ---------- + grad_baseline: tensor, shape (dim,) + """ + n_dim, _, _ = kernel.shape + + if n_dim > 1: + cst1 = end_time * baseline * square_int_marks + cst2 = 0.5 * n_ground_events.sum() + + dot_kernel = torch.einsum('kju,ju->kj', kernel, zG) + dot_kernel_ = (dot_kernel * alpha).sum(1) * square_int_marks + + grad_baseline = (dot_kernel_ * delta + cst1 - sum_marks) / cst2 + else: + grad_baseline_ = torch.zeros(n_dim) + for k in range(n_dim): + temp = 0 + for j in range(n_dim): + temp += alpha[k, j] * (zG[j] @ kernel[k, j]) + grad_baseline_[k] = delta * temp * square_int_marks[k] + grad_baseline_[k] += end_time * baseline[k] * square_int_marks[k] + grad_baseline_[k] -= 2 * sum_marks[k] # add of joint MarkedFaDIn + grad_baseline = 2 * grad_baseline_ / n_ground_events.sum() + + return grad_baseline + + +def get_grad_alpha_joint(zG, zN, ztzG, baseline, alpha, kernel, delta, + square_int_marks, n_ground_events): + """Return the gradient of the discrete l2 loss w.r.t. alpha. + + .. math:: + N_T\\frac{\\partial\\mathcal{L}_G}{\\partial \\alpha_{ml}} = + 2\\Delta \\mu_m \\sum_{\\tau=1}^{L} \\frac{\\partial + \\phi_{m,l}^\\Delta[\\tau]}{\\partial \\alpha_{m,l}} \\Phi_l(\\tau; G) + + 2 \\Delta \\sum_{k=1}^{p} \\sum_{\\tau=1}^{L} \\sum_{\\tau'=1}^{L} + \\phi_{mk}^\\Delta[\\tau'] \\frac{\\partial \\phi_{m,l}^\\Delta[\\tau]} + {\\partial \\alpha_{m,l}} \\Psi_{l,k}(\\tau, \\tau'; G) + + Parameters + ---------- + zG : tensor, shape (n_dim, L) + + zN : tensor, shape (n_dim, L) + + ztzG : tensor, shape (n_dim, n_dim, L, L) + + baseline : tensor, shape (n_dim,) + Baseline parameter of the intensity of the Hawkes process. + + alpha : tensor, shape (n_dim, n_dim) + Alpha parameter of the intensity of the Hawkes process. + + kernel : tensor, shape (n_dim, n_dim, L) + Kernel values on the discretization. + + delta : float + Step size of the discretization grid. + + n_ground_events : tensor, shape (n_dim,) + Number of events for each dimension. + + Returns + ---------- + grad_alpha : tensor, shape (n_dim, n_dim) + """ + n_dim, _, _ = kernel.shape + + if n_dim > 1: + cst1 = delta * baseline.view(n_dim, 1) + + dot_kernel = torch.einsum('njuv,knu->njkv', ztzG, kernel) + ker_ztzg = torch.einsum('kju,njku->knj', kernel, dot_kernel) + + term1 = torch.einsum('knj,kj->kn', ker_ztzg, alpha) + term2 = torch.einsum('knu,nu->kn', kernel, zG) + term3 = torch.einsum('knu,knu->kn', zN, kernel) + + grad_alpha_ = square_int_marks.unsqueeze(1) * term1 * delta + \ + square_int_marks.unsqueeze(1) * cst1 * term2 - term3 + else: + grad_alpha_ = torch.zeros(n_dim, n_dim) + for k in range(n_dim): + dk = delta * baseline[k] + for n in range(n_dim): + temp = 0 + for j in range(n_dim): + temp += alpha[k, j] * (torch.outer(kernel[k, n], kernel[k, j]) * + ztzG[n, j]).sum() + grad_alpha_[k, n] += delta * temp * square_int_marks[k] + grad_alpha_[k, n] += dk * kernel[k, n] @ zG[n] * square_int_marks[k] + grad_alpha_[k, n] -= zN[k, n] @ kernel[k, n] + + grad_alpha = 2 * grad_alpha_ / n_ground_events.sum() + + return grad_alpha + + +def get_grad_eta_joint(zG, zN, ztzG, baseline, alpha, kernel, + grad_kernel, delta, square_int_marks, n_ground_events): + """Return the gradient of the discrete l2 loss w.r.t. + one kernel parameters. + + .. math:: + N_T\\frac{\\partial\\mathcal{L}_G}{\\partial \\eta{ml}} = + 2\\Delta \\mu_m \\sum_{\\tau=1}^{L} \\frac{\\partial + \\phi_{m,l}^\\Delta[\\tau]}{\\partial \\eta{m,l}} \\Phi_l(\\tau; G) + + 2 \\Delta \\sum_{k=1}^{p} \\sum_{\\tau=1}^{L} \\sum_{\\tau'=1}^{L} + \\phi_{mk}^\\Delta[\\tau'] \\frac{\\partial \\phi_{m,l}^\\Delta[\\tau]} + {\\partial \\eta{m,l}} \\Psi_{l,k}(\\tau, \\tau'; G) + + Parameters + ---------- + zG : tensor, shape (n_dim, L) + + zN : tensor, shape (n_dim, n_dim, L) + + ztzG : tensor, shape (n_dim, n_dim, L, L) + + baseline : tensor, shape (n_dim,) + Baseline parameter of the intensity of the Hawkes process. + + alpha : tensor, shape (n_dim, n_dim) + Alpha parameter of the intensity of the Hawkes process. + + kernel : tensor, shape (n_dim, n_dim, L) + Kernel values on the discretization. + + grad_kernel : list of tensor of shape (n_dim, n_dim, L) + Gradient values on the discretization. + + delta : float + Step size of the discretization grid. + + n_ground_events : tensor, shape (n_dim,) + Number of events for each dimension. + + Returns + ---------- + grad_theta : tensor, shape (n_dim, n_dim) + """ + n_dim, _, L = kernel.shape + grad_theta_ = torch.zeros(n_dim, n_dim) + + if n_dim > 1: + + cst1 = 2 * alpha + cst2 = delta * baseline.view(n_dim, 1) * cst1 + cst3 = torch.einsum('mn,mk->mkn', alpha, alpha) + temp1 = torch.einsum('mnk,nk->mn', grad_kernel, zG) + temp2 = torch.einsum('mnk,mnk->mn', grad_kernel, zN) + temp3 = torch.einsum('nkuv,mnu->nkmv', ztzG, grad_kernel) + temp4 = torch.einsum('mkv,nkmv->mknv', kernel, temp3) + + grad_theta_ = square_int_marks.unsqueeze(1) * cst2 * temp1 - cst1 * temp2 + \ + square_int_marks.unsqueeze(1) * 2 * delta * (cst3 * temp4.sum(3)).sum(1) + else: + for m in range(n_dim): + cst = 2 * delta * baseline[m] + for n in range(n_dim): + grad_theta_[m, n] = cst * alpha[m, n] * (grad_kernel[m, n] @ zG[n]) * \ + square_int_marks[m] + grad_theta_[m, n] -= 2 * alpha[m, n] * (grad_kernel[m, n] @ zN[m, n]) + temp = 0 + for k in range(n_dim): + cst2 = alpha[m, n] * alpha[m, k] + temp_ = 0 + temp_ += 2 * (kernel[m, k].view(1, L) + * (ztzG[n, k] * grad_kernel[m, n].view(L, 1)).sum(0)) + + temp += cst2 * temp_.sum() + + grad_theta_[m, n] += delta * temp * square_int_marks[m] + + grad_theta = grad_theta_ / n_ground_events.sum() + + return grad_theta + + +class JointFaDIn(object): + + get_zN = staticmethod(get_zN_joint) + + def __init__(self, n_dim, kernel, kernel_params_init=None, + baseline_init=None, baseline_mask=None, + alpha_init=None, alpha_mask=None, moment_matching=False, + kernel_length=1, delta=0.01, optim='RMSprop', + params_optim=dict(), max_iter=2000, optimize_kernel=True, + precomputations=True, ztzG_approx=True, device='cpu', + log=False, grad_kernel=None, criterion='l2', tol=10e-5, + random_state=None, optimize_alpha=True): + # Set discretization parameters + self.delta = delta + self.W = kernel_length + self.L = int(self.W / delta) + self.ztzG_approx = ztzG_approx + + # Set optimizer parameters + self.solver = optim + self.max_iter = max_iter + self.log = log + self.tol = tol + + # Set model parameters + self.moment_matching = moment_matching + self.n_dim = n_dim + if baseline_init is None: + self.baseline = torch.rand(self.n_dim) + else: + self.baseline = baseline_init.float() + if baseline_mask is None: + self.baseline_mask = torch.ones([n_dim]) + else: + assert baseline_mask.shape == self.baseline.shape, \ + "Invalid baseline_mask shape, must be (n_dim,)" + self.baseline_mask = baseline_mask + self.baseline = ( + self.baseline * self.baseline_mask + ).requires_grad_(True) + if alpha_init is None: + self.alpha = torch.rand(self.n_dim, self.n_dim) + else: + self.alpha = alpha_init.float() + if alpha_mask is None: + self.alpha_mask = torch.ones([self.n_dim, self.n_dim]) + else: + assert alpha_mask.shape == self.alpha.shape, \ + "Invalid alpha_mask shape, must be (n_dim, n_dim)" + self.alpha_mask = alpha_mask + self.alpha = (self.alpha * self.alpha_mask).requires_grad_(True) + + if kernel_params_init is None: + kernel_params_init = [] + if kernel == 'raised_cosine': + temp = 0.5 * self.W * torch.rand(self.n_dim, self.n_dim) + temp2 = 0.5 * self.W * torch.rand(self.n_dim, self.n_dim) + kernel_params_init.append(temp) + kernel_params_init.append(temp2) + elif kernel == 'truncated_gaussian': + temp = 0.25 * self.W * torch.rand(self.n_dim, self.n_dim) + temp2 = 0.5 * self.W * torch.rand(self.n_dim, self.n_dim) + kernel_params_init.append(temp) + kernel_params_init.append(temp2) + elif kernel == 'truncated_exponential': + kernel_params_init.append(2 * torch.rand(self.n_dim, + self.n_dim)) + else: + raise NotImplementedError( + 'kernel initial parameters of not \ + implemented kernel have to be given' + ) + + self.kernel_params_fixed = kernel_params_init + + self.n_kernel_params = len(kernel_params_init) + + self.kernel_model = DiscreteKernelFiniteSupport(self.delta, self.n_dim, + kernel, + self.W, 0, grad_kernel) + self.kernel = kernel + # Set optimizer + if optimize_alpha: + self.params_intens = [self.baseline, self.alpha] + else: + self.params_intens = [self.baseline] + self.alpha_fixed = self.alpha.detach() + + self.optimize_kernel = optimize_kernel + self.optimize_alpha = optimize_alpha + + if self.optimize_kernel: + for i in range(self.n_kernel_params): + self.params_intens.append( + kernel_params_init[i].float().clip(1e-4).requires_grad_( + True) + ) + + self.precomputations = precomputations + # If the learning rate is not given, fix it to 1e-3 + if 'lr' not in params_optim.keys(): + params_optim['lr'] = 1e-3 + + self.params_solver = params_optim + # self.opt = optimizer(self.params_intens, params_optim, solver=optim) + + if criterion == 'll': + self.precomputations = False + + self.criterion = criterion + # device and seed + if random_state is None: + torch.manual_seed(0) + else: + torch.manual_seed(random_state) + + if torch.cuda.is_available() and device == 'cuda': + self.device = 'cuda' + else: + self.device = 'cpu' + + def fit(self, events, end_time): + """Learn the parameters of the Hawkes processes on a discrete grid. + + Parameters + ---------- + events : list of array of size number of timestamps, + size of the list is dim. + Each array element is a list of size two : + [event timestamp, event mark]. + + end_time : int + The end time of the Hawkes process. + + Returns + ------- + self : object + Fitted parameters. + """ + n_grid = int(1 / self.delta) * end_time + 1 + discretization = torch.linspace(0, self.W, self.L) + events_grid, events_grid_wm = projected_grid_marked( + events, self.delta, n_grid + ) + + sum_marks = events_grid.sum(1) + print('sum of marks per channel:', sum_marks) + n_ground_events = [events[i].shape[0] for i in range(len(events))] + print('number of events is:', n_ground_events) + n_ground_events = torch.tensor(n_ground_events) + self.sum_marks = sum_marks + self.n_ground_events = n_ground_events +# useless in the solver since kernel[i, j, 0] = 0. + max_mark = torch.tensor([ + events[i][:, 1].max() for i in range(len(events))]).float() + min_mark = torch.tensor([ + events[i][:, 1].min() for i in range(len(events))]).float() + + self.square_int_marks = 4 * ((max_mark ** 3) / 3 - (min_mark ** 3) / 3) + + if self.moment_matching: + baseline, alpha, kernel_params_init = momentmatching_fadin( + self, events, n_ground_events, end_time + ) + self.baseline = baseline + self.alpha = alpha + # Set optimizer with moment_matching parameters + if self.optimize_alpha: + self.params_intens = [self.baseline, self.alpha] + else: + self.params_intens = [baseline] + self.alpha_fixed = alpha.detach() + + if self.optimize_kernel: + for i in range(2): # range(self.n_params_kernel) + self.params_intens.append( + kernel_params_init[i].float().clip( + 1e-4).requires_grad_(True) + ) + self.opt = optimizer_fadin( + self.params_intens, + self.params_solver, + solver=self.solver + ) + + #################################################### + # Precomputations + #################################################### + if self.precomputations: + + start = time.time() + zG = get_zG(events_grid.double().numpy(), self.L) + zN = self.get_zN(events_grid.double().numpy(), + events_grid_wm.double().numpy(), + self.L) + + if self.ztzG_approx: + ztzG = get_ztzG_approx(events_grid.double().numpy(), self.L) + else: + ztzG = get_ztzG(events_grid.double().numpy(), self.L) + + zG = torch.tensor(zG).float() + zN = torch.tensor(zN).float() + ztzG = torch.tensor(ztzG).float() + print('precomput:', time.time() - start) + + #################################################### + # save results + #################################################### + self.v_loss = torch.zeros(self.max_iter) + + self.param_baseline = torch.zeros(self.max_iter + 1, self.n_dim) + self.param_alpha = torch.zeros(self.max_iter + 1, + self.n_dim, self.n_dim) + self.param_kernel = torch.zeros(self.n_kernel_params, + self.max_iter + 1, + self.n_dim, self.n_dim) + + self.param_baseline[0] = self.params_intens[0].detach() + + if self.optimize_alpha: + self.param_alpha[0] = self.params_intens[1].detach() + + # If kernel parameters are optimized + if self.optimize_kernel: + for i in range(self.n_kernel_params): + self.param_kernel[i, 0] = \ + self.params_intens[2 + i].detach() + + #################################################### + start = time.time() + for i in range(self.max_iter): + print(f"Fitting model... {i/self.max_iter:6.1%}\r", end='', + flush=True) + + self.opt.zero_grad() + if self.precomputations: + if self.optimize_kernel: + # Update kernel + kernel = self.kernel_model.kernel_eval( + self.params_intens[2:], + discretization + ) + grad_theta = self.kernel_model.grad_eval( + self.params_intens[2:], + discretization + ) + else: + kernel = self.kernel_model.kernel_eval( + self.kernel_params_fixed, + discretization + ) + + if self.log: + self.v_loss[i] = discrete_l2_loss_precomputation( + zG, zN, ztzG, *self.params_intens[:2], kernel, + sum_marks, self.delta, end_time, n_ground_events + ).detach() + + if self.optimize_alpha: + # Update baseline + self.params_intens[0].grad = get_grad_baseline_joint( + zG, + self.params_intens[0], + self.params_intens[1], + kernel, + self.delta, + end_time, + sum_marks, + self.square_int_marks, + n_ground_events) + else: + self.params_intens[0].grad = get_grad_baseline_joint( + zG, + self.params_intens[0], + self.alpha_fixed, + kernel, + self.delta, + end_time, + sum_marks, + self.square_int_marks, + n_ground_events) + if self.optimize_alpha: + # Update alpha + self.params_intens[1].grad = get_grad_alpha_joint( + zG, + zN, + ztzG, + self.params_intens[0], + self.params_intens[1], + kernel, + self.delta, + self.square_int_marks, + n_ground_events) + if self.optimize_kernel: + # Update kernel + for j in range(self.n_kernel_params): + self.params_intens[2 + j].grad = \ + get_grad_eta_joint(zG, + zN, + ztzG, + self.params_intens[0], + self.params_intens[1], + kernel, + grad_theta[j], + self.delta, + self.square_int_marks, + n_ground_events) + + else: + intens = self.kernel_model.intensity_eval( + self.params_intens[0], + self.params_intens[1], + self.params_intens[2:], + events_grid, + discretization + ) + if self.criterion == 'll': + loss = discrete_ll_loss_conv(intens, + events_grid, + self.delta) + else: + loss = discrete_l2_loss_conv(intens, + events_grid, + self.delta) + loss.backward() + + self.opt.step() + + # Save parameters + self.params_intens[0].data = self.params_intens[0].data.clip(0) * \ + self.baseline_mask + if self.optimize_alpha: + self.params_intens[1].data = self.params_intens[1].data.clip(0) * \ + self.alpha_mask + self.param_alpha[i + 1] = self.params_intens[1].detach() + + self.param_baseline[i + 1] = self.params_intens[0].detach() + + # If kernel parameters are optimized + if self.optimize_kernel: + for j in range(self.n_kernel_params): + self.params_intens[2 + j].data = \ + self.params_intens[2 + j].data.clip(1e-3) + self.param_kernel[j, i + 1] = \ + self.params_intens[2 + j].detach() + + # Early stopping + if i % 100 == 0: + error_b = torch.abs(self.param_baseline[i + 1] - + self.param_baseline[i]).max() + error_al = torch.abs(self.param_alpha[i + 1] - + self.param_alpha[i]).max() + error_k = torch.abs(self.param_kernel[0, i + 1] - + self.param_kernel[0, i]).max() + + if error_b < self.tol and error_al < self.tol and error_k < self.tol: + print('early stopping at iteration:', i) + self.param_baseline = self.param_baseline[:i + 1] + self.param_alpha = self.param_alpha[:i + 1] + for j in range(self.n_kernel_params): + self.param_kernel[j] = self.param_kernel[j, i + 1] + break + print('iterations in ', time.time() - start) + + return self + + +def get_grad_baseline_joint_noise2(zG, baseline, baseline_noise, alpha, kernel, + delta, end_time, sum_marks, + square_int_marks, n_ground_events): + """Return the gradient of the discrete l2 loss w.r.t. the baseline. + + .. math:: + N_T\\frac{\\partial\\mathcal{L}_G}{\\partial \\mu_{m}} = + 2 T \\mu_m - 2 N_T^m + 2 \\Delta\\sum_{j=1}^{p} \\sum_{\\tau=1}^{L} + \\phi_{mj}^\\Delta[\\tau]\\Phi_{j}(\\tau; G) + + Parameters + ---------- + zG : tensor, shape (n_dim, L) + + baseline : tensor, shape (n_dim,) + Baseline parameter of the intensity of the Hawkes process. + + alpha : tensor, shape (n_dim, n_dim) + Alpha parameter of the intensity of the Hawkes process. + + kernel : tensor, shape (n_dim, n_dim, L) + Kernel values on the discretization. + + delta : float + Step size of the discretization grid. + + sum_marks : tensor, shape (n_dim,) + Sum of marks for each dimension. + + end_time : float + The end time of the Hawkes process. + + n_ground_events : tensor, shape (n_dim,) + Number of events for each dimension. + + Returns + ---------- + grad_baseline: tensor, shape (dim,) + """ + n_dim, _, _ = kernel.shape + + if n_dim > 1: + cst1 = end_time * baseline * square_int_marks + cst2 = 0.5 * n_ground_events.sum() + cst3 = end_time * baseline_noise + + dot_kernel = torch.einsum('kju,ju->kj', kernel, zG) + dot_kernel_ = (dot_kernel * alpha).sum(1) * square_int_marks + + grad_baseline = (dot_kernel_ * delta + cst1 - sum_marks + cst3) / cst2 + else: + grad_baseline_ = torch.zeros(n_dim) + for k in range(n_dim): + temp = 0 + for j in range(n_dim): + temp += alpha[k, j] * (zG[j] @ kernel[k, j]) + grad_baseline_[k] = delta * temp * square_int_marks[k] + grad_baseline_[k] += end_time * baseline[k] * square_int_marks[k] + grad_baseline_[k] -= 2 * sum_marks[k] # add of joint MarkedFaDIn + grad_baseline_[k] += end_time * baseline_noise[k] + grad_baseline = 2 * grad_baseline_ / n_ground_events.sum() + + return grad_baseline + + +def get_grad_baseline_noise_joint_noise2(zG, baseline, baseline_noise, alpha, kernel, + delta, end_time, sum_marks, square_int_marks, + n_ground_events): + """Return the gradient of the discrete l2 loss w.r.t. the baseline. + + .. math:: + N_T\\frac{\\partial\\mathcal{L}_G}{\\partial \\mu_{m}} = + 2 T \\mu_m - 2 N_T^m + 2 \\Delta\\sum_{j=1}^{p} \\sum_{\\tau=1}^{L} + \\phi_{mj}^\\Delta[\\tau]\\Phi_{j}(\\tau; G) + + Parameters + ---------- + zG : tensor, shape (n_dim, L) + + baseline : tensor, shape (n_dim,) + Baseline parameter of the intensity of the Hawkes process. + + alpha : tensor, shape (n_dim, n_dim) + Alpha parameter of the intensity of the Hawkes process. + + kernel : tensor, shape (n_dim, n_dim, L) + Kernel values on the discretization. + + delta : float + Step size of the discretization grid. + + sum_marks : tensor, shape (n_dim,) + Sum of marks for each dimension. + + end_time : float + The end time of the Hawkes process. + + n_ground_events : tensor, shape (n_dim,) + Number of events for each dimension. + + Returns + ---------- + grad_baseline: tensor, shape (dim,) + """ + n_dim, _, _ = kernel.shape + + if n_dim > 1: + cst1 = end_time * baseline_noise + cst2 = 0.5 * n_ground_events.sum() + cst3 = end_time * baseline + + dot_kernel = torch.einsum('kju,ju->kj', kernel, zG) + dot_kernel_ = (dot_kernel * alpha).sum(1) + + grad_baseline_noise = ( + dot_kernel_ * delta + cst1 - n_ground_events + cst3) / cst2 + else: + grad_baseline_noise_ = torch.zeros(n_dim) + for k in range(n_dim): + temp = 0 + for j in range(n_dim): + temp += alpha[k, j] * (zG[j] @ kernel[k, j]) + grad_baseline_noise_[k] = delta * temp + grad_baseline_noise_[k] += end_time * baseline_noise[k] + grad_baseline_noise_[k] -= n_ground_events[k] + grad_baseline_noise_[k] += end_time * baseline[k] + grad_baseline_noise = 2 * grad_baseline_noise_ / n_ground_events.sum() + + return grad_baseline_noise + + +def get_grad_baseline_noise_joint_noise2(zG, baseline, baseline_noise, alpha, kernel, + delta, end_time, sum_marks, square_int_marks, + n_ground_events): + """Return the gradient of the discrete l2 loss w.r.t. the baseline. + + .. math:: + N_T\\frac{\\partial\\mathcal{L}_G}{\\partial \\mu_{m}} = + 2 T \\mu_m - 2 N_T^m + 2 \\Delta\\sum_{j=1}^{p} \\sum_{\\tau=1}^{L} + \\phi_{mj}^\\Delta[\\tau]\\Phi_{j}(\\tau; G) + + Parameters + ---------- + zG : tensor, shape (n_dim, L) + + baseline : tensor, shape (n_dim,) + Baseline parameter of the intensity of the Hawkes process. + + alpha : tensor, shape (n_dim, n_dim) + Alpha parameter of the intensity of the Hawkes process. + + kernel : tensor, shape (n_dim, n_dim, L) + Kernel values on the discretization. + + delta : float + Step size of the discretization grid. + + sum_marks : tensor, shape (n_dim,) + Sum of marks for each dimension. + + end_time : float + The end time of the Hawkes process. + + n_ground_events : tensor, shape (n_dim,) + Number of events for each dimension. + + Returns + ---------- + grad_baseline: tensor, shape (dim,) + """ + n_dim, _, _ = kernel.shape + + if n_dim > 1: + cst1 = end_time * baseline_noise + cst2 = 0.5 * n_ground_events.sum() + cst3 = end_time * baseline + + dot_kernel = torch.einsum('kju,ju->kj', kernel, zG) + dot_kernel_ = (dot_kernel * alpha).sum(1) + + grad_baseline_noise = ( + dot_kernel_ * delta + cst1 - n_ground_events + cst3) / cst2 + else: + grad_baseline_noise_ = torch.zeros(n_dim) + for k in range(n_dim): + temp = 0 + for j in range(n_dim): + temp += alpha[k, j] * (zG[j] @ kernel[k, j]) + grad_baseline_noise_[k] = delta * temp + grad_baseline_noise_[k] += end_time * baseline_noise[k] + grad_baseline_noise_[k] -= n_ground_events[k] + grad_baseline_noise_[k] += end_time * baseline[k] + grad_baseline_noise = 2 * grad_baseline_noise_ / n_ground_events.sum() + + return grad_baseline_noise + + +def get_grad_alpha_joint_noise2(zG, zN, ztzG, baseline, baseline_noise, alpha, kernel, + delta, square_int_marks, n_ground_events): + """Return the gradient of the discrete l2 loss w.r.t. alpha. + + .. math:: + N_T\\frac{\\partial\\mathcal{L}_G}{\\partial \\alpha_{ml}} = + 2\\Delta \\mu_m \\sum_{\\tau=1}^{L} \\frac{\\partial + \\phi_{m,l}^\\Delta[\\tau]}{\\partial \\alpha_{m,l}} \\Phi_l(\\tau; G) + + 2 \\Delta \\sum_{k=1}^{p} \\sum_{\\tau=1}^{L} \\sum_{\\tau'=1}^{L} + \\phi_{mk}^\\Delta[\\tau'] \\frac{\\partial \\phi_{m,l}^\\Delta[\\tau]} + {\\partial \\alpha_{m,l}} \\Psi_{l,k}(\\tau, \\tau'; G) + + Parameters + ---------- + zG : tensor, shape (n_dim, L) + + zN : tensor, shape (n_dim, L) + + ztzG : tensor, shape (n_dim, n_dim, L, L) + + baseline : tensor, shape (n_dim,) + Baseline parameter of the intensity of the Hawkes process. + + alpha : tensor, shape (n_dim, n_dim) + Alpha parameter of the intensity of the Hawkes process. + + kernel : tensor, shape (n_dim, n_dim, L) + Kernel values on the discretization. + + delta : float + Step size of the discretization grid. + + n_ground_events : tensor, shape (n_dim,) + Number of events for each dimension. + + Returns + ---------- + grad_alpha : tensor, shape (n_dim, n_dim) + """ + n_dim, _, _ = kernel.shape + + if n_dim > 1: + cst1 = delta * (baseline.view(n_dim, 1) * square_int_marks.unsqueeze(1) + + baseline_noise.view(n_dim, 1)) + + dot_kernel = torch.einsum('njuv,knu->njkv', ztzG, kernel) + ker_ztzg = torch.einsum('kju,njku->knj', kernel, dot_kernel) + + term1 = torch.einsum('knj,kj->kn', ker_ztzg, alpha) + term2 = torch.einsum('knu,nu->kn', kernel, zG) + term3 = torch.einsum('knu,knu->kn', zN, kernel) + + grad_alpha_ = square_int_marks.unsqueeze(1) * term1 * delta + \ + cst1 * term2 - term3 + else: + grad_alpha_ = torch.zeros(n_dim, n_dim) + for k in range(n_dim): + dk = delta * (baseline[k] * square_int_marks[k] + baseline_noise[k]) + for n in range(n_dim): + temp = 0 + for j in range(n_dim): + temp += alpha[k, j] * (torch.outer(kernel[k, n], kernel[k, j]) * + ztzG[n, j]).sum() + grad_alpha_[k, n] += delta * temp * (square_int_marks[k]) + grad_alpha_[k, n] += dk * kernel[k, n] @ zG[n] + grad_alpha_[k, n] -= zN[k, n] @ kernel[k, n] + + grad_alpha = 2 * grad_alpha_ / n_ground_events.sum() + + return grad_alpha + + +def get_grad_eta_joint_noise2(zG, zN, ztzG, baseline, baseline_noise, alpha, kernel, + grad_kernel, delta, square_int_marks, n_ground_events): + """Return the gradient of the discrete l2 loss w.r.t. + one kernel parameters. + + .. math:: + N_T\\frac{\\partial\\mathcal{L}_G}{\\partial \\eta{ml}} = + 2\\Delta \\mu_m \\sum_{\\tau=1}^{L} \\frac{\\partial + \\phi_{m,l}^\\Delta[\\tau]}{\\partial \\eta{m,l}} \\Phi_l(\\tau; G) + + 2 \\Delta \\sum_{k=1}^{p} \\sum_{\\tau=1}^{L} \\sum_{\\tau'=1}^{L} + \\phi_{mk}^\\Delta[\\tau'] \\frac{\\partial \\phi_{m,l}^\\Delta[\\tau]} + {\\partial \\eta{m,l}} \\Psi_{l,k}(\\tau, \\tau'; G) + + Parameters + ---------- + zG : tensor, shape (n_dim, L) + + zN : tensor, shape (n_dim, n_dim, L) + + ztzG : tensor, shape (n_dim, n_dim, L, L) + + baseline : tensor, shape (n_dim,) + Baseline parameter of the intensity of the Hawkes process. + + alpha : tensor, shape (n_dim, n_dim) + Alpha parameter of the intensity of the Hawkes process. + + kernel : tensor, shape (n_dim, n_dim, L) + Kernel values on the discretization. + + grad_kernel : list of tensor of shape (n_dim, n_dim, L) + Gradient values on the discretization. + + delta : float + Step size of the discretization grid. + + n_ground_events : tensor, shape (n_dim,) + Number of events for each dimension. + + Returns + ---------- + grad_theta : tensor, shape (n_dim, n_dim) + """ + n_dim, _, L = kernel.shape + grad_theta_ = torch.zeros(n_dim, n_dim) + + if n_dim > 1: + + cst1 = 2 * alpha + cst2 = delta * (baseline.view(n_dim, 1) * square_int_marks.unsqueeze(1) + + baseline_noise.view(n_dim, 1)) * cst1 + cst3 = torch.einsum('mn,mk->mkn', alpha, alpha) + temp1 = torch.einsum('mnk,nk->mn', grad_kernel, zG) + temp2 = torch.einsum('mnk,mnk->mn', grad_kernel, zN) + temp3 = torch.einsum('nkuv,mnu->nkmv', ztzG, grad_kernel) + temp4 = torch.einsum('mkv,nkmv->mknv', kernel, temp3) + + grad_theta_ = cst2 * temp1 - cst1 * temp2 + ( + square_int_marks.unsqueeze(1) + 1) * 2 * delta * ( + cst3 * temp4.sum(3)).sum(1) + else: + for m in range(n_dim): + cst = 2 * delta * (baseline[m] * square_int_marks[m] + baseline_noise[m]) + for n in range(n_dim): + grad_theta_[m, n] = cst * alpha[m, n] * (grad_kernel[m, n] @ zG[n]) + + grad_theta_[m, n] -= 2 * alpha[m, n] * (grad_kernel[m, n] @ zN[m, n]) + temp = 0 + for k in range(n_dim): + cst2 = alpha[m, n] * alpha[m, k] + temp_ = 0 + temp_ += 2 * (kernel[m, k].view(1, L) + * (ztzG[n, k] * grad_kernel[m, n].view(L, 1)).sum(0)) + + temp += cst2 * temp_.sum() + + grad_theta_[m, n] += delta * temp * (square_int_marks[m]) + + grad_theta = grad_theta_ / n_ground_events.sum() + + return grad_theta + + +class JointFaDInDenoising(object): + """Joint inference of the Hawkes process and the noise, without unmixing. + + This class is used for ablation study and should not be used in practice. + """ + + get_zN = staticmethod(get_zN_joint) + + def __init__(self, n_dim, kernel, kernel_params_init=None, + baseline_init=None, baseline_mask=None, + alpha_init=None, alpha_mask=None, moment_matching=True, + kernel_length=1, delta=0.01, optim='RMSprop', + params_optim=dict(), max_iter=2000, optimize_kernel=True, + precomputations=True, ztzG_approx=True, device='cpu', + log=False, grad_kernel=None, criterion='l2', tol=10e-5, + random_state=None, optimize_alpha=True): + # Set discretization parameters + self.delta = delta + self.W = kernel_length + self.L = int(self.W / delta) + self.ztzG_approx = ztzG_approx + + # Set optimizer parameters + self.solver = optim + self.max_iter = max_iter + self.log = log + self.tol = tol + + # Set model parameters + self.moment_matching = moment_matching + self.n_dim = n_dim + if baseline_init is None: + self.baseline = torch.rand(self.n_dim) + else: + self.baseline = baseline_init.float() + if baseline_mask is None: + self.baseline_mask = torch.ones([n_dim]) + else: + assert baseline_mask.shape == self.baseline.shape, \ + "Invalid baseline_mask shape, must be (n_dim,)" + self.baseline_mask = baseline_mask + self.baseline = (self.baseline * self.baseline_mask).requires_grad_(True) + + self.baseline_noise = torch.rand(self.n_dim).requires_grad_(True) + + if alpha_init is None: + self.alpha = torch.rand(self.n_dim, self.n_dim) + else: + self.alpha = alpha_init.float() + if alpha_mask is None: + self.alpha_mask = torch.ones([self.n_dim, self.n_dim]) + else: + assert alpha_mask.shape == self.alpha.shape, \ + "Invalid alpha_mask shape, must be (n_dim, n_dim)" + self.alpha_mask = alpha_mask + self.alpha = (self.alpha * self.alpha_mask).requires_grad_(True) + + if kernel_params_init is None: + kernel_params_init = [] + if kernel == 'raised_cosine': + temp = 0.5 * self.W * torch.rand(self.n_dim, self.n_dim) + temp2 = 0.5 * self.W * torch.rand(self.n_dim, self.n_dim) + kernel_params_init.append(temp) + kernel_params_init.append(temp2) + elif kernel == 'truncated_gaussian': + temp = 0.25 * self.W * torch.rand(self.n_dim, self.n_dim) + temp2 = 0.5 * self.W * torch.rand(self.n_dim, self.n_dim) + kernel_params_init.append(temp) + kernel_params_init.append(temp2) + elif kernel == 'truncated_exponential': + kernel_params_init.append(2 * torch.rand(self.n_dim, + self.n_dim)) + else: + raise NotImplementedError( + 'kernel initial parameters of not \ + implemented kernel have to be given' + ) + + self.kernel_params_fixed = kernel_params_init + + self.n_kernel_params = len(kernel_params_init) + + self.kernel_model = DiscreteKernelFiniteSupport(self.delta, self.n_dim, + kernel, + self.W, 0, grad_kernel) + self.kernel = kernel + # Set optimizer + if optimize_alpha: + self.params_intens = [ + self.baseline, + self.baseline_noise, + self.alpha + ] + else: + self.params_intens = [self.baseline, self.baseline_noise] + self.alpha_fixed = self.alpha.detach() + + self.optimize_kernel = optimize_kernel + self.optimize_alpha = optimize_alpha + + if self.optimize_kernel: + for i in range(self.n_kernel_params): + self.params_intens.append( + kernel_params_init[i].float().clip(1e-4).requires_grad_( + True + ) + ) + + self.precomputations = precomputations + # If the learning rate is not given, fix it to 1e-3 + if 'lr' not in params_optim.keys(): + params_optim['lr'] = 1e-3 + + self.params_solver = params_optim + self.opt = optimizer_fadin( + self.params_intens, + params_optim, + solver=optim + ) + + if criterion == 'll': + self.precomputations = False + + self.criterion = criterion + # device and seed + if random_state is None: + torch.manual_seed(0) + else: + torch.manual_seed(random_state) + + if torch.cuda.is_available() and device == 'cuda': + self.device = 'cuda' + else: + self.device = 'cpu' + + def fit(self, events, end_time): + """Learn the parameters of the Hawkes processes on a discrete grid. + + Parameters + ---------- + events : list of array of size number of timestamps, + size of the list is dim. + Each array element is a list of size two : + [event timestamp, event mark]. + + end_time : int + The end time of the Hawkes process. + + Returns + ------- + self : object + Fitted parameters. + """ + n_grid = int(1 / self.delta) * end_time + 1 + discretization = torch.linspace(0, self.W, self.L) + events_grid, events_grid_wm = projected_grid_marked( + events, self.delta, n_grid + ) + + sum_marks = events_grid.sum(1) + print('sum of marks per channel:', sum_marks) + n_ground_events = [events[i].shape[0] for i in range(len(events))] + print('number of events is:', n_ground_events) + n_ground_events = torch.tensor(n_ground_events) + self.sum_marks = sum_marks + self.n_ground_events = n_ground_events + + self.square_int_marks = torch.tensor([4/3 for _ in range(self.n_dim)]) + + if self.moment_matching: + # Smart initialization of solver parameters + baseline, bl_noise, alpha, kernel_params_init = momentmatching_unhap( + self, events, events_grid, n_ground_events, end_time + ) + # Initial baseline and alpha divided by n_dim+2 instead of n_dim+1 + self.baseline = baseline + self.bl_noise = bl_noise + self.alpha = alpha + # Set optimizer with moment_matching parameters + if self.optimize_alpha: + self.params_intens = [self.baseline, self.baseline_noise, + self.alpha] + else: + self.params_intens = [self.baseline, self.baseline_noise] + self.alpha_fixed = self.alpha.detach() + + if self.optimize_kernel: + for i in range(2): # range(self.n_params_kernel) + self.params_intens.append( + kernel_params_init[i].float().clip(1e-4).requires_grad_(True) + ) + self.opt = optimizer_fadin( + self.params_intens, + self.params_solver, + solver=self.solver + ) + + #################################################### + # Precomputations + #################################################### + if self.precomputations: + + start = time.time() + zG = get_zG(events_grid.double().numpy(), self.L) + zN = self.get_zN(events_grid.double().numpy(), + events_grid_wm.double().numpy(), + self.L) + + if self.ztzG_approx: + ztzG = get_ztzG_approx(events_grid.double().numpy(), self.L) + else: + ztzG = get_ztzG(events_grid.double().numpy(), self.L) + + zG = torch.tensor(zG).float() + zN = torch.tensor(zN).float() + ztzG = torch.tensor(ztzG).float() + print('precomput:', time.time() - start) + + #################################################### + # save results + #################################################### + self.v_loss = torch.zeros(self.max_iter) + + self.param_baseline = torch.zeros(self.max_iter + 1, self.n_dim) + self.param_baseline_noise = torch.zeros(self.max_iter + 1, self.n_dim) + self.param_alpha = torch.zeros(self.max_iter + 1, + self.n_dim, self.n_dim) + self.param_kernel = torch.zeros(self.n_kernel_params, + self.max_iter + 1, + self.n_dim, self.n_dim) + + self.param_baseline[0] = self.params_intens[0].detach() + self.param_baseline_noise[0] = self.params_intens[1].detach() + + if self.optimize_alpha: + self.param_alpha[0] = self.params_intens[2].detach() + + # If kernel parameters are optimized + if self.optimize_kernel: + for i in range(self.n_kernel_params): + self.param_kernel[i, 0] = self.params_intens[3 + i].detach() + + #################################################### + start = time.time() + for i in range(self.max_iter): + print(f"Fitting model... {i/self.max_iter:6.1%}\r", end='', + flush=True) + + self.opt.zero_grad() + if self.precomputations: + if self.optimize_kernel: + # Update kernel + kernel = self.kernel_model.kernel_eval( + self.params_intens[3:], + discretization + ) + grad_theta = self.kernel_model.grad_eval( + self.params_intens[3:], + discretization + ) + else: + kernel = self.kernel_model.kernel_eval( + self.kernel_params_fixed, + discretization + ) + + if self.log: + # Not done + self.v_loss[i] = \ + discrete_l2_loss_precomputation( + zG, zN, ztzG, + self.params_intens[0], + self.params_intens[1], + kernel, sum_marks, + self.delta, + end_time, + n_ground_events + ).detach() + + self.params_intens[0].grad = get_grad_baseline_joint_noise2( + zG, + self.params_intens[0], + self.params_intens[1], + self.params_intens[2], + kernel, + self.delta, + end_time, + sum_marks, + self.square_int_marks, + n_ground_events) + + self.params_intens[1].grad = get_grad_baseline_noise_joint_noise2( + zG, + self.params_intens[0], + self.params_intens[1], + self.params_intens[2], + kernel, + self.delta, + end_time, + sum_marks, + self.square_int_marks, + n_ground_events) + + if self.optimize_alpha: + # Update alpha + self.params_intens[2].grad = get_grad_alpha_joint_noise2( + zG, + zN, + ztzG, + self.params_intens[0], + self.params_intens[1], + self.params_intens[2], + kernel, + self.delta, + self.square_int_marks, + n_ground_events) + if self.optimize_kernel: + # Update kernel + for j in range(self.n_kernel_params): + self.params_intens[3 + j].grad = \ + get_grad_eta_joint_noise2(zG, + zN, + ztzG, + self.params_intens[0], + self.params_intens[1], + self.params_intens[2], + kernel, + grad_theta[j], + self.delta, + self.square_int_marks, + n_ground_events) + + else: + # Not done + intens = self.kernel_model.intensity_eval( + self.params_intens[0], + self.params_intens[1], + self.params_intens[2:], + events_grid, + discretization + ) + + # Not done + if self.criterion == 'll': + loss = discrete_ll_loss_conv(intens, + events_grid, + self.delta) + else: + loss = discrete_l2_loss_conv(intens, + events_grid, + self.delta) + loss.backward() + + self.opt.step() + + # Save parameters + self.params_intens[0].data = self.params_intens[0].data.clip(0) * \ + self.baseline_mask + self.params_intens[1].data = self.params_intens[1].data.clip(0) + if self.optimize_alpha: + self.params_intens[2].data = self.params_intens[2].data.clip(0) * \ + self.alpha_mask + self.param_alpha[i + 1] = self.params_intens[2].detach() + + self.param_baseline[i + 1] = self.params_intens[0].detach() + self.param_baseline_noise[i + 1] = self.params_intens[1].detach() + + # If kernel parameters are optimized + if self.optimize_kernel: + for j in range(self.n_kernel_params): + self.params_intens[3 + j].data = \ + self.params_intens[3 + j].data.clip(1e-3) + self.param_kernel[j, i + 1] = self.params_intens[3 + j].detach() + + # Early stopping + if i % 100 == 0: + error_b = torch.abs(self.param_baseline[i + 1] - + self.param_baseline[i]).max() + error_al = torch.abs(self.param_alpha[i + 1] - + self.param_alpha[i]).max() + error_k = torch.abs(self.param_kernel[0, i + 1] - + self.param_kernel[0, i]).max() + + if error_b < self.tol and error_al < self.tol and error_k < self.tol: + print('early stopping at iteration:', i) + self.param_baseline = self.param_baseline[:i + 1] + self.param_baseline_noise = self.param_baseline_noise[:i + 1] + self.param_alpha = self.param_alpha[:i + 1] + for j in range(self.n_kernel_params): + self.param_kernel[j] = self.param_kernel[j, i + 1] + break + print('iterations in ', time.time() - start) + + return self + + +# %% Define utilitary functions +def simulate_data(baseline, baseline_noise, alpha, end_time, seed=0): + n_dim = len(baseline) + + def identity(x, **param): + return x + + def linear_zero_one(x, **params): + temp = 2 * x + mask = x > 1 + temp[mask] = 0. + return temp + + def truncated_gaussian(x, **params): + rc = DiscreteKernelFiniteSupport(delta=0.01, n_dim=1, + kernel='truncated_gaussian') + mu = params['mu'] + sigma = params['sigma'] + kernel_values = rc.kernel_eval( + [torch.Tensor(mu), torch.Tensor(sigma)], torch.tensor(x)) + + return kernel_values.double().numpy() + + marks_kernel = identity + marks_density = linear_zero_one + time_kernel = truncated_gaussian + + params_marks_density = dict() + # params_marks_density = dict(scale=1) + params_marks_kernel = dict(slope=1.2) + params_time_kernel = dict(mu=mu, sigma=sigma) + + marked_events, _ = simu_marked_hawkes_cluster( + end_time, baseline, alpha, time_kernel, marks_kernel, marks_density, + params_marks_kernel=params_marks_kernel, + params_marks_density=params_marks_density, + time_kernel_length=None, marks_kernel_length=None, + params_time_kernel=params_time_kernel, random_state=seed) + + noisy_events_ = simu_multi_poisson(end_time, [baseline_noise]) + + random_marks = [ + np.random.rand(noisy_events_[i].shape[0]) for i in range(n_dim)] + noisy_events = [ + np.concatenate((noisy_events_[i].reshape(-1, 1), + random_marks[i].reshape(-1, 1)), axis=1) for i in range(n_dim)] + + events = [ + np.concatenate( + (noisy_events[i], marked_events[i]), axis=0) for i in range(n_dim)] + + events_cat = [events[i][events[i][:, 0].argsort()] for i in range(n_dim)] + # put the mark to one to test the impact of the marks + # events_cat[0][:, 1] = 1. + + return events_cat + + +def run_experiment(baseline, baseline_noise, alpha, end_time, delta, seed): + + events = simulate_data(baseline, baseline_noise, alpha, + end_time=end_time, seed=seed) + + max_iter = 10000 + start = time.time() + solver = UNHaP( + n_dim=1, + kernel="truncated_gaussian", + kernel_length=1., + delta=delta, optim="RMSprop", + params_optim={'lr': 1e-3}, + max_iter=max_iter, + batch_rho=200 + ) + + solver.fit(events, end_time) + + comp_time_mix = time.time() - start + results = dict(param_baseline_mix=solver.param_baseline_[-10:].mean().item(), + param_baseline_noise_mix=solver.param_baseline_noise_[ + -10:].mean().item(), + param_alpha_mix=solver.param_alpha_[-10:].mean().item(), + param_mu_mix=solver.param_kernel_[0][-10:].mean().item(), + param_sigma_mix=solver.param_kernel_[1][-10:].mean().item()) + + start = time.time() + solver = JointFaDIn( + n_dim=1, + kernel="truncated_gaussian", + kernel_length=1., + delta=delta, optim="RMSprop", + params_optim={'lr': 1e-3}, + max_iter=max_iter, criterion='l2', + optimize_kernel=True, + optimize_alpha=True + ) + + solver.fit(events, end_time) + + comp_time_joint = time.time() - start + results["param_baseline_joint"] = solver.param_baseline[-10:].mean().item() + results["param_alpha_joint"] = solver.param_alpha[-10:].mean().item() + results["param_mu_joint"] = solver.param_kernel[0][-10:].mean().item() + results["param_sigma_joint"] = solver.param_kernel[1][-10:].mean().item() + + start = time.time() + solver = JointFaDInDenoising( + n_dim=1, + kernel="truncated_gaussian", + kernel_length=1., + delta=delta, optim="RMSprop", + params_optim={'lr': 1e-3}, + max_iter=max_iter, criterion='l2', + optimize_kernel=True, + optimize_alpha=True + ) + + solver.fit(events, end_time) + comp_time_denoise = time.time() - start + + results["param_baseline_denoise"] = solver.param_baseline[-10:].mean().item() + results["param_alpha_denoise"] = solver.param_alpha[-10:].mean().item() + results["param_mu_denoise"] = solver.param_kernel[0][-10:].mean().item() + results["param_sigma_denoise"] = solver.param_kernel[1][-10:].mean().item() + + results["seed"] = seed + results["end_time"] = end_time + results["delta"] = delta + results["noise"] = baseline_noise.item() + results["comp_time_mix"] = comp_time_mix + results["comp_time_joint"] = comp_time_joint + results["comp_time_denoise"] = comp_time_denoise + + return results + + +# Parameters +setting = 'high' # high +baseline = np.array([.8]) +mu = np.array([[0.5]]) +sigma = np.array([[0.1]]) + +# low structure +if setting == 'low': + alpha = np.array([[1.3]]) +else: + # high structure + alpha = np.array([[1.45]]) + + +delta = 0.01 +end_time_list = [100, 10_00, 100_00] +baseline_noise_list = np.linspace(0.1, 1.5, 10) +seeds = np.arange(100) + +n_jobs = 10 +all_results = Parallel(n_jobs=n_jobs, verbose=10)( + delayed(run_experiment)(baseline, baseline_noise, + alpha, end_time, + delta, seed=seed) + for end_time, baseline_noise, seed in itertools.product( + end_time_list, baseline_noise_list, seeds + ) +) + +# save results +df = pd.DataFrame(all_results) +true_param = {'baseline': baseline.item(), 'alpha': alpha.item(), + 'mu': mu.item(), 'sigma': sigma.item()} +for param, value in true_param.items(): + df[param] = value + + +def compute_norm2_error_joint(s): + return np.sqrt(np.array([(s[param] - s[f'param_{param}_joint'])**2 + for param in ['baseline', 'alpha', 'mu', 'sigma']]).sum()) + + +def compute_norm2_error_mix(s): + return np.sqrt(np.array([(s[param] - s[f'param_{param}_mix'])**2 + for param in ['baseline', 'alpha', 'mu', 'sigma']]).sum()) + + +def compute_norm2_error_denoise(s): + return np.sqrt(np.array([(s[param] - s[f'param_{param}_denoise'])**2 + for param in ['baseline', 'alpha', 'mu', 'sigma']]).sum()) + + +df['err_norm2_jointfadin'] = df.apply( + lambda x: compute_norm2_error_joint(x), axis=1) +df['err_norm2_mixture'] = df.apply( + lambda x: compute_norm2_error_mix(x), axis=1) +df['err_norm2_denoise'] = df.apply( + lambda x: compute_norm2_error_denoise(x), axis=1) + + +df.to_csv(f'results/error_denoising_infer_{setting}.csv', index=False) diff --git a/experiments/unhap/simulated_data/run_rho_error.py b/experiments/unhap/simulated_data/run_rho_error.py new file mode 100644 index 0000000..d162209 --- /dev/null +++ b/experiments/unhap/simulated_data/run_rho_error.py @@ -0,0 +1,163 @@ +# %% import stuff +import numpy as np +import itertools +from joblib import Parallel, delayed +import pandas as pd +import torch +import time +from sklearn.metrics import precision_score, recall_score + +from fadin.kernels import DiscreteKernelFiniteSupport + +from unhap.utils.utils import smooth_projection_marked +from unhap.utils.utils_simu import simu_marked_hawkes_cluster, simu_multi_poisson +from unhap.solver import UNHaP + + +def simulate_data(baseline, baseline_noise, alpha, end_time, seed=0): + n_dim = len(baseline) + + def identity(x, **param): + return x + + def linear_zero_one(x, **params): + temp = 2 * x + mask = x > 1 + temp[mask] = 0. + return temp + + def truncated_gaussian(x, **params): + rc = DiscreteKernelFiniteSupport(delta=0.01, n_dim=1, + kernel='truncated_gaussian') + mu = params['mu'] + sigma = params['sigma'] + kernel_values = rc.kernel_eval( + [torch.Tensor(mu), torch.Tensor(sigma)], torch.tensor(x)) + + return kernel_values.double().numpy() + + marks_kernel = identity + marks_density = linear_zero_one + time_kernel = truncated_gaussian + + params_marks_density = dict() + # params_marks_density = dict(scale=1) + params_marks_kernel = dict(slope=1.2) + params_time_kernel = dict(mu=mu, sigma=sigma) + + marked_events, y_true = simu_marked_hawkes_cluster( + end_time, baseline, alpha, time_kernel, marks_kernel, marks_density, + params_marks_kernel=params_marks_kernel, + params_marks_density=params_marks_density, + time_kernel_length=None, marks_kernel_length=None, upper_bound=None, + params_time_kernel=params_time_kernel, random_state=seed) + + noisy_events_ = simu_multi_poisson(end_time, baseline_noise) + + random_marks = [ + np.random.rand(noisy_events_[i].shape[0]) for i in range(n_dim)] + noisy_events = [ + np.concatenate((noisy_events_[i].reshape(-1, 1), + random_marks[i].reshape(-1, 1)), axis=1) for i in range(n_dim)] + + events = [ + np.concatenate( + (noisy_events[i], marked_events[i]), axis=0) for i in range(n_dim)] + + events_cat = [events[i][events[i][:, 0].argsort()] for i in range(n_dim)] + n_grid = int(1 / 0.01) * end_time + 1 + events_grid, _, marked_events_smooth, index_unique = \ + smooth_projection_marked(events_cat, 0.01, n_grid) + + ########################################## + # Compute rho for the original events + labels = [] + false_events = np.where(y_true == 0)[0] + for i in range(n_dim): + a = marked_events[i].shape[0] + b = noisy_events_[i].shape[0] + labels_i = np.zeros(a+b) + labels_i[b:] = 1. + labels_i[false_events] = 0. + labels.append(labels_i) + true_rho = [labels[i][events[i][:, 0].argsort()] for i in range(n_dim)] + + # select the rho for the smoothed events + true_rho_smooth = [torch.tensor( + true_rho[i][index_unique]).float() for i in range(n_dim)] + loc_events = torch.where(events_grid > 0) + rho_star = events_grid.clone() + rho_star[loc_events] = torch.vstack(true_rho_smooth) + + return marked_events_smooth, rho_star, loc_events, labels, a + + +def run_experiment(baseline, baseline_noise, alpha, end_time, delta, seed): + + marked_events_smooth, rho_star, loc_events, _, _ = simulate_data( + baseline, baseline_noise, alpha, end_time=end_time, seed=seed) + start = time.time() + max_iter = 10000 + solver = UNHaP( + n_dim=1, + kernel="truncated_gaussian", + kernel_length=1., + delta=delta, optim="RMSprop", + params_optim={'lr': 1e-3}, + max_iter=max_iter, criterion='l2', + optimize_kernel=True, + optimize_alpha=True, optimize_rho=True, + batch_rho=200, + ) + + solver.fit(marked_events_smooth, end_time) + comp_time = time.time() - start + mean_rho_error = torch.abs( + solver.param_rho[loc_events] - rho_star[loc_events]).mean() + + y_pred = solver.param_rho[loc_events] + + y_true = rho_star[loc_events] + + rec_score = recall_score(y_true, y_pred) + pr_score = precision_score(y_true, y_pred) + + results = dict(mean_rho_error=mean_rho_error.item(), + rec_score=rec_score, + pr_score=pr_score) + + results["seed"] = seed + results['alpha'] = alpha.item() + results["end_time"] = end_time + results["delta"] = delta + results["noise"] = baseline_noise.item() + results["n_events"] = marked_events_smooth[0].shape[0] + results["comp_time"] = comp_time + + return results + + +baseline = np.array([.4]) +mu = np.array([[0.5]]) +sigma = np.array([[0.1]]) +delta = 0.01 + +end_time_list = [100, 1000, 10_000] +alpha_list = [0, 0.2, 0.3, 0.5, 0.7, 0.8, 0.9, 1.1, 1.2, 1.3, 1.47] +baseline_noise_list = [np.array([0.1]), np.array([0.5]), np.array([1.])] +seeds = np.arange(100) + +n_jobs = 70 +all_results = Parallel(n_jobs=n_jobs, verbose=10)( + delayed(run_experiment)(baseline, baseline_noise, + np.array([[alpha]]), end_time, + delta, seed=seed) + for end_time, alpha, baseline_noise, seed in itertools.product( + end_time_list, alpha_list, baseline_noise_list, seeds + ) +) + +# save results +df = pd.DataFrame(all_results) + +df.to_csv('results/error_rho_infer_5000.csv', index=False) From 464db70b8ffaacd3203a14e715c7ca0bfcaea88a Mon Sep 17 00:00:00 2001 From: vloison Date: Fri, 11 Jul 2025 19:13:08 +0200 Subject: [PATCH 29/41] linter fix --- .../unhap/simulated_data/run_infer_noise.py | 196 ++++++------------ 1 file changed, 65 insertions(+), 131 deletions(-) diff --git a/experiments/unhap/simulated_data/run_infer_noise.py b/experiments/unhap/simulated_data/run_infer_noise.py index 63db47d..9f5b721 100644 --- a/experiments/unhap/simulated_data/run_infer_noise.py +++ b/experiments/unhap/simulated_data/run_infer_noise.py @@ -623,137 +623,6 @@ def fit(self, events, end_time): return self -def get_grad_baseline_joint_noise2(zG, baseline, baseline_noise, alpha, kernel, - delta, end_time, sum_marks, - square_int_marks, n_ground_events): - """Return the gradient of the discrete l2 loss w.r.t. the baseline. - - .. math:: - N_T\\frac{\\partial\\mathcal{L}_G}{\\partial \\mu_{m}} = - 2 T \\mu_m - 2 N_T^m + 2 \\Delta\\sum_{j=1}^{p} \\sum_{\\tau=1}^{L} - \\phi_{mj}^\\Delta[\\tau]\\Phi_{j}(\\tau; G) - - Parameters - ---------- - zG : tensor, shape (n_dim, L) - - baseline : tensor, shape (n_dim,) - Baseline parameter of the intensity of the Hawkes process. - - alpha : tensor, shape (n_dim, n_dim) - Alpha parameter of the intensity of the Hawkes process. - - kernel : tensor, shape (n_dim, n_dim, L) - Kernel values on the discretization. - - delta : float - Step size of the discretization grid. - - sum_marks : tensor, shape (n_dim,) - Sum of marks for each dimension. - - end_time : float - The end time of the Hawkes process. - - n_ground_events : tensor, shape (n_dim,) - Number of events for each dimension. - - Returns - ---------- - grad_baseline: tensor, shape (dim,) - """ - n_dim, _, _ = kernel.shape - - if n_dim > 1: - cst1 = end_time * baseline * square_int_marks - cst2 = 0.5 * n_ground_events.sum() - cst3 = end_time * baseline_noise - - dot_kernel = torch.einsum('kju,ju->kj', kernel, zG) - dot_kernel_ = (dot_kernel * alpha).sum(1) * square_int_marks - - grad_baseline = (dot_kernel_ * delta + cst1 - sum_marks + cst3) / cst2 - else: - grad_baseline_ = torch.zeros(n_dim) - for k in range(n_dim): - temp = 0 - for j in range(n_dim): - temp += alpha[k, j] * (zG[j] @ kernel[k, j]) - grad_baseline_[k] = delta * temp * square_int_marks[k] - grad_baseline_[k] += end_time * baseline[k] * square_int_marks[k] - grad_baseline_[k] -= 2 * sum_marks[k] # add of joint MarkedFaDIn - grad_baseline_[k] += end_time * baseline_noise[k] - grad_baseline = 2 * grad_baseline_ / n_ground_events.sum() - - return grad_baseline - - -def get_grad_baseline_noise_joint_noise2(zG, baseline, baseline_noise, alpha, kernel, - delta, end_time, sum_marks, square_int_marks, - n_ground_events): - """Return the gradient of the discrete l2 loss w.r.t. the baseline. - - .. math:: - N_T\\frac{\\partial\\mathcal{L}_G}{\\partial \\mu_{m}} = - 2 T \\mu_m - 2 N_T^m + 2 \\Delta\\sum_{j=1}^{p} \\sum_{\\tau=1}^{L} - \\phi_{mj}^\\Delta[\\tau]\\Phi_{j}(\\tau; G) - - Parameters - ---------- - zG : tensor, shape (n_dim, L) - - baseline : tensor, shape (n_dim,) - Baseline parameter of the intensity of the Hawkes process. - - alpha : tensor, shape (n_dim, n_dim) - Alpha parameter of the intensity of the Hawkes process. - - kernel : tensor, shape (n_dim, n_dim, L) - Kernel values on the discretization. - - delta : float - Step size of the discretization grid. - - sum_marks : tensor, shape (n_dim,) - Sum of marks for each dimension. - - end_time : float - The end time of the Hawkes process. - - n_ground_events : tensor, shape (n_dim,) - Number of events for each dimension. - - Returns - ---------- - grad_baseline: tensor, shape (dim,) - """ - n_dim, _, _ = kernel.shape - - if n_dim > 1: - cst1 = end_time * baseline_noise - cst2 = 0.5 * n_ground_events.sum() - cst3 = end_time * baseline - - dot_kernel = torch.einsum('kju,ju->kj', kernel, zG) - dot_kernel_ = (dot_kernel * alpha).sum(1) - - grad_baseline_noise = ( - dot_kernel_ * delta + cst1 - n_ground_events + cst3) / cst2 - else: - grad_baseline_noise_ = torch.zeros(n_dim) - for k in range(n_dim): - temp = 0 - for j in range(n_dim): - temp += alpha[k, j] * (zG[j] @ kernel[k, j]) - grad_baseline_noise_[k] = delta * temp - grad_baseline_noise_[k] += end_time * baseline_noise[k] - grad_baseline_noise_[k] -= n_ground_events[k] - grad_baseline_noise_[k] += end_time * baseline[k] - grad_baseline_noise = 2 * grad_baseline_noise_ / n_ground_events.sum() - - return grad_baseline_noise - - def get_grad_baseline_noise_joint_noise2(zG, baseline, baseline_noise, alpha, kernel, delta, end_time, sum_marks, square_int_marks, n_ground_events): @@ -975,6 +844,71 @@ def get_grad_eta_joint_noise2(zG, zN, ztzG, baseline, baseline_noise, alpha, ker return grad_theta +def get_grad_baseline_joint_noise2(zG, baseline, baseline_noise, alpha, kernel, + delta, end_time, sum_marks, + square_int_marks, n_ground_events): + """Return the gradient of the discrete l2 loss w.r.t. the baseline. + + .. math:: + N_T\\frac{\\partial\\mathcal{L}_G}{\\partial \\mu_{m}} = + 2 T \\mu_m - 2 N_T^m + 2 \\Delta\\sum_{j=1}^{p} \\sum_{\\tau=1}^{L} + \\phi_{mj}^\\Delta[\\tau]\\Phi_{j}(\\tau; G) + + Parameters + ---------- + zG : tensor, shape (n_dim, L) + + baseline : tensor, shape (n_dim,) + Baseline parameter of the intensity of the Hawkes process. + + alpha : tensor, shape (n_dim, n_dim) + Alpha parameter of the intensity of the Hawkes process. + + kernel : tensor, shape (n_dim, n_dim, L) + Kernel values on the discretization. + + delta : float + Step size of the discretization grid. + + sum_marks : tensor, shape (n_dim,) + Sum of marks for each dimension. + + end_time : float + The end time of the Hawkes process. + + n_ground_events : tensor, shape (n_dim,) + Number of events for each dimension. + + Returns + ---------- + grad_baseline: tensor, shape (dim,) + """ + n_dim, _, _ = kernel.shape + + if n_dim > 1: + cst1 = end_time * baseline * square_int_marks + cst2 = 0.5 * n_ground_events.sum() + cst3 = end_time * baseline_noise + + dot_kernel = torch.einsum('kju,ju->kj', kernel, zG) + dot_kernel_ = (dot_kernel * alpha).sum(1) * square_int_marks + + grad_baseline = (dot_kernel_ * delta + cst1 - sum_marks + cst3) / cst2 + else: + grad_baseline_ = torch.zeros(n_dim) + for k in range(n_dim): + temp = 0 + for j in range(n_dim): + temp += alpha[k, j] * (zG[j] @ kernel[k, j]) + grad_baseline_[k] = delta * temp * square_int_marks[k] + grad_baseline_[k] += end_time * baseline[k] * square_int_marks[k] + grad_baseline_[k] -= 2 * sum_marks[k] # add of joint MarkedFaDIn + grad_baseline_[k] += end_time * baseline_noise[k] + grad_baseline = 2 * grad_baseline_ / n_ground_events.sum() + + return grad_baseline + + class JointFaDInDenoising(object): """Joint inference of the Hawkes process and the noise, without unmixing. From 34a60aecc3281da3e2e4009fbd21d4362531cf9f Mon Sep 17 00:00:00 2001 From: Virginie Loison <48970001+vloison@users.noreply.github.com> Date: Tue, 15 Jul 2025 12:44:26 +0200 Subject: [PATCH 30/41] Update README.md Co-authored-by: Thomas Moreau --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index a5adc58..0a425f1 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ ![license](https://img.shields.io/github/license/GuillaumeStaermanML/FaDIn?style=for-the-badge) ![code style](https://img.shields.io/badge/code_style-black-black?style=for-the-badge) -This package implements FaDIn and UNHaP. FaDIn and UNHaP are solvers for inference of Hawkes Processes with finite-support kernels on simulated or real data, with the following features: +This package implements FaDIn and UNHaP. FaDIn and UNHaP are inference methods for parametric Hawkes Processes (HP) with finite-support kernels, with the following features: - Computation time is low compared to other methods. - Compatible in an univariate setting as well as a multivariate setting. - Classical kernels (exponential truncated gaussian, raised cosine) are implemented. The user can add their own kernel for inference. From 8ad23c5235ebb6cb89e553dc1485c41ff6183d9a Mon Sep 17 00:00:00 2001 From: Virginie Loison <48970001+vloison@users.noreply.github.com> Date: Tue, 15 Jul 2025 12:44:51 +0200 Subject: [PATCH 31/41] Update README.md Co-authored-by: Thomas Moreau --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 0a425f1..54d93ce 100644 --- a/README.md +++ b/README.md @@ -7,7 +7,7 @@ This package implements FaDIn and UNHaP. FaDIn and UNHaP are inference methods for parametric Hawkes Processes (HP) with finite-support kernels, with the following features: - Computation time is low compared to other methods. -- Compatible in an univariate setting as well as a multivariate setting. +- Compatible in univariate and multivariate settings. - Classical kernels (exponential truncated gaussian, raised cosine) are implemented. The user can add their own kernel for inference. - Masking: if only a few Hawkes Parameters need to be inferred, the user can mask the other parameters. - Smart initialization of parameters before optimization: the user can choose between `random` (purely random), `moment_matching_max` (moment matching with maximum mode) and `moment_matching_mean` (moment matching with mean mode). From 9f7e9cc622663d7b987e8a71a1336becbcab0552 Mon Sep 17 00:00:00 2001 From: Virginie Loison <48970001+vloison@users.noreply.github.com> Date: Tue, 15 Jul 2025 12:45:06 +0200 Subject: [PATCH 32/41] Update README.md Co-authored-by: Thomas Moreau --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 54d93ce..f4b8115 100644 --- a/README.md +++ b/README.md @@ -8,7 +8,7 @@ This package implements FaDIn and UNHaP. FaDIn and UNHaP are inference methods for parametric Hawkes Processes (HP) with finite-support kernels, with the following features: - Computation time is low compared to other methods. - Compatible in univariate and multivariate settings. -- Classical kernels (exponential truncated gaussian, raised cosine) are implemented. The user can add their own kernel for inference. +- *Flexible:* Various kernel choices are implemented, with classical ones (exponential, truncated Gaussian, raised cosine) and an API to add custom kernels for inference. - Masking: if only a few Hawkes Parameters need to be inferred, the user can mask the other parameters. - Smart initialization of parameters before optimization: the user can choose between `random` (purely random), `moment_matching_max` (moment matching with maximum mode) and `moment_matching_mean` (moment matching with mean mode). From 55b95b98e72d4f0eeec8674727413fde0bfe1c5e Mon Sep 17 00:00:00 2001 From: Virginie Loison <48970001+vloison@users.noreply.github.com> Date: Tue, 15 Jul 2025 12:45:21 +0200 Subject: [PATCH 33/41] Update README.md Co-authored-by: Thomas Moreau --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index f4b8115..e07893b 100644 --- a/README.md +++ b/README.md @@ -72,7 +72,7 @@ If this package was useful to you, please cite it in your work: organization={PMLR} } -@improceedings{loison2025unhap, +@inproceedings{loison2025unhap, title={UNHaP: Unmixing Noise from Hawkes Process}, author={Loison, Virginie and Staerman, Guillaume and Moreau, Thomas}, booktitle={International Conference on Artificial Intelligence and Statistics}, From 522090379e1e059d81c3f5b635daffbb274744e1 Mon Sep 17 00:00:00 2001 From: Virginie Loison <48970001+vloison@users.noreply.github.com> Date: Tue, 15 Jul 2025 12:45:47 +0200 Subject: [PATCH 34/41] Update fadin/utils/utils_simu.py Co-authored-by: Thomas Moreau --- fadin/utils/utils_simu.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/fadin/utils/utils_simu.py b/fadin/utils/utils_simu.py index cd2dad0..04e219a 100644 --- a/fadin/utils/utils_simu.py +++ b/fadin/utils/utils_simu.py @@ -248,8 +248,7 @@ def simu_hawkes_cluster(end_time, baseline, alpha, kernel, References: Møller, J., & Rasmussen, J. G. (2006). Approximate simulation of Hawkes - processes. - Methodology and Computing in Applied Probability, 8, 53-64. + processes. Methodology and Computing in Applied Probability, 8, 53-64. Møller, J., & Rasmussen, J. G. (2005). Perfect simulation of Hawkes processes. From bbd4eddc78c8290ccf2d975c773134143aea040b Mon Sep 17 00:00:00 2001 From: Virginie Loison <48970001+vloison@users.noreply.github.com> Date: Tue, 15 Jul 2025 12:46:09 +0200 Subject: [PATCH 35/41] Update fadin/utils/utils_simu.py Co-authored-by: Thomas Moreau --- fadin/utils/utils_simu.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/fadin/utils/utils_simu.py b/fadin/utils/utils_simu.py index 04e219a..d670f71 100644 --- a/fadin/utils/utils_simu.py +++ b/fadin/utils/utils_simu.py @@ -182,10 +182,9 @@ def simu_multi_poisson(end_time, intensity, upper_bound=None, random_state=None): """Simulate multivariate Poisson processes on [0, end_time] with the Ogata's modified thinning algorithm by superposition of univariate - processes. - If the intensity is a numerical value, simulate a Homegenous Poisson - Process. - If the intensity is a function, simulate an Inhomogenous Poisson Process. + processes. If the intensity is a numerical value, simulate a Homegenous + Poisson Process. If the intensity is a function, simulate an + Inhomogenous Poisson Process. Parameters ---------- From a1fc636f0e0a935d806bc87a2c200151d3005663 Mon Sep 17 00:00:00 2001 From: Virginie Loison <48970001+vloison@users.noreply.github.com> Date: Tue, 15 Jul 2025 12:46:27 +0200 Subject: [PATCH 36/41] Update README.md Co-authored-by: Thomas Moreau --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index e07893b..644084f 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ This package implements FaDIn and UNHaP. FaDIn and UNHaP are inference methods f - Computation time is low compared to other methods. - Compatible in univariate and multivariate settings. - *Flexible:* Various kernel choices are implemented, with classical ones (exponential, truncated Gaussian, raised cosine) and an API to add custom kernels for inference. -- Masking: if only a few Hawkes Parameters need to be inferred, the user can mask the other parameters. +- *Masking:* if some parameters can be fixed, the user can mask them easily. - Smart initialization of parameters before optimization: the user can choose between `random` (purely random), `moment_matching_max` (moment matching with maximum mode) and `moment_matching_mean` (moment matching with mean mode). From c762783fb490700a2746241e00601fcad300924f Mon Sep 17 00:00:00 2001 From: Virginie Loison <48970001+vloison@users.noreply.github.com> Date: Tue, 15 Jul 2025 12:46:57 +0200 Subject: [PATCH 37/41] Update fadin/utils/utils_simu.py Co-authored-by: Thomas Moreau --- fadin/utils/utils_simu.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fadin/utils/utils_simu.py b/fadin/utils/utils_simu.py index d670f71..8e6954b 100644 --- a/fadin/utils/utils_simu.py +++ b/fadin/utils/utils_simu.py @@ -241,8 +241,8 @@ def simu_hawkes_cluster(end_time, baseline, alpha, kernel, params_kernel=dict(), kernel_length=None, upper_bound=None, random_state=None): """ Simulate a multivariate Hawkes process following an immigration-birth - procedure. - Edge effects may be reduced according to the second references below. + procedure. Edge effects may be reduced according to the second + references below. References: From dc19cf58525dfb367153551ecbc745f54171b491 Mon Sep 17 00:00:00 2001 From: Virginie Loison <48970001+vloison@users.noreply.github.com> Date: Tue, 15 Jul 2025 12:47:23 +0200 Subject: [PATCH 38/41] Update fadin/utils/utils_simu.py Co-authored-by: Thomas Moreau --- fadin/utils/utils_simu.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/fadin/utils/utils_simu.py b/fadin/utils/utils_simu.py index 8e6954b..0e1c3ae 100644 --- a/fadin/utils/utils_simu.py +++ b/fadin/utils/utils_simu.py @@ -250,8 +250,7 @@ def simu_hawkes_cluster(end_time, baseline, alpha, kernel, processes. Methodology and Computing in Applied Probability, 8, 53-64. Møller, J., & Rasmussen, J. G. (2005). Perfect simulation of Hawkes - processes. - Advances in applied probability, 37(3), 629-646. + processes. Advances in applied probability, 37(3), 629-646. Parameters ---------- From 50c2538d8a20a4b673a5e1b5f505449113d13bf4 Mon Sep 17 00:00:00 2001 From: vloison Date: Tue, 15 Jul 2025 15:44:58 +0200 Subject: [PATCH 39/41] minor changes to readme --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 644084f..43a7360 100644 --- a/README.md +++ b/README.md @@ -6,11 +6,11 @@ ![code style](https://img.shields.io/badge/code_style-black-black?style=for-the-badge) This package implements FaDIn and UNHaP. FaDIn and UNHaP are inference methods for parametric Hawkes Processes (HP) with finite-support kernels, with the following features: -- Computation time is low compared to other methods. +- *Fast:* computation time is low compared to other methods. - Compatible in univariate and multivariate settings. -- *Flexible:* Various kernel choices are implemented, with classical ones (exponential, truncated Gaussian, raised cosine) and an API to add custom kernels for inference. +- *Flexible:* various kernel choices are implemented, with classical ones (exponential, truncated Gaussian, raised cosine) and an API to add custom kernels for inference. - *Masking:* if some parameters can be fixed, the user can mask them easily. -- Smart initialization of parameters before optimization: the user can choose between `random` (purely random), `moment_matching_max` (moment matching with maximum mode) and `moment_matching_mean` (moment matching with mean mode). +- Smart initialization of parameters before optimization: the user can choose between `random` (purely random), `moment_matching_max` (moment matching with maximum mode) and `moment_matching_mean` (moment matching with mean mode). The moment matching options are implemented for UNHaP. [FaDIn](https://proceedings.mlr.press/v202/staerman23a/staerman23a.pdf) does classical Hawkes inference with gradient descent. From 409eb6561c4bb7dced897cbe11cbd55b89e5301e Mon Sep 17 00:00:00 2001 From: vloison Date: Tue, 15 Jul 2025 15:45:18 +0200 Subject: [PATCH 40/41] add doc for events --- examples/plot_unhap.py | 15 +++++++++++++++ examples/plot_univariate_fadin.py | 11 +++++++++++ fadin/solver.py | 20 +++++++++++++++----- 3 files changed, 41 insertions(+), 5 deletions(-) diff --git a/examples/plot_unhap.py b/examples/plot_unhap.py index 01eed23..bf3cdba 100644 --- a/examples/plot_unhap.py +++ b/examples/plot_unhap.py @@ -44,6 +44,21 @@ ev, noisy_marks, true_rho = simulate_marked_data( baseline, baseline_noise.item(), alpha, end_time, mu, sigma, seed=0 ) + +# %% Let's take a closer look at the events +print('Type of ev object', type(ev)) +# ev is a list of numpy arrays, one for each dimension +print('Number of events', len(ev[0])) +print('Shape of first event array', ev[0].shape) +# Each dimension is stored as a numpy array of shape (n_events, 2). +print('First 10 events timestamps and marks', ev[0][:10]) +# Each event is stored as [timestamp, mark]. +# This is the expected data format for UNHaP. +print('First event timestamp', ev[0][0][0]) +print('First event mark', ev[0][0][1]) +print('Second event timestamp', ev[0][1][0]) +print('Second event mark', ev[0][1][1]) + # %% Initiate and fit UNHAP to the simulated events solver = UNHaP( diff --git a/examples/plot_univariate_fadin.py b/examples/plot_univariate_fadin.py index 45b1d17..5536132 100644 --- a/examples/plot_univariate_fadin.py +++ b/examples/plot_univariate_fadin.py @@ -48,6 +48,17 @@ events = simu_hawkes_cluster(T, baseline, alpha, kernel, params_kernel={'scale': 1 / beta}) +############################################################################### +# Let's take a closer look at the events +print('Type of events object', type(events)) +# events is a list of numpy arrays, one for each dimension +print('Number of events', len(events[0])) +print('First 10 events timestamps', events[0][:10]) +# For each event, its occurence time (timestamp) is stored in the numpy array. + +# `events`` is a list. Its elements are numpy arrays containing the timestamps +# of the events. This is the expected data format for FaDIn. + ############################################################################### # Here, we apply FaDIn. diff --git a/fadin/solver.py b/fadin/solver.py index 84b879c..e6c5ae8 100644 --- a/fadin/solver.py +++ b/fadin/solver.py @@ -221,8 +221,10 @@ def fit(self, events, end_time): Parameters ---------- - events : list of array of size number of timestamps, - list size is self.n_dim. + events : list of `self.n_dim` numpy arrays. + One numpy array = one dimension. One numpy array + has shape (n_events, 1), where n_events is the number of events in this + dimension. The timestamp of each event is stored. end_time : int The end time of the Hawkes process. @@ -359,7 +361,13 @@ class UNHaP(object): """Define the UNHaP framework for estimated mixture of Hawkes and Poisson processes. - Unhap minimizes the discretized L2 mixture loss of a mixture of Hawkes and + The framework is detailed in: + + Virginie Loison, Guillaume Staerman, Thomas Moreau + UNHaP: Unmixing Noise from Hawkes Processes + https://proceedings.mlr.press/v258/loison25a.html + + UNHaP minimizes the discretized L2 mixture loss of a mixture of Hawkes and Poisson processes. Parameters @@ -578,8 +586,10 @@ def fit(self, events, end_time): Parameters ---------- - events : pandas DataFrame with three columns: timestamps, marks values - and type of events. + events : list of `self.n_dim` numpy arrays. + One numpy array = one dimension. One numpy array + has shape (n_events, 2), where n_events is the number of events in this + dimension. Each event is stored as (timestamp, mark). end_time : int The end time of the Hawkes process. From ddcb368c2f11381c8688c01ac7e1006cae8352b8 Mon Sep 17 00:00:00 2001 From: vloison Date: Tue, 15 Jul 2025 17:41:17 +0200 Subject: [PATCH 41/41] update pyproject.toml and readme with experiment dependencies --- README.md | 6 ++++++ pyproject.toml | 8 ++++++++ 2 files changed, 14 insertions(+) diff --git a/README.md b/README.md index 43a7360..f7b800d 100644 --- a/README.md +++ b/README.md @@ -52,6 +52,12 @@ pip install -e ".[dev]" pre-commit install ``` +Before running the experiments of the FaDIn and UNHaP papers located in the `experiments` directory, please install the corresponding dependencies beforehand. + +```bash +pip install -e ".[experiments]" +``` + ## Short examples A few illustrative examples are provided in the `examples` folder of this repository, in particular: - `plot_univariate_fadin`: simulate an univariate unmarked Hawkes process, infer Hawkes Process parameters using FaDIn, and plot inferred kernel. diff --git a/pyproject.toml b/pyproject.toml index c2124cc..680d80a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,7 +43,15 @@ test = [ experiments = [ "dicodile>=0.3", + "vitaldb>=1.4.5", + "alphacsc>=0.4.0", + "pyhrv>=0.4.1", + "biosppy>=2.2.3", + "neurokit2>=0.2.7", + "plotly>=6.2.0", + "peakutils>=1.3.5", ] + doc = [ "furo", "matplotlib",