ENH add UNHaP experiments and documentation (#21)

Co-authored-by: Thomas Moreau <thomas.moreau.2010@gmail.com>

Author: vloison

README.md

@@ -1,17 +1,28 @@
-# FaDIn: Fast Discretized Inference For Hawkes Processes with General Parametric Kernels
+# 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)
 ![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 inference methods for parametric Hawkes Processes (HP) with finite-support kernels, with the following features:
+- *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. 
+- *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). The moment matching options are implemented for UNHaP.
 
-## Installation
 
-**To install this package, make sure you have an up-to-date version of** `pip`.
+[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`.
+```bash
+python3 -m pip install --upgrade pip
+```
 ### From PyPI (coming soon)
 
 In a dedicated Python env, run:
@@ -41,6 +52,18 @@ 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.
+- `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 +77,12 @@ If this package was useful to you, please cite it in your work:
   year={2023},
   organization={PMLR}
 }
+
+@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},
+pages={1342--1350},
+year={2025}
+}
 ```

examples/plot_multivariate_fadin.py

@@ -49,8 +49,8 @@
 events = simu_hawkes_cluster(T, baseline, alpha, kernel)
 
 ###############################################################################
-# Here, we apply FaDIn.
-
+# 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,25 +60,20 @@
     max_iter=10000
 )
 solver.fit(events, T)
+print("FaDIn solver fitted.")
+# We can now access the estimated parameters of the model.
 
-# We average on the 10 last values of the optimization.
-
-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):

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,120 +34,74 @@
 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_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
-
-
-ev, noisy_marks, true_rho = simulate_data(
-    baseline, baseline_noise.item(), alpha, end_time, seed=0
+# %% 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, mu, sigma, seed=0
 )
-# %% Apply UNHAP
+
+# %% 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(
     n_dim=1,
     kernel="truncated_gaussian",
     kernel_length=1.0,
+    init='moment_matching_mean',
     delta=delta,
     optim="RMSprop",
     params_optim={"lr": 1e-3},
     max_iter=max_iter,
     batch_rho=batch_rho,
     density_hawkes="linear",
     density_noise="uniform",
-    moment_matching=True,
 )
 solver.fit(ev, end_time)
 
 # %% 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 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)
-# %%

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.
 
@@ -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.
 
@@ -60,24 +71,20 @@
                )
 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_)
 
 
 ###############################################################################
 # Here, we plot the values of the estimated kernel with FaDIn.
 
 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')