Update documentation

Author: FaDIn GitHub Action

_downloads/07fcc19ba03226cd3d83d4e40ec44385/auto_examples_python.zip

@@ -33,7 +33,7 @@
       },
       "outputs": [],
       "source": [
-        "import numpy as np\nimport matplotlib.pyplot as plt\n\nfrom fadin.utils.utils_simu import simu_marked_hawkes_cluster, custom_density\nfrom fadin.utils.utils_simu import simu_multi_poisson\nfrom fadin.solver import UNHaP\nfrom fadin.utils.functions import identity, linear_zero_one\nfrom fadin.utils.functions import reverse_linear_zero_one, truncated_gaussian\nfrom fadin.utils.vis import plot"
+        "import numpy as np\nimport matplotlib.pyplot as plt\n\nfrom fadin.utils.utils_simu import simulate_marked_data\nfrom fadin.solver import UNHaP\nfrom fadin.utils.vis import plot"
       ]
     },
     {
@@ -44,7 +44,7 @@
       },
       "outputs": [],
       "source": [
-        "baseline = np.array([0.3])\nbaseline_noise = np.array([0.05])\nalpha = np.array([[1.45]])\nmu = np.array([[0.4]])\nsigma = np.array([[0.1]])\n\ndelta = 0.01\nend_time = 10000\nseed = 0\nmax_iter = 20000\nbatch_rho = 200"
+        "baseline = np.array([0.3])\nbaseline_noise = np.array([0.05])\nalpha = np.array([[1.45]])\nmu = np.array([[0.4]])\nsigma = np.array([[0.1]])\n\ndelta = 0.01\nend_time = 1000\nseed = 0\nmax_iter = 2000\nbatch_rho = 200"
       ]
     },
     {
@@ -55,7 +55,7 @@
       },
       "outputs": [],
       "source": [
-        "def simulate_data(baseline, baseline_noise, alpha, end_time, seed=0):\n    n_dim = len(baseline)\n\n    marks_kernel = identity\n    marks_density = linear_zero_one\n    time_kernel = truncated_gaussian\n\n    params_marks_density = dict()\n    # params_marks_density = dict(scale=1)\n    params_marks_kernel = dict(slope=1.2)\n    params_time_kernel = dict(mu=mu, sigma=sigma)\n\n    marked_events, _ = simu_marked_hawkes_cluster(\n        end_time,\n        baseline,\n        alpha,\n        time_kernel,\n        marks_kernel,\n        marks_density,\n        params_marks_kernel=params_marks_kernel,\n        params_marks_density=params_marks_density,\n        time_kernel_length=None,\n        marks_kernel_length=None,\n        params_time_kernel=params_time_kernel,\n        random_state=seed,\n    )\n\n    noisy_events_ = simu_multi_poisson(end_time, [baseline_noise])\n\n    random_marks = [np.random.rand(noisy_events_[i].shape[0]) for i in range(n_dim)]\n    noisy_marks = [\n        custom_density(\n            reverse_linear_zero_one,\n            dict(),\n            size=noisy_events_[i].shape[0],\n            kernel_length=1.0,\n        )\n        for i in range(n_dim)\n    ]\n    noisy_events = [\n        np.concatenate(\n            (noisy_events_[i].reshape(-1, 1), random_marks[i].reshape(-1, 1)), axis=1\n        )\n        for i in range(n_dim)\n    ]\n\n    events = [\n        np.concatenate((noisy_events[i], marked_events[i]), axis=0)\n        for i in range(n_dim)\n    ]\n\n    events_cat = [events[i][events[i][:, 0].argsort()] for i in range(n_dim)]\n\n    labels = [\n        np.zeros(marked_events[i].shape[0] + noisy_events_[i].shape[0])\n        for i in range(n_dim)\n    ]\n    labels[0][-marked_events[0].shape[0] :] = 1.0\n    true_rho = [labels[i][events[i][:, 0].argsort()] for i in range(n_dim)]\n    # put the mark to one to test the impact of the marks\n    # events_cat[0][:, 1] = 1.\n\n    return events_cat, noisy_marks, true_rho\n\n\nev, noisy_marks, true_rho = simulate_data(\n    baseline, baseline_noise.item(), alpha, end_time, seed=0\n)"
+        "ev, noisy_marks, true_rho = simulate_marked_data(\n    baseline, baseline_noise.item(), alpha, end_time, mu, sigma, seed=0\n)"
       ]
     },
     {
@@ -66,7 +66,7 @@
       },
       "outputs": [],
       "source": [
-        "solver = UNHaP(\n    n_dim=1,\n    kernel=\"truncated_gaussian\",\n    kernel_length=1.0,\n    delta=delta,\n    optim=\"RMSprop\",\n    params_optim={\"lr\": 1e-3},\n    max_iter=max_iter,\n    batch_rho=batch_rho,\n    density_hawkes=\"linear\",\n    density_noise=\"uniform\",\n    moment_matching=True,\n)\nsolver.fit(ev, end_time)"
+        "print('Type of ev object', type(ev))\n# ev is a list of numpy arrays, one for each dimension\nprint('Number of events', len(ev[0]))\nprint('Shape of first event array', ev[0].shape)\n# Each dimension is stored as a numpy array of shape (n_events, 2).\nprint('First 10 events timestamps and marks', ev[0][:10])\n# Each event is stored as [timestamp, mark].\n# This is the expected data format for UNHaP.\nprint('First event timestamp', ev[0][0][0])\nprint('First event mark', ev[0][0][1])\nprint('Second event timestamp', ev[0][1][0])\nprint('Second event mark', ev[0][1][1])"
       ]
     },
     {
@@ -77,7 +77,7 @@
       },
       "outputs": [],
       "source": [
-        "print(\"Estimated baseline is: \", solver.param_baseline[-10:].mean().item())\nprint(\"Estimated alpha is: \", solver.param_alpha[-10:].mean().item())\nprint(\"Estimated kernel mean is: \", (solver.param_kernel[0][-10:].mean().item()))\nprint(\"Estimated kernel sd is: \", solver.param_kernel[1][-10:].mean().item())\nprint(\"Estimated noise baseline is: \", solver.param_baseline_noise[-10:].mean().item())\n# error on params\nerror_baseline = (solver.param_baseline[-10:].mean().item() - baseline.item()) ** 2\nerror_baseline_noise = (\n    solver.param_baseline_noise[-10:].mean().item() - baseline_noise.item()\n) ** 2\nerror_alpha = (solver.param_alpha[-10:].mean().item() - alpha.item()) ** 2\nerror_mu = (solver.param_kernel[0][-10:].mean().item() - 0.5) ** 2\nerror_sigma = (solver.param_kernel[1][-10:].mean().item() - 0.1) ** 2\nsum_error = error_baseline + error_baseline_noise + error_alpha + error_mu + error_sigma\nerror_params = np.sqrt(sum_error)\n\nprint(\"L2 square errors of the vector of parameters is:\", error_params)"
+        "solver = UNHaP(\n    n_dim=1,\n    kernel=\"truncated_gaussian\",\n    kernel_length=1.0,\n    init='moment_matching_mean',\n    delta=delta,\n    optim=\"RMSprop\",\n    params_optim={\"lr\": 1e-3},\n    max_iter=max_iter,\n    batch_rho=batch_rho,\n    density_hawkes=\"linear\",\n    density_noise=\"uniform\",\n)\nsolver.fit(ev, end_time)"
       ]
     },
     {
@@ -88,7 +88,18 @@
       },
       "outputs": [],
       "source": [
-        "fig, axs = plot(solver, plotfig=False, bl_noise=True, title=\"UNHaP fit\", savefig=None)\nplt.show(block=True)"
+        "print(\"Estimated baseline is: \", solver.baseline_.item())\nprint(\"Estimated alpha is: \", solver.alpha_.item())\nprint(\"Estimated kernel mean is: \", solver.kernel_[0].item())\nprint(\"Estimated kernel sd is: \", solver.kernel_[1].item())\nprint(\"Estimated noise baseline is: \", solver.baseline_noise_.item())\n# error on params\nerror_bl = (solver.baseline_.item() - baseline.item()) ** 2\nerror_bl_noise = (\n    solver.baseline_noise_.item() - baseline_noise.item()\n) ** 2\nerror_alpha = (solver.alpha_.item() - alpha.item()) ** 2\nerror_mu = (solver.kernel_[0].item() - mu.item()) ** 2\nerror_sigma = (solver.kernel_[1].item() - sigma.item()) ** 2\nsum_error = error_bl + error_bl_noise + error_alpha + error_mu + error_sigma\nerror_params = np.sqrt(sum_error)\n\nprint(\"L2 square error of the vector of parameters is:\", error_params)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "fig, axs = plot(\n    solver,\n    plotfig=False,\n    bl_noise=True,\n    title=\"UNHaP fit\",\n    savefig=None\n)\nplt.show(block=True)"
       ]
     }
   ],

