Merge pull request #147 from X-DataInitiative/implement-nphc

Implement Hawkes Cumulant non parametric learner

doc/modules/api.rst

@@ -50,6 +50,7 @@ Learners
    HawkesBasisKernels
    HawkesSumGaussians
    HawkesConditionalLaw
+   HawkesCumulantMatching
 
 Simulation
 ----------

doc/modules/hawkes.rst

@@ -96,6 +96,7 @@ These learners might be used to infer more exotic kernels.
    HawkesEM
    HawkesBasisKernels
    HawkesConditionalLaw
+   HawkesCumulantMatching
 
 Example
 *******

examples/plot_hawkes_cumulants_matching.py

@@ -0,0 +1,55 @@
+"""
+=======================================
+Fit Hawkes kernel norms using cumulants
+=======================================
+
+This non parametric Hawkes cumulants matching
+(`tick.hawkes.HawkesCumulantMatching`) algorithm estimates directly
+kernels norms without making any assumption on kernel shapes.
+
+It has been originally described in this paper:
+
+Achab, M., Bacry, E., Gaiffas, S., Mastromatteo, I., & Muzy, J. F.
+(2017, July). Uncovering causality from multivariate Hawkes integrated
+cumulants.
+`In International Conference on Machine Learning (pp. 1-10)`_.
+
+.. _In International Conference on Machine Learning (pp. 1-10): http://proceedings.mlr.press/v70/achab17a.html
+"""
+
+import numpy as np
+
+from tick.hawkes import (HawkesCumulantMatching, SimuHawkesExpKernels,
+                         SimuHawkesMulti)
+from tick.plot import plot_hawkes_kernel_norms
+
+np.random.seed(7168)
+
+n_nodes = 3
+baselines = 0.3 * np.ones(n_nodes)
+decays = 0.5 + np.random.rand(n_nodes, n_nodes)
+adjacency = np.array([
+    [1, 1, -0.5],
+    [0, 1, 0],
+    [0, 0, 2],
+], dtype=float)
+
+adjacency /= 4
+
+end_time = 1e5
+integration_support = 5
+n_realizations = 5
+
+simu_hawkes = SimuHawkesExpKernels(
+    baseline=baselines, adjacency=adjacency, decays=decays,
+    end_time=end_time, verbose=False, seed=7168)
+simu_hawkes.threshold_negative_intensity(True)
+
+multi = SimuHawkesMulti(simu_hawkes, n_simulations=n_realizations, n_threads=-1)
+multi.simulate()
+
+nphc = HawkesCumulantMatching(integration_support, cs_ratio=.15, tol=1e-10,
+                              step=0.3)
+
+nphc.fit(multi.timestamps)
+plot_hawkes_kernel_norms(nphc)

examples/plot_hawkes_em.py

@@ -7,7 +7,7 @@
 piecewise constant. Hence it can fit basically any kernel form. However it
 doesn't scale very well.
 
-I has been originally described in this paper:
+It has been originally described in this paper:
 
 Lewis, E., & Mohler, G. (2011).
 A nonparametric EM algorithm for multiscale Hawkes processes.

lib/cpp/hawkes/inference/CMakeLists.txt

