<jupyter_start><jupyter_text>Experiments on Synthetic Data Import libs<jupyter_code># External libs
from pprint import pprint
import functools
import itertools
import numpy as np
import scipy.stats
import pandas as pd
from matplotlib import pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.colors as colors
%matplotlib inline

import sys
if '../' not in sys.path:
    sys.path.append('../')

from tick.hawkes.simulation import SimuHawkesExpKernels
# lib for: ADM4, NPHC, WH
from tick.hawkes.inference import HawkesConditionalLaw, HawkesADM4, HawkesCumulantMatching
from tick.dataset import fetch_hawkes_bund_data
# lib for: Desync-MLE 
from desync_mhp.lib.inference import HawkesExpKernConditionalMLE

# Internal lib
import lib

%matplotlib inline<jupyter_output><empty_output><jupyter_text>--- 1. Define model and generate data<jupyter_code># Define the ground-truth parameters
adjacency = np.array([ [0.23, 0.23, 0.23, 0.23, 0.23, 0.  , 0.  , 0.  , 0.  , 0.23],
                       [0.  , 0.23, 0.23, 0.23, 0.23, 0.  , 0.  , 0.  , 0.23, 0.  ],
                       [0.  , 0.  , 0.23, 0.23, 0.23, 0.  , 0.  , 0.  , 0.  , 0.  ],
                       [0.  , 0.  , 0.  , 0.23, 0.23, 0.  , 0.  , 0.  , 0.  , 0.  ],
                       [0.  , 0.  , 0.  , 0.  , 0.23, 0.  , 0.  , 0.  , 0.  , 0.  ],
                       [0.  , 0.  , 0.  , 0.  , 0.  , 0.23, 0.  , 0.  , 0.  , 0.  ],
                       [0.  , 0.23, 0.  , 0.  , 0.  , 0.23, 0.23, 0.  , 0.  , 0.  ],
                       [0.23, 0.  , 0.  , 0.  , 0.  , 0.23, 0.23, 0.23, 0.  , 0.  ],
                       [0.  , 0.  , 0.  , 0.  , 0.  , 0.23, 0.23, 0.23, 0.23, 0.  ],
                       [0.  , 0.  , 0.  , 0.  , 0.  , 0.23, 0.23, 0.23, 0.23, 0.23] ])
decays = 1.0
baseline = np.array([0.01] * len(adjacency))
# Compute the ground-truth cumulants
L_true, C_true, Kc_true = lib.utils.cumulants.compute_cumulants(G=adjacency, mus=baseline,)


def simulate(noise_scale, seed):
    """Simulate a randomly translated realization of the process"""
    # Define the (noiseless) MHP simulation object
    simu_hawkes = SimuHawkesExpKernels(adjacency=adjacency, decays=decays, baseline=baseline, max_jumps=0, verbose=False)
    # Define noise distributions
    noise_rand_state = np.random.RandomState(seed=None)
    noise_dist_arr =  ['gaussian' for _ in range(simu_hawkes.n_nodes)]
    noise_scale_arr = [noise_scale for _ in range(simu_hawkes.n_nodes)]
    # Build noisy simulation object
    simu_noisy_hawkes = lib.simulation.noisy_hawkes.SimulatorNoisyHawkesCustomKernels(
        simu_obj=simu_hawkes,
        noise_dist=noise_dist_arr,
        noise_scale=noise_scale_arr,
        burn_in_quantile=0.99,
        num_real=1,
        num_jumps=100000,
        seed=seed,
        no_multi=False)
    # Simulate data
    noisy_events, orig_events = simu_noisy_hawkes.simulate(return_original_events = True)
    return noisy_events


