diff --git a/README.md b/README.md index 5490ea4..f7b800d 100644 --- a/README.md +++ b/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} +} ``` diff --git a/examples/plot_multivariate_fadin.py b/examples/plot_multivariate_fadin.py index 8e1f886..e01b745 100644 --- a/examples/plot_multivariate_fadin.py +++ b/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): diff --git a/examples/plot_unhap.py b/examples/plot_unhap.py index d8bdf60..bf3cdba 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,88 +34,38 @@ 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}, @@ -126,31 +73,35 @@ 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) # %% 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) -# %% diff --git a/examples/plot_univariate_fadin.py b/examples/plot_univariate_fadin.py index 9c47e61..5536132 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. @@ -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,15 +71,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 +84,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/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/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_ecg.py b/experiments/unhap/ecg_gait/utils_ecg.py new file mode 100644 index 0000000..29b6c4b --- /dev/null +++ b/experiments/unhap/ecg_gait/utils_ecg.py @@ -0,0 +1,782 @@ +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.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.kernel_ + + if param_kernel[0] == 0: + abs_error = np.nan + else: + 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_ + else: + intens_bl = estimated_baseline + _, events_grid = projected_grid_marked(events, delta, L * T + 1) + intens = evaluate_intensity( + [intens_bl], + solver.alpha_, + solver.kernel_[0], + solver.kernel_[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].item(), decimals=4) + ) + print( + "Estimated std of gaussian kernel:", + np.round(param_kernel[1].item(), 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 diff --git a/experiments/unhap/ecg_gait/utils_gait.py b/experiments/unhap/ecg_gait/utils_gait.py new file mode 100644 index 0000000..02e1071 --- /dev/null +++ b/experiments/unhap/ecg_gait/utils_gait.py @@ -0,0 +1,104 @@ + + +import numpy as np + + +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 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/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..9f5b721 --- /dev/null +++ b/experiments/unhap/simulated_data/run_infer_noise.py @@ -0,0 +1,1504 @@ +# %% 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_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 + + +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. + + 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) diff --git a/fadin/init.py b/fadin/init.py index 271e1d6..25588e1 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,9 @@ 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: 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 @@ -268,6 +271,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 +316,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/loss_and_gradient.py b/fadin/loss_and_gradient.py index 9620108..1364b6f 100644 --- a/fadin/loss_and_gradient.py +++ b/fadin/loss_and_gradient.py @@ -27,52 +27,54 @@ 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 ) - 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._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, @@ -148,6 +150,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 @@ -548,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/solver.py b/fadin/solver.py index 8767cf5..e6c5ae8 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,33 @@ 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. - - param_alpha : `tensor`, shape (n_dim, n_dim) - Weight parameter of the Hawkes process. - param_kernel : `list` of `tensor` - list containing tensor array of kernels parameters. + 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) - If `log=True`, compute the loss accross iterations. + v_loss_ : `tensor`, shape (n_iter) + loss accross iterations. If no early stopping, `n_iter` is equal to `max_iter`. """ compute_gradient = staticmethod(compute_gradient_fadin) @@ -136,8 +151,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 @@ -146,10 +160,10 @@ 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 - self.log = log self.tol = tol self.n_dim = n_dim self.kernel_model = DiscreteKernelFiniteSupport( @@ -207,15 +221,16 @@ 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. Returns ------- - TODO: attributes self : object Fitted parameters. """ @@ -228,7 +243,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, @@ -238,7 +254,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 ) @@ -257,18 +273,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() @@ -284,45 +305,70 @@ 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) - + 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 Poisson processes. - Unhap minimizes the discretized L2 mixture loss of Hawkes and Poisson processes. + 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 ---------- @@ -347,18 +393,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,40 +422,42 @@ 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. + 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'. @@ -410,26 +465,41 @@ class UNHaP(object): 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. - param_kernel : `list` of `tensor` - list containing tensor array of kernels parameters. + 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. + 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. + 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`. + """ 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, + density_noise='uniform', stoc_classif=False, random_state=None): # Set discretization parameters @@ -439,24 +509,30 @@ 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 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 - self.moment_matching = moment_matching - - # Set model parameters - self.n_dim = n_dim - - self.baseline = torch.rand(self.n_dim) + self.stoc_classif = stoc_classif # 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]) @@ -478,7 +554,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): @@ -493,22 +568,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(): @@ -527,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. @@ -538,7 +599,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, _ = \ @@ -572,19 +633,24 @@ 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 + # _params_intens: [baseline, baseline_noise, alpha, kernel_params] + 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.rho_ = self.params_mixture[0].detach() #################################################### # Precomputations @@ -595,7 +661,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) @@ -603,22 +669,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() @@ -630,7 +698,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 ) @@ -639,25 +707,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.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, @@ -675,7 +743,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, @@ -683,32 +751,70 @@ 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 = torch.round(self.params_mixture[0].detach()) + self.rho_, marked_quantities[1], self.L) + if not self.stoc_classif: + # vanilla UNHaP + self.rho_ = torch.round( + self.params_mixture[0].detach() + ) + else: + # StocUNHaP + random_rho = torch.rand(self.n_dim, n_grid) + self.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) - + # 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)] 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" diff --git a/fadin/tests/test_init.py b/fadin/tests/test_init.py new file mode 100644 index 0000000..06e0106 --- /dev/null +++ b/fadin/tests/test_init.py @@ -0,0 +1,52 @@ +import numpy as np + +from fadin.utils.utils_simu import simulate_marked_data +from fadin.solver import UNHaP + + +def test_unhap_init(): + + 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 + ) + unhap_mmmax = UNHaP( + n_dim=1, + kernel="truncated_gaussian", + init="moment_matching_max", + max_iter=max_iter, + batch_rho=batch_rho + ) + unhap_mmmax.fit(ev, end_time) + + unhap_mmtimean = UNHaP( + n_dim=1, + kernel="truncated_gaussian", + init="moment_matching_mean", + max_iter=max_iter, + batch_rho=batch_rho + ) + unhap_mmtimean.fit(ev, end_time) + + unhap_random = UNHaP( + n_dim=1, + kernel="truncated_gaussian", + init="random", + 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" + assert unhap_mmtimean is not None, "UNHaP moment matching mean failed" + assert unhap_random is not None, "UNHaP random initialization failed" + return None diff --git a/fadin/tests/test_mask.py b/fadin/tests/test_mask.py index 5dfb23a..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] @@ -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)) 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 " ) 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.""" diff --git a/fadin/utils/utils_simu.py b/fadin/utils/utils_simu.py index 6aae9e9..0e1c3ae 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,11 +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, - If the intensity is a function, simulate an Inhomogenous 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,16 +240,17 @@ 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. - Edge effects may be reduced according to the second references below. + """ 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. - Methodology and Computing in Applied Probability, 8, 53-64. + 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. - Advances in applied probability, 37(3), 629-646. + Møller, J., & Rasmussen, J. G. (2005). Perfect simulation of Hawkes + processes. Advances in applied probability, 37(3), 629-646. Parameters ---------- @@ -251,7 +265,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 +314,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 +426,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 +531,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 +541,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 +585,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 diff --git a/fadin/utils/vis.py b/fadin/utils/vis.py index 07627bf..0faee76 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_[0].item(), 2) axs[i, j].vlines( x=mean, ymin=0, diff --git a/pyproject.toml b/pyproject.toml index e6f5458..680d80a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,6 +41,17 @@ test = [ "pytest>=7.2", ] +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",