@@ -6,11 +6,13 @@ add_library(tick_hawkes_inference EXCLUDE_FROM_ALL
         ${TICK_HAWKES_INFERENCE_INCLUDE_DIR}/hawkes_adm4.h
         ${TICK_HAWKES_INFERENCE_INCLUDE_DIR}/hawkes_basis_kernels.h
         ${TICK_HAWKES_INFERENCE_INCLUDE_DIR}/hawkes_sumgaussians.h
+        ${TICK_HAWKES_INFERENCE_INCLUDE_DIR}/hawkes_cumulant.h
         hawkes_adm4.cpp
         hawkes_basis_kernels.cpp
         hawkes_conditional_law.cpp
         hawkes_em.cpp
         hawkes_sumgaussians.cpp
+        hawkes_cumulant.cpp
         )
 
 target_link_libraries(tick_hawkes_inference

lib/cpp/hawkes/inference/hawkes_cumulant.cpp

@@ -0,0 +1,122 @@
+#include "tick/hawkes/inference/hawkes_cumulant.h"
+
+HawkesCumulant::HawkesCumulant(double integration_support)
+  : integration_support(integration_support), are_cumulants_ready(false) { }
+
+SArrayDoublePtr HawkesCumulant::compute_A_and_I_ij(ulong r, ulong i, ulong j,
+                                                   double mean_intensity_j) {
+  auto timestamps_i = timestamps_list[r][i];
+  auto timestamps_j = timestamps_list[r][j];
+
+  ulong n_i = timestamps_i->size();
+  ulong n_j = timestamps_j->size();
+  double res_C = 0;
+  double res_J = 0;
+  double width = 2 * integration_support;
+  double trend_C_j = mean_intensity_j * width;
+  double trend_J_j = mean_intensity_j * width * width;
+
+  ulong last_l = 0;
+  for (ulong k = 0; k < n_i; ++k) {
+    double t_i_k = (*timestamps_i)[k];
+    double t_i_k_minus_half_width = t_i_k - integration_support;
+    double t_i_k_minus_width = t_i_k - width;
+
+    if (t_i_k_minus_half_width < 0) continue;
+
+    // Find next t_j_l that occurs width before t_i_k
+    while (last_l < n_j) {
+      if ((*timestamps_j)[last_l] <= t_i_k_minus_width) ++last_l;
+      else
+        break;
+    }
+
+    ulong l = last_l;
+    ulong timestamps_in_interval = 0;
+
+    double sub_res = 0.;
+
+    while (l < n_j) {
+      double t_j_l_minus_t_i_k = (*timestamps_j)[l] - t_i_k;
+      double abs_t_j_l_minus_t_i_k = fabs(t_j_l_minus_t_i_k);
+
+      if (abs_t_j_l_minus_t_i_k < width) {
+        sub_res += width - abs_t_j_l_minus_t_i_k;
+
+        if (abs_t_j_l_minus_t_i_k < integration_support) timestamps_in_interval++;
+      } else {
+        break;
+      }
+      l += 1;
+    }
+
+    if (l == n_j) continue;
+    res_C += timestamps_in_interval - trend_C_j;
+    res_J += sub_res - trend_J_j;
+  }
+
+  res_C /= (*end_times)[r];
+  res_J /= (*end_times)[r];
+
+  ArrayDouble return_array{res_C, res_J};
+  return return_array.as_sarray_ptr();
+}
+
+double HawkesCumulant::compute_E_ijk(ulong r, ulong i, ulong j, ulong k,
+                                     double mean_intensity_i, double mean_intensity_j,
+                                     double J_ij) {
+  auto timestamps_i = timestamps_list[r][i];
+  auto timestamps_j = timestamps_list[r][j];
+  auto timestamps_k = timestamps_list[r][k];
+
+  double L_i = mean_intensity_i;
+  double L_j = mean_intensity_j;
+
+  double res = 0;
+  ulong last_l = 0;
+  ulong last_m = 0;
+  ulong n_i = timestamps_i->size();
+  ulong n_j = timestamps_j->size();
+  ulong n_k = timestamps_k->size();
+
+  double trend_i = L_i * 2 * integration_support;
+  double trend_j = L_j * 2 * integration_support;
+
+  for (ulong t = 0; t < n_k; ++t) {
+    double tau = (*timestamps_k)[t];
+
+    if (tau - integration_support < 0) continue;
+
+    while (last_l < n_i) {
+      if ((*timestamps_i)[last_l] <= tau - integration_support) last_l += 1;
+      else
+        break;
+    }
+    ulong l = last_l;
+
+    while (l < n_i) {
+      if ((*timestamps_i)[l] < tau + integration_support) l += 1;
+      else
+        break;
+    }
+
+    while (last_m < n_j) {
+      if ((*timestamps_j)[last_m] <= tau - integration_support) last_m += 1;
+      else
+        break;
+    }
+    ulong m = last_m;
+
+    while (m < n_j) {
+      if ((*timestamps_j)[m] < tau + integration_support) m += 1;
+      else
+        break;
+    }
+
+    if ((m == n_j) || (l == n_i)) continue;
+
+    res += (l - last_l - trend_i) * (m - last_m - trend_j) - J_ij;
+  }
+  res /= (*end_times)[r];
+  return res;
+}

lib/cpp/hawkes/model/base/model_hawkes_list.cpp

@@ -6,7 +6,7 @@
 ModelHawkesList::ModelHawkesList(
   const int max_n_threads,
   const unsigned int optimization_level)
-  : ModelHawkes(max_n_threads, optimization_level), n_realizations(0) {
+  : ModelHawkes(max_n_threads, optimization_level), n_realizations(0), timestamps_list(0) {
   n_jumps_per_realization = VArrayULong::new_ptr(n_realizations);
   end_times = VArrayDouble::new_ptr(n_realizations);
 }

lib/include/tick/hawkes/inference/hawkes_cumulant.h

@@ -0,0 +1,41 @@
+#ifndef LIB_INCLUDE_TICK_HAWKES_INFERENCE_HAWKES_CUMULANT_H_
+#define LIB_INCLUDE_TICK_HAWKES_INFERENCE_HAWKES_CUMULANT_H_
+
+// License: BSD 3 clause
+
+#include "tick/base/base.h"
+#include "tick/hawkes/model/base/model_hawkes_list.h"
+
+class HawkesCumulant : public ModelHawkesList {
+  double integration_support;
+  bool are_cumulants_ready;
+
+ public:
+  explicit HawkesCumulant(double integration_support);
+
+  SArrayDoublePtr compute_A_and_I_ij(ulong r, ulong i, ulong j, double mean_intensity_j);
+
+  double compute_E_ijk(ulong r, ulong i, ulong j, ulong k,
+                       double mean_intensity_i, double mean_intensity_j,
+                       double J_ij);
+
+  double get_integration_support() const {
+    return integration_support;
+  }
+
+  void set_integration_support(const double integration_support) {
+    if (integration_support <= 0) TICK_ERROR("Kernel support must be positive");
+    this->integration_support = integration_support;
+    are_cumulants_ready = false;
+  }
+
+  bool get_are_cumulants_ready() const {
+    return are_cumulants_ready;
+  }
+
+  void set_are_cumulants_ready(const bool are_cumulants_ready) {
+    this->are_cumulants_ready = are_cumulants_ready;
+  }
+};
+
+#endif  // LIB_INCLUDE_TICK_HAWKES_INFERENCE_HAWKES_CUMULANT_H_

lib/include/tick/hawkes/model/base/model_hawkes_list.h

@@ -20,7 +20,7 @@ class DLL_PUBLIC ModelHawkesList : public ModelHawkes {
   SArrayDoublePtrList2D timestamps_list;
 
   //! @brief Ending time of the realization
-  VArrayDoublePtr end_times;
+  VArrayDoublePtr end_times = nullptr;
 
   //! @brief Number of jumps of the process per realization (size=n_realizations)
   VArrayULongPtr n_jumps_per_realization;
@@ -33,7 +33,8 @@ class DLL_PUBLIC ModelHawkesList : public ModelHawkes {
   ModelHawkesList(const int max_n_threads = 1,
                   const unsigned int optimization_level = 0);
 
-  void set_data(const SArrayDoublePtrList2D &timestamps_list, const VArrayDoublePtr end_times);
+  virtual void set_data(const SArrayDoublePtrList2D &timestamps_list,
+                        const VArrayDoublePtr end_times);
 
   //! @brief returns the number of jumps per realization
   SArrayULongPtr get_n_jumps_per_realization() const {
@@ -46,6 +47,10 @@ class DLL_PUBLIC ModelHawkesList : public ModelHawkes {
 
   virtual unsigned int get_n_threads() const;
 
+  SArrayDoublePtrList2D get_timestamps_list() const {
+    return timestamps_list;
+  }
+
  public:
   template<class Archive>
   void serialize(Archive &ar) {

lib/swig/hawkes/inference/hawkes_cumulant.i

@@ -0,0 +1,26 @@
+
+
+%include std_shared_ptr.i
+%shared_ptr(HawkesCumulant);
+
+%{
+#include "tick/hawkes/inference/hawkes_cumulant.h"
+%}
+
+
+class HawkesCumulant : public ModelHawkesList {
+
+public:
+  HawkesCumulant(double integration_support);
+
+  SArrayDoublePtr compute_A_and_I_ij(ulong r, ulong i, ulong j, double mean_intensity_j);
+
+  double compute_E_ijk(ulong r, ulong i, ulong j, ulong k,
+                       double mean_intensity_i, double mean_intensity_j,
+                       double J_ij);
+
+  double get_integration_support() const;
+  void set_integration_support(const double integration_support);
+  bool get_are_cumulants_ready() const;
+  void set_are_cumulants_ready(const bool recompute_cumulants);
+};

lib/swig/hawkes/inference/hawkes_inference_module.i

@@ -23,3 +23,4 @@
 %include hawkes_adm4.i
 %include hawkes_basis_kernels.i
 %include hawkes_sumgaussians.i
+%include hawkes_cumulant.i

lib/swig/hawkes/model/base/model_hawkes_list.i

@@ -18,5 +18,7 @@ class ModelHawkesList : public ModelHawkes {
   ulong get_n_threads() const;
   SArrayULongPtr get_n_jumps_per_realization() const;
 
+  SArrayDoublePtrList2D get_timestamps_list() const;
+
   void set_n_threads(const int max_n_threads);
 };

tick/base/base.py

@@ -249,7 +249,8 @@ def find_documented_attributes(class_name, attrs):
 
         current_class_doc = inspect.cleandoc(attrs['__doc__'])
         parsed_doc = nd.docscrape.ClassDoc(None, doc=current_class_doc)
-        attr_docs = parsed_doc['Parameters'] + parsed_doc['Attributes']
+        attr_docs = parsed_doc['Parameters'] + parsed_doc['Attributes'] + \
+            parsed_doc['Other Parameters']
 
         attr_and_doc = []
 

tick/hawkes/__init__.py

@@ -12,7 +12,7 @@
 )
 from .inference import (
     HawkesADM4, HawkesExpKern, HawkesSumExpKern, HawkesBasisKernels,
-    HawkesConditionalLaw, HawkesEM, HawkesSumGaussians
+    HawkesConditionalLaw, HawkesEM, HawkesSumGaussians, HawkesCumulantMatching
 )
 
 __all__ = [
@@ -23,5 +23,5 @@
     "SimuPoissonProcess", "SimuInhomogeneousPoisson", "SimuHawkes",
     "SimuHawkesMulti", "SimuHawkesExpKernels", "SimuHawkesSumExpKernels",
     "HawkesKernel0", "HawkesKernelExp", "HawkesKernelPowerLaw",
-    "HawkesKernelSumExp", "HawkesKernelTimeFunc"
+    "HawkesKernelSumExp", "HawkesKernelTimeFunc", "HawkesCumulantMatching",
 ]

tick/hawkes/inference/__init__.py

@@ -7,6 +7,7 @@
 from .hawkes_expkern_fixeddecay import HawkesExpKern
 from .hawkes_sumexpkern_fixeddecay import HawkesSumExpKern
 from .hawkes_sumgaussians import HawkesSumGaussians
+from .hawkes_cumulant_matching import HawkesCumulantMatching
 
 __all__ = [
     "HawkesExpKern",
@@ -15,5 +16,6 @@
     "HawkesEM",
     "HawkesADM4",
     "HawkesBasisKernels",
-    "HawkesSumGaussians"
+    "HawkesSumGaussians",
+    "HawkesCumulantMatching",
 ]