def get_best_integration_support(events, max_iter=20, initial_simplex=[[10.0], [50.0]], verbose=False):
    """Hyper-parameter tuning for the bandwidht of the cumulant estimator in NPHC"""
    def int_support_loss(H, events):
        nphc = HawkesCumulantMatching(integration_support=float(H), max_iter=0, verbose=False)
        nphc.fit(events)
        skew_loss = np.linalg.norm(nphc.skewness - Kc_true, ord=2)
        cov_loss = np.linalg.norm(nphc.covariance - C_true, ord=2)
        norm_sq_K_c = np.linalg.norm(nphc.skewness, ord=2)**2 
        norm_sq_C = np.linalg.norm(nphc.covariance, ord=2)**2 
        cs_ratio = norm_sq_K_c / (norm_sq_K_c + norm_sq_C)
        loss = (1-cs_ratio) * skew_loss + cs_ratio * cov_loss
        if verbose:
            print(f"{float(H):>6.2f}, loss={loss:.2e}, skew_loss={skew_loss:.2e}, cov_loss={cov_loss:.2e}")
        return skew_loss
    res = scipy.optimize.minimize(
        int_support_loss, 
        x0=20.0, 
        args=(events,), 
        options={'max_iter': max_iter, 
                 'maxfev': max_iter, 
                 'initial_simplex': initial_simplex}, 
        method='Nelder-Mead')
    return float(res.x)<jupyter_output><empty_output><jupyter_text>--- 2. Run the Proof-Of-Concept for a fixed noise scale $\sigma^2 = 5$<jupyter_code>print('Simulate data...')
noisy_events = simulate(noise_scale=5.0, seed=534543)
print('done.')

# ADM4
adm4 = HawkesADM4(decay=1.0, verbose=True)
adm4.fit(noisy_events)
print("ADM4: L2-dist: done.")

# WH
wh = HawkesConditionalLaw(delta_lag=0.1, min_lag=0.0001, max_lag=100.0,
                          n_quad=20, max_support=10.0)
wh.fit(noisy_events)
wh_adj = wh.get_kernel_norms()
print("WH: L2-dist: done.")

# NPHC
nphc = HawkesCumulantMatching(integration_support=15.0, solver='adam', max_iter=10000, 
                              penalty='none', verbose=True)
nphc.fit(noisy_events)
print("NPHC: done.")

# Desync-MLE
end_time = max([max(map(max, real)) for real in noisy_events])
dim = len(baseline)
desyncmle = HawkesExpKernConditionalMLE(
    decay=1.0,
    noise_penalty='l2', noise_C=1e3,
    hawkes_penalty='l1', hawkes_base_C=1e2, hawkes_adj_C=1e5,
    solver='sgd', tol=1e-4, max_iter=1000,
    verbose=False
)
desyncmle.fit(noisy_events[0], end_time=end_time,
            z_start=np.zeros(dim),
            theta_start=np.hstack((
                0.01*np.ones(dim),
                np.random.uniform(0.0, 0.1, size=dim**2)
            )),
            callback=None)
desyncmle_adj = np.reshape(desyncmle.coeffs[2*dim:], (dim, dim))
print("Desync-MLE: done.")
fig = lib.utils.plotting.plotmat_sidebyside(
    {'(a) Ground-truth': adjacency, 
     '(b) ADM4': adm4.adjacency,
     '(c) Desync-MLE': desyncmle_adj,
     '(d) WH': wh_adj,
     '(e) NPHC': nphc.adjacency, 
    }, grid=(1, 5), figsize=(10, 3.25), ytitle=1.0, ticks=[1, 5, 10])
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>--- 2. Estimate the excitation matrix for varying noise levels  2.1. Run the simulationsBecause the experiments are time consuming, there are run in a separate script.<jupyter_code>%run script_run_synthetic_experiments.py<jupyter_output><empty_output><jupyter_text>--- 3. Visualize results Load the results from the experiment script<jupyter_code>df = pd.read_pickle('res-synthetic.pkl')
df.head()<jupyter_output><empty_output><jupyter_text>Compute the evaluation metrics<jupyter_code>method_queries = [
    ('adm4', 'ADM4'),
    ('nphc', 'NPHC'),
    ('wh',   'WH'),
    ('desyncmle', 'Desync-MLE'),
]

