Add UNHAP solver

examples/plot_multivariate_fadin.py

@@ -35,7 +35,8 @@
 discretization = torch.linspace(0, kernel_length, L)
 
 ###############################################################################
-# Here, we set the parameters of a Hawkes process with a Exponential(1) distributions.
+# Here, we set the parameters of a Hawkes process with a Exponential(1)
+# distribution.
 
 baseline = np.array([.1, .5])
 alpha = np.array([[0.6, 0.3], [0.25, 0.7]])
@@ -50,13 +51,14 @@
 ###############################################################################
 # Here, we apply FaDIn.
 
-solver = FaDIn(n_dim=n_dim,
-               kernel="truncated_exponential",
-               kernel_length=kernel_length,
-               delta=dt, optim="RMSprop",
-               params_optim={'lr': 1e-3},
-               max_iter=10000, criterion='l2'
-               )
+solver = FaDIn(
+    n_dim=n_dim,
+    kernel="truncated_exponential",
+    kernel_length=kernel_length,
+    delta=dt, optim="RMSprop",
+    params_optim={'lr': 1e-3},
+    max_iter=10000
+)
 solver.fit(events, T)
 
 # We average on the 10 last values of the optimization.

examples/plot_unhap.py

@@ -0,0 +1,108 @@
+# %% import stuff
+import numpy as np
+
+from fadin.utils.utils_simu import simu_marked_hawkes_cluster, custom_density, \
+    simu_multi_poisson
+from fadin.solver import UNHaP
+from fadin.utils.functions import identity, linear_zero_one, reverse_linear_zero_one, \
+    truncated_gaussian
+
+
+# %% Fixing the parameter of the simulation setting
+
+baseline = np.array([.3])
+baseline_noise = np.array([.05])
+alpha = np.array([[1.45]])
+mu = np.array([[0.5]])
+sigma = np.array([[0.1]])
+
+delta = 0.01
+end_time = 10000
+seed = 0
+max_iter = 20000
+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.) for i in range(n_dim)]
+    noisy_events = [
+        np.concatenate((noisy_events_[i].reshape(-1, 1),
+                        noisy_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.
+    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)
+# %% Apply UNHAP
+
+
+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=batch_rho,
+               density_hawkes='linear',
+               density_noise='reverse_linear',
+               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())
+
+# error on params
+error_baseline = (solver.param_baseline[-10:].mean().item() - baseline.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_alpha + error_mu + error_sigma
+error_params = np.sqrt(sum_error)
+
+print('L2 square errors of the vector of parameters is:', error_params)
+# %%

examples/plot_univariate_fadin.py

@@ -56,7 +56,7 @@
                kernel_length=kernel_length,
                delta=dt, optim="RMSprop",
                params_optim={'lr': 1e-3},
-               max_iter=2000, criterion='l2'
+               max_iter=2000
                )
 solver.fit(events, T)
 

experiments/benchmark/run_benchmark_exp.py

@@ -25,8 +25,13 @@ def simulate_data(baseline, alpha, decay, T, dt, seed=0):
     L = int(1 / dt)
     discretization = torch.linspace(0, 1, L)
     n_dim = decay.shape[0]
-    EXP = DiscreteKernelFiniteSupport(dt, n_dim=n_dim, kernel='truncated_exponential',
-                                      lower=0, upper=1)
+    EXP = DiscreteKernelFiniteSupport(
+        dt,
+        n_dim=n_dim,
+        kernel='truncated_exponential',
+        lower=0,
+        upper=1
+    )
 
     kernel_values = EXP.kernel_eval([torch.Tensor(decay)],
                                     discretization)
@@ -62,14 +67,16 @@ def simulate_data(baseline, alpha, decay, T, dt, seed=0):
 def run_fadin(events, decay_init, baseline_init, alpha_init, T, dt, seed=0):
     start = time.time()
     max_iter = 2000
+    init = {
+        'alpha': torch.tensor(alpha_init),
+        'baseline': torch.tensor(baseline_init),
+        'kernel': [torch.tensor(decay_init)]
+    }
     solver = FaDIn(2,
                    "truncated_exponential",
-                   [torch.tensor(decay_init)],
-                   torch.tensor(baseline_init),
-                   torch.tensor(alpha_init),
+                   init=init,
                    delta=dt, optim="RMSprop",
                    step_size=1e-3, max_iter=max_iter,
-                   optimize_kernel=True, precomputations=True,
                    ztzG_approx=True, device='cpu', log=False
                    )
 
@@ -258,14 +265,17 @@ def run_experiment(baseline, alpha, decay, T, dt, seed=0):
     alpha_hat = results['param_alpha']
     decay_hat = results['param_kernel'][0]
 
-    RC = DiscreteKernelFiniteSupport(dt, n_dim=2, kernel='truncated_exponential',
+    RC = DiscreteKernelFiniteSupport(dt, n_dim=2,
+                                     kernel='truncated_exponential',
                                      lower=0, upper=1)
     intens_fadin = RC.intensity_eval(torch.tensor(baseline_hat),
                                      torch.tensor(alpha_hat),
                                      [torch.Tensor(decay_hat)],
                                      events_grid, torch.linspace(0, 1, L))
 
-    res['err_fadin'] = np.absolute(intens.numpy() - intens_fadin.numpy()).mean()
+    res['err_fadin'] = np.absolute(
+        intens.numpy() - intens_fadin.numpy()
+    ).mean()
     res['time_fadin'] = results['time']
 
     results = run_gibbs(S, size_grid, dt, seed=seed)

experiments/benchmark/run_benchmark_rc.py

@@ -63,15 +63,16 @@ def simulate_data(baseline, alpha, mu, sigma, T, dt, seed=0):
 def run_fadin(events, u_init, sigma_init, baseline_init, alpha_init, T, dt, seed=0):
     start = time.time()
     max_iter = 2000
+    init = {
+        'kernel': [torch.tensor(u_init), torch.tensor(sigma_init)],
+        'baseline': torch.tensor(baseline_init),
+        'alpha': torch.tensor(alpha_init)
+    }
     solver = FaDIn(2,
                    "raised_cosine",
-                   [torch.tensor(u_init),
-                    torch.tensor(sigma_init)],
-                   torch.tensor(baseline_init),
-                   torch.tensor(alpha_init),
+                   init=init,
                    delta=dt, optim="RMSprop",
                    step_size=1e-3, max_iter=max_iter,
-                   optimize_kernel=True, precomputations=True,
                    ztzG_approx=True, device='cpu', log=False
                    )
 

experiments/benchmark/run_benchmark_tg.py

@@ -62,15 +62,16 @@ def simulate_data(baseline, alpha, m, sigma, T, dt, seed=0):
 def run_fadin(events, m_init, sigma_init, baseline_init, alpha_init, T, dt, seed=0):
     start = time.time()
     max_iter = 2000
+    init = {
+        'alpha': torch.tensor(alpha_init),
+        'baseline': torch.tensor(baseline_init),
+        'kernel': [torch.tensor(m_init), torch.tensor(sigma_init)]
+    }
     solver = FaDIn(2,
                    "truncated_gaussian",
-                   [torch.tensor(m_init),
-                    torch.tensor(sigma_init)],
-                   torch.tensor(baseline_init),
-                   torch.tensor(alpha_init),
+                   init=init,
                    delta=dt, optim="RMSprop",
                    step_size=1e-3, max_iter=max_iter,
-                   optimize_kernel=True, precomputations=True,
                    ztzG_approx=True, device='cpu', log=False
                    )
 

experiments/benchmark/run_comparison_ll.py

@@ -11,7 +11,55 @@
 from fadin.kernels import DiscreteKernelFiniteSupport
 from fadin.solver import FaDIn
 
-# %% simulate data
+# %%
+################################
+# Define solver with loglikelihood criterion
+################################
+
+
+def discrete_ll_loss_conv(intensity, events_grid, delta):
+    """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() / events_grid.sum()
+
+
+def compute_gradient_loglikelihood(solver, events_grid, discretization,
+                                   i, n_events, end_time):
+    """One optimizer iteration of FaDIn_loglikelihood solver,
+    with loglikelihood loss.
+    """
+    intens = solver.kernel_model.intensity_eval(
+        solver.params_intens[0],
+        solver.params_intens[1],
+        solver.params_intens[2:],
+        events_grid,
+        discretization
+    )
+    loss = discrete_ll_loss_conv(intens, events_grid, solver.delta)
+    loss.backward()
+
+
+class FaDInLogLikelihood(FaDIn):
+    """Define the FaDIn framework for estimated Hawkes processes *with
+    loglikelihood criterion instead of l2 loss*."""
+    compute_gradient = staticmethod(compute_gradient_loglikelihood)
+    precomputations = False
+
+################################
 # Simulated data
 ################################
 
@@ -56,23 +104,40 @@ def simulate_data(baseline, alpha, kernel_params,
 
 
 def run_solver(criterion, events, kernel_params_init,
-               baseline_init, alpha_init, T, dt, seed=0, kernel='raised_cosine'):
-    k_params_init = [torch.tensor(a) for a in kernel_params_init]
+               baseline_init, alpha_init, T, dt, seed=0,
+               kernel='raised_cosine'):
     max_iter = 2000
-    solver = FaDIn(1,
-                   kernel,
-                   k_params_init,
-                   torch.tensor(baseline_init),
-                   torch.tensor(alpha_init),
-                   delta=dt,
-                   optim="RMSprop",
-                   step_size=1e-3,
-                   max_iter=max_iter,
-                   log=False,
-                   random_state=seed,
-                   device="cpu",
-                   optimize_kernel=True,
-                   criterion=criterion)
+    init = {
+        'alpha': torch.tensor(alpha_init),
+        'baseline': torch.tensor(baseline_init),
+        'kernel': [torch.tensor(a) for a in kernel_params_init]
+    }
+    if criterion == 'l2':
+        solver = FaDIn(
+            1,
+            kernel,
+            init=init,
+            delta=dt,
+            optim="RMSprop",
+            step_size=1e-3,
+            max_iter=max_iter,
+            log=False,
+            random_state=seed,
+            device="cpu"
+        )
+    elif criterion == 'll':
+        solver = FaDInLogLikelihood(
+            1,
+            kernel,
+            init=init,
+            delta=dt,
+            optim="RMSprop",
+            step_size=1e-3,
+            max_iter=max_iter,
+            log=False,
+            random_state=seed,
+            device="cpu"
+        )
     results = solver.fit(events, T)
     if kernel == 'truncated_exponential':
         results_ = dict(param_baseline=results['param_baseline'][-10:].mean().item(),

experiments/benchmark/run_nonparam_exp.py

@@ -61,19 +61,22 @@ def simulate_data(baseline, alpha, decay, T, dt, seed=0):
 #  @mem.cache
 def run_solver(events, decay_init, baseline_init, alpha_init, T, dt, seed=0):
     max_iter = 800
+    init = {
+        'alpha': torch.tensor(alpha_init),
+        'baseline': torch.tensor(baseline_init),
+        'kernel': [torch.tensor(decay_init)]
+    }
     solver = FaDIn(1,
                    "truncated_exponential",
-                   [torch.tensor(decay_init)],
-                   torch.tensor(baseline_init),
-                   torch.tensor(alpha_init),
+                   init=init,
                    delta=dt,
                    optim="RMSprop",
                    step_size=1e-3,
                    max_iter=max_iter,
                    log=False,
                    random_state=seed,
                    device="cpu",
-                   optimize_kernel=True)
+                   )
     results = solver.fit(events, T)
     return results