_downloads/1c3b32a31630f6e84287c7300eec1f0b/plot_unhap.ipynb

@@ -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):

_downloads/249f41bcacce2f11424e25d9f61ef920/plot_multivariate_fadin.py

@@ -37,8 +37,8 @@
 # Here, we set the parameters of a Hawkes process with an Exponential(1) distribution.
 
 baseline = np.array([.4])
-alpha = np.array([[0.8]])
-beta = 2.
+alpha = np.array([[0.9]])
+beta = 0.8
 ###############################################################################
 # Here, we simulate the data.
 
@@ -48,6 +48,17 @@
 events = simu_hawkes_cluster(T, baseline, alpha, kernel,
                              params_kernel={'scale': 1 / beta})
 
+###############################################################################
+# Let's take a closer look at the events
+print('Type of events object', type(events))
+# events is a list of numpy arrays, one for each dimension
+print('Number of events', len(events[0]))
+print('First 10 events timestamps', events[0][:10])
+# For each event, its occurence time (timestamp) is stored in the numpy array.
+
+# `events`` is a list. Its elements are numpy arrays containing the timestamps
+# of the events. This is the expected data format for FaDIn.
+
 ###############################################################################
 # Here, we apply FaDIn.
 
@@ -60,24 +71,20 @@
                )
 solver.fit(events, T)
 