metric_queries = [
    ('norm',     'Norm',           lambda y_test, y_true: np.linalg.norm(y_test.ravel(), ord=2)**2                             ),
    ('relerr',   'Relative Error', lambda y_test, y_true: lib.utils.metrics.relerr(y_test, y_true, null_norm='min')       ),
    ('precAt5',  'Precison@5',     lambda y_test, y_true: lib.utils.metrics.precision_at_n(y_test.ravel(), y_true.ravel(), n=5)           ),
    ('precAt10', 'Precison@10',    lambda y_test, y_true: lib.utils.metrics.precision_at_n(y_test.ravel(), y_true.ravel(), n=10)          ),
    ('precAt20', 'Precison@20',    lambda y_test, y_true: lib.utils.metrics.precision_at_n(y_test.ravel(), y_true.ravel(), n=20)          ),
    ('pr-auc',   'PR-AUC',         lambda y_test, y_true: lib.utils.metrics.pr_auc_score(y_test.ravel(), y_true.ravel())                  ),
]

for method, _ in method_queries:
    for suffix, _, func in metric_queries:
        df[f'{method}_{suffix}'] = df[f'{method}_adj'].apply(lambda adj: func(y_test=adj, y_true=adjacency))
        
df['cov_l2'] = df['cov'].apply(lambda val: np.linalg.norm(C_true.T - val, ord=2))
df['skew_l2'] = df['skew'].apply(lambda val: np.linalg.norm(Kc_true.T - val, ord=2))

df_plot = df.groupby('noise_scale').agg(['mean', 'std', 'count'])
df_plot<jupyter_output><empty_output><jupyter_text>3.1. Plot the Cumulant Stability<jupyter_code>plt.figure(figsize=(6.75, 6.75 * 2/3))
plt.errorbar(df_plot.index, 
             df_plot.cov_l2['mean'],
             yerr=df_plot.cov_l2['std']/np.sqrt(df_plot.cov_l2['count']), 
             label=r'$2^{\mathrm{nd}}$-order cumulants (Covariance)')
plt.errorbar(df_plot.index, 
             df_plot.skew_l2['mean'],
             yerr=df_plot.skew_l2['std']/np.sqrt(df_plot.skew_l2['count']), 
             label=r'$3^{\mathrm{rd}}$-order cumulants(Skewness)')
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'Noise scale $\sigma^2$')
plt.ylabel(('Average L2-distance to\n'
            'ground-truth cumulants\n'));
plt.legend(loc='upper left')
plt.ylim(top=0.3);<jupyter_output><empty_output><jupyter_text>3.2. Plot the Estimation Performance w.r.t. Varying Noise Levels<jupyter_code>for suffix, metric_label, _ in metric_queries:    
    print(metric_label, flush=True)
    plt.figure(figsize=(6.75/4, 3.25*0.5*0.85))
    plt.grid()
    for method, method_label in method_queries:
        yseries = df_plot[f'{method}_{suffix}']
        y = yseries['mean']
        #yerr = yseries['std']
        yerr = yseries['std'] / np.sqrt(yseries['count'])
        plt.errorbar(df_plot.index, y, yerr=yerr, label=method_label, )
    if suffix.startswith('relerr'):
        # plt.legend();
        pass
    if suffix.startswith('prec'):
        plt.legend();
        pass
    plt.xscale('log')
    plt.xlabel(f'Noise scale $\sigma^2$')
    plt.ylabel(metric_label)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15285530270502362, 'fraction_numerical': 0.03885787891799056, 'mean_word_length': 3.0211480362537766, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '192216', 'license': 'MIT', 'n_tokens_mistral': 3926, 'n_tokens_neox': 3594, 'n_words': 653, 'path': 'notebooks/1-synthetic-experiments.ipynb', 'repo_name': 'trouleau/noisy-hawkes-cumulants', 'size': '9316'}