Merge pull request #123 from X-DataInitiative/simulate-pp-threshold-i…

…ntensity

Allow simulation of point processes with thresholded negative intensi…

tick/simulation/base/simu_point_process.py

@@ -183,3 +183,14 @@ def reset(self):
         """Reset the process, so that is is ready for a brand new simulation
         """
         self._pp.reset()
+
+    def threshold_negative_intensity(self, allow=True):
+        """Threshold intensity to 0 if it becomes negative.
+        This allows simulation with negative kernels
+
+        Parameters
+        ----------
+        allow : `bool`
+            Flag to allow negative intensity thresholding
+        """
+        self._pp.set_threshold_negative_intensity(allow)

tick/simulation/src/hawkes.cpp

@@ -45,7 +45,10 @@ bool Hawkes::update_time_shift_(double delay,
       if (total_intensity_bound1) {
         *total_intensity_bound1 += bound;
       }
-      if (intensity[i] < 0) flag_negative_intensity1 = true;
+      if (intensity[i] < 0) {
+        if (threshold_negative_intensity) intensity[i] = 0;
+        flag_negative_intensity1 = true;
+      }
     }
   }
   return flag_negative_intensity1;

tick/simulation/src/hawkes_kernels/hawkes_kernel_sum_exp.cpp

@@ -106,10 +106,13 @@ double HawkesKernelSumExp::get_convolution(const double time, const ArrayDouble
 
   if (bound) {
     if (!intensities_all_positive) {
-      TICK_ERROR("All intensities of HawkesKernelSumExp must be positive "
-                     "to compute kernel bound");
+      *bound = 0;
+      for (ulong u = 0; u < n_decays; ++u) {
+        if (intensities[u] > 0) *bound += last_convolution_values[u];
+      }
+    } else {
+      *bound = value;
     }
-    *bound = value;
   }
 
   return value;

tick/simulation/src/pp.cpp

@@ -116,31 +116,32 @@ void PP::simulate(double end_time, ulong n_points) {
   }
 
   // We loop till we reach the endTime
-  while (time < end_time && n_total_jumps < n_points && !flag_negative_intensity) {
+  while (time < end_time && n_total_jumps < n_points &&
+    (!flag_negative_intensity || threshold_negative_intensity) ) {
     // We compute the time of the potential next random jump
-    const double timeOfNextJump = time + rand.exponential(total_intensity_bound);
+    const double time_of_next_jump = time + rand.exponential(total_intensity_bound);
 
     // If we must track record the intensities we perform a loop
     if (itr_on()) {
-      while (itr_time + itr_time_step < std::min(timeOfNextJump, end_time)) {
+      while (itr_time + itr_time_step < std::min(time_of_next_jump, end_time)) {
         update_time_shift(itr_time_step + itr_time - time, false, true);
-        if (flag_negative_intensity) break;
+        if (flag_negative_intensity && !threshold_negative_intensity) break;
         itr_time = itr_time + itr_time_step;
       }
-      if (flag_negative_intensity) break;
+      if (flag_negative_intensity && !threshold_negative_intensity) break;
     }
 
     // Are we done ?
-    if (timeOfNextJump >= end_time) {
+    if (time_of_next_jump >= end_time) {
       time = end_time;
       break;
     }
 
     // We go to timeOfNextJump
     // Do not compute intensity bound here as we want new intensities but old bound that we
     // have used for the exponential law
-    update_time_shift(timeOfNextJump - time, false, false);
-    if (flag_negative_intensity) break;
+    update_time_shift(time_of_next_jump - time, false, false);
+    if (flag_negative_intensity && !threshold_negative_intensity) break;
 
     // We need to understand which component jumps
     double temp = rand.uniform() * total_intensity_bound;
@@ -154,7 +155,7 @@ void PP::simulate(double end_time, ulong n_points) {
     // Case we discard the jump, we should recompute max intensity ?????????  !
     if (i == n_nodes) {
       update_time_shift(0, true, false);
-      if (flag_negative_intensity) break;
+      if (flag_negative_intensity && !threshold_negative_intensity) break;
       continue;
     }
 
@@ -164,16 +165,17 @@ void PP::simulate(double end_time, ulong n_points) {
     // We compute the new intensities (taking into account eventually the new jump)
     update_time_shift(0, true, true);
 
-    if (flag_negative_intensity) break;
+    if (flag_negative_intensity && !threshold_negative_intensity) break;
   }
 
   // This causes deadlock, see MLPP-334 - Investigate deadlock in PP
   // #ifdef PYTHON_LINK
   //    Py_END_ALLOW_THREADS;
   // #endif
 
-  if (flag_negative_intensity) TICK_ERROR(
-      "Stopped because intensity went negative (you could set the field ``thresholdNegativeIntensity`` to True)");
+  if (flag_negative_intensity && !threshold_negative_intensity) TICK_ERROR(
+      "Simulation stopped because intensity went negative (you could call "
+        "``threshold_negative_intensity`` to allow it)");
 }
 
 // Update the process component 'index' with current time

tick/simulation/src/pp.h

@@ -60,9 +60,9 @@ class PP {
   // Keeps track of maximum total intensity bound
   double max_total_intensity_bound;
 
- public:
+ protected:
   /// @brief If set then it thresholds negative intensities
-  bool flag_threshold_negative_intensity;
+  bool threshold_negative_intensity = false;
 
 // Fields to deal with intensity track recording (itr)
  private :
@@ -233,6 +233,14 @@ class PP {
   /// @brief Gets Maximimum Total intensity bound that wwas encountered during realization
   inline double get_max_total_intensity_bound() { return max_total_intensity_bound; }
 
+  bool get_threshold_negative_intensity() const {
+    return threshold_negative_intensity;
+  }
+
+  void set_threshold_negative_intensity(const bool threshold_negative_intensity) {
+    this->threshold_negative_intensity = threshold_negative_intensity;
+  }
+
 ////////////////////////////////////////////////////////////////////////////////
 //                            Serialization
 ////////////////////////////////////////////////////////////////////////////////
@@ -248,7 +256,7 @@ class PP {
     ar(CEREAL_NVP(intensity));
     ar(CEREAL_NVP(flag_negative_intensity));
     ar(CEREAL_NVP(max_total_intensity_bound));
-    ar(CEREAL_NVP(flag_threshold_negative_intensity));
+    ar(CEREAL_NVP(threshold_negative_intensity));
     ar(CEREAL_NVP(itr_time));
     ar(CEREAL_NVP(itr_time_step));
     ar(CEREAL_NVP(itr));
@@ -271,7 +279,7 @@ class PP {
     ar(CEREAL_NVP(intensity));
     ar(CEREAL_NVP(flag_negative_intensity));
     ar(CEREAL_NVP(max_total_intensity_bound));
-    ar(CEREAL_NVP(flag_threshold_negative_intensity));
+    ar(CEREAL_NVP(threshold_negative_intensity));
     ar(CEREAL_NVP(itr_time));
     ar(CEREAL_NVP(itr_time_step));
     ar(CEREAL_NVP(itr));

tick/simulation/swig/simulation_module.i

@@ -41,6 +41,9 @@ class PP {
   void reseed_random_generator(int seed);
 
   SArrayDoublePtrList1D get_timestamps();
+
+  bool get_threshold_negative_intensity() const;
+  void set_threshold_negative_intensity(const bool threshold_negative_intensity);
 };
 
 

tick/simulation/tests/hawkes_test.py

@@ -178,6 +178,43 @@ def test_simu_hawkes_force_simulation(self):
             hawkes.force_simulation = True
             hawkes.simulate()
 
+    def test_hawkes_negative_intensity_fail(self):
+        """...Test simulation with negative kernel without threshold_negative_intensity
+        """
+        run_time = 40
+
+        hawkes = SimuHawkes(n_nodes=1, end_time=run_time, verbose=False,
+                            seed=1398)
+        kernel = HawkesKernelExp(-1.3, .8)
+        hawkes.set_kernel(0, 0, kernel)
+        hawkes.set_baseline(0, 0.3)
+
+        msg = 'Simulation stopped because intensity went negative ' \
+              '\(you could call ``threshold_negative_intensity`` to allow it\)'
+        with self.assertRaisesRegex(RuntimeError, msg):
+            hawkes.simulate()
+
+    def test_hawkes_negative_intensity(self):
+        """...Test simulation with negative kernel
+        """
+        run_time = 40
+
+        hawkes = SimuHawkes(n_nodes=1, end_time=run_time, verbose=False,
+                            seed=1398)
+        kernel = HawkesKernelExp(-1.3, .8)
+        hawkes.set_kernel(0, 0, kernel)
+        hawkes.set_baseline(0, 0.3)
+        hawkes.threshold_negative_intensity()
+
+        dt = 0.1
+        hawkes.track_intensity(dt)
+        hawkes.simulate()
+
+        self.assertAlmostEqual(hawkes.tracked_intensity[0].min(), 0)
+        self.assertAlmostEqual(hawkes.tracked_intensity[0].max(),
+                               hawkes.baseline[0])
+        self.assertGreater(hawkes.n_total_jumps, 1)
+
 
 if __name__ == "__main__":
     unittest.main()

tick/simulation/tests/src/hawkes_kernel_sumexp_gtest.cpp

@@ -70,17 +70,6 @@ TEST_F(HawkesKernelSumExpTest, is_zero) {
   EXPECT_FALSE(hawkes_kernel_sum_exp->is_zero());
 }
 
-TEST_F(HawkesKernelSumExpTest, constructor_negative_intensities_error) {
-  intensities[1] = -1;
-  auto kernel_1 = HawkesKernelSumExp(intensities, decays);
-  // this should work as we don't ask for the bound
-  kernel_1.get_convolution(1., timestamps, nullptr);
-  // this should throw an error as we ask to compute the future bound
-  auto kernel_2 = HawkesKernelSumExp(intensities, decays);
-  double bound;
-  ASSERT_THROW(kernel_2.get_convolution(1., timestamps, &bound), std::runtime_error);
-}
-
 TEST_F(HawkesKernelSumExpTest, constructor_negative_decays_error) {
   decays[1] = -1;
   ASSERT_THROW(HawkesKernelSumExp(intensities, decays), std::invalid_argument);