-# We average on the 10 last values of the optimization.
-
-estimated_baseline = solver.param_baseline[-10:].mean().item()
-estimated_alpha = solver.param_alpha[-10:].mean().item()
-param_kernel = [solver.param_kernel[0][-10:].mean().item()]
+###############################################################################
+# Here, we print the estimated parameters of the Hawkes process.
 
-print('Estimated baseline is:', estimated_baseline)
-print('Estimated alpha is:', estimated_alpha)
-print('Estimated beta parameter of the exponential kernel is:', param_kernel[0])
+print('Estimated baseline is:', solver.baseline_.item())
+print('Estimated alpha is:', solver.alpha_.item())
+print('Estimated beta parameter of the exponential kernel is:', solver.kernel_)
 
 
 ###############################################################################
 # Here, we plot the values of the estimated kernel with FaDIn.
 
 kernel = DiscreteKernelFiniteSupport(dt, n_dim, kernel='truncated_exponential',
                                      kernel_length=kernel_length)
-kernel_values = kernel.kernel_eval([torch.Tensor([param_kernel])],
-                                   discretization)
+kernel_values = kernel.kernel_eval(solver.kernel_, discretization)
 
 plt.plot(discretization[1:], kernel_values.squeeze()[1:]/kernel_length,
          label='FaDIn\' estimated kernel')

_downloads/25f9018d86375d7894bb79acc78a8375/plot_univariate_fadin.py

@@ -20,15 +20,12 @@
 import numpy as np
 import matplotlib.pyplot as plt
 
-from fadin.utils.utils_simu import simu_marked_hawkes_cluster, custom_density
-from fadin.utils.utils_simu import simu_multi_poisson
+from fadin.utils.utils_simu import simulate_marked_data
 from fadin.solver import UNHaP
-from fadin.utils.functions import identity, linear_zero_one
-from fadin.utils.functions import reverse_linear_zero_one, truncated_gaussian
 from fadin.utils.vis import plot
 
 
-# %% Fixing the parameter of the simulation setting
+# %% Fix the simulation and solver parameters
 
 baseline = np.array([0.3])
 baseline_noise = np.array([0.05])
@@ -37,120 +34,74 @@
 sigma = np.array([[0.1]])
 
 delta = 0.01
-end_time = 10000
+end_time = 1000
 seed = 0
-max_iter = 20000
+max_iter = 2000
 batch_rho = 200
 
-# %% Create the simulating function
-
-
-def simulate_data(baseline, baseline_noise, alpha, end_time, seed=0):
-    n_dim = len(baseline)
-
-    marks_kernel = identity
-    marks_density = linear_zero_one
-    time_kernel = truncated_gaussian
-
-    params_marks_density = dict()
-    # params_marks_density = dict(scale=1)
-    params_marks_kernel = dict(slope=1.2)
-    params_time_kernel = dict(mu=mu, sigma=sigma)
-
-    marked_events, _ = simu_marked_hawkes_cluster(
-        end_time,
-        baseline,
-        alpha,
-        time_kernel,
-        marks_kernel,
-        marks_density,
-        params_marks_kernel=params_marks_kernel,
-        params_marks_density=params_marks_density,
-        time_kernel_length=None,
-        marks_kernel_length=None,
-        params_time_kernel=params_time_kernel,
-        random_state=seed,
-    )
-
-    noisy_events_ = simu_multi_poisson(end_time, [baseline_noise])
-
-    random_marks = [np.random.rand(noisy_events_[i].shape[0]) for i in range(n_dim)]
-    noisy_marks = [
-        custom_density(
-            reverse_linear_zero_one,
-            dict(),
-            size=noisy_events_[i].shape[0],
-            kernel_length=1.0,
-        )
-        for i in range(n_dim)
-    ]
-    noisy_events = [
-        np.concatenate(
-            (noisy_events_[i].reshape(-1, 1), random_marks[i].reshape(-1, 1)), axis=1
-        )
-        for i in range(n_dim)
-    ]
-
-    events = [
-        np.concatenate((noisy_events[i], marked_events[i]), axis=0)
-        for i in range(n_dim)
-    ]
-
-    events_cat = [events[i][events[i][:, 0].argsort()] for i in range(n_dim)]
-
-    labels = [
-        np.zeros(marked_events[i].shape[0] + noisy_events_[i].shape[0])
-        for i in range(n_dim)
-    ]
-    labels[0][-marked_events[0].shape[0] :] = 1.0
-    true_rho = [labels[i][events[i][:, 0].argsort()] for i in range(n_dim)]
-    # put the mark to one to test the impact of the marks
-    # events_cat[0][:, 1] = 1.
-
-    return events_cat, noisy_marks, true_rho
-
-
-ev, noisy_marks, true_rho = simulate_data(
-    baseline, baseline_noise.item(), alpha, end_time, seed=0
+# %% Simulate Hawkes Process with truncated Gaussian kernel and Poisson noise
+
+ev, noisy_marks, true_rho = simulate_marked_data(
+    baseline, baseline_noise.item(), alpha, end_time, mu, sigma, seed=0
 )
-# %% Apply UNHAP
+
+# %% Let's take a closer look at the events
+print('Type of ev object', type(ev))
+# ev is a list of numpy arrays, one for each dimension
+print('Number of events', len(ev[0]))
+print('Shape of first event array', ev[0].shape)
+# Each dimension is stored as a numpy array of shape (n_events, 2).
+print('First 10 events timestamps and marks', ev[0][:10])
+# Each event is stored as [timestamp, mark].
+# This is the expected data format for UNHaP.
+print('First event timestamp', ev[0][0][0])
+print('First event mark', ev[0][0][1])
+print('Second event timestamp', ev[0][1][0])
+print('Second event mark', ev[0][1][1])
+
+# %% Initiate and fit UNHAP to the simulated events
 
 solver = UNHaP(
     n_dim=1,
     kernel="truncated_gaussian",
     kernel_length=1.0,
+    init='moment_matching_mean',
     delta=delta,
     optim="RMSprop",
     params_optim={"lr": 1e-3},
     max_iter=max_iter,
     batch_rho=batch_rho,
     density_hawkes="linear",
     density_noise="uniform",
-    moment_matching=True,
 )
 solver.fit(ev, end_time)
 
 # %% Print estimated parameters
 
-print("Estimated baseline is: ", solver.param_baseline[-10:].mean().item())
-print("Estimated alpha is: ", solver.param_alpha[-10:].mean().item())
-print("Estimated kernel mean is: ", (solver.param_kernel[0][-10:].mean().item()))
-print("Estimated kernel sd is: ", solver.param_kernel[1][-10:].mean().item())
-print("Estimated noise baseline is: ", solver.param_baseline_noise[-10:].mean().item())
+print("Estimated baseline is: ", solver.baseline_.item())
+print("Estimated alpha is: ", solver.alpha_.item())
+print("Estimated kernel mean is: ", solver.kernel_[0].item())
+print("Estimated kernel sd is: ", solver.kernel_[1].item())
+print("Estimated noise baseline is: ", solver.baseline_noise_.item())
 # error on params
-error_baseline = (solver.param_baseline[-10:].mean().item() - baseline.item()) ** 2
-error_baseline_noise = (
-    solver.param_baseline_noise[-10:].mean().item() - baseline_noise.item()
+error_bl = (solver.baseline_.item() - baseline.item()) ** 2
+error_bl_noise = (
+    solver.baseline_noise_.item() - baseline_noise.item()
 ) ** 2
-error_alpha = (solver.param_alpha[-10:].mean().item() - alpha.item()) ** 2
-error_mu = (solver.param_kernel[0][-10:].mean().item() - 0.5) ** 2
-error_sigma = (solver.param_kernel[1][-10:].mean().item() - 0.1) ** 2
-sum_error = error_baseline + error_baseline_noise + error_alpha + error_mu + error_sigma
+error_alpha = (solver.alpha_.item() - alpha.item()) ** 2
+error_mu = (solver.kernel_[0].item() - mu.item()) ** 2
+error_sigma = (solver.kernel_[1].item() - sigma.item()) ** 2
+sum_error = error_bl + error_bl_noise + error_alpha + error_mu + error_sigma
 error_params = np.sqrt(sum_error)
 
-print("L2 square errors of the vector of parameters is:", error_params)
+print("L2 square error of the vector of parameters is:", error_params)
 
 # %% Plot estimated parameters
-fig, axs = plot(solver, plotfig=False, bl_noise=True, title="UNHaP fit", savefig=None)
+fig, axs = plot(
+    solver,
+    plotfig=False,
+    bl_noise=True,
+    title="UNHaP fit",
+    savefig=None
+)
 plt.show(block=True)
-# %%