From 885eb59f5a4d5e9478571b90c1d8b89065b2f15b Mon Sep 17 00:00:00 2001 From: Martin B Date: Mon, 23 Oct 2017 14:42:39 +0200 Subject: [PATCH] Allow simulation of point processes with thresholded negative intensities --- tick/simulation/base/simu_point_process.py | 11 ++++++ tick/simulation/src/hawkes.cpp | 5 ++- .../hawkes_kernels/hawkes_kernel_sum_exp.cpp | 9 +++-- tick/simulation/src/pp.cpp | 26 +++++++------ tick/simulation/src/pp.h | 16 ++++++-- tick/simulation/swig/simulation_module.i | 3 ++ tick/simulation/tests/hawkes_test.py | 37 +++++++++++++++++++ .../tests/src/hawkes_kernel_sumexp_gtest.cpp | 11 ------ 8 files changed, 87 insertions(+), 31 deletions(-) diff --git a/tick/simulation/base/simu_point_process.py b/tick/simulation/base/simu_point_process.py index e81217714..3fa58fd42 100644 --- a/tick/simulation/base/simu_point_process.py +++ b/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) diff --git a/tick/simulation/src/hawkes.cpp b/tick/simulation/src/hawkes.cpp index 54187ac95..6598c30fe 100644 --- a/tick/simulation/src/hawkes.cpp +++ b/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; diff --git a/tick/simulation/src/hawkes_kernels/hawkes_kernel_sum_exp.cpp b/tick/simulation/src/hawkes_kernels/hawkes_kernel_sum_exp.cpp index aa1592305..4c19abf8c 100644 --- a/tick/simulation/src/hawkes_kernels/hawkes_kernel_sum_exp.cpp +++ b/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; diff --git a/tick/simulation/src/pp.cpp b/tick/simulation/src/pp.cpp index 2d380ed20..de9c68615 100644 --- a/tick/simulation/src/pp.cpp +++ b/tick/simulation/src/pp.cpp @@ -116,22 +116,23 @@ 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; } @@ -139,8 +140,8 @@ void PP::simulate(double end_time, ulong n_points) { // 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,7 +165,7 @@ 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 @@ -172,8 +173,9 @@ void PP::simulate(double end_time, ulong n_points) { // 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 diff --git a/tick/simulation/src/pp.h b/tick/simulation/src/pp.h index 93bc01520..055787782 100644 --- a/tick/simulation/src/pp.h +++ b/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)); diff --git a/tick/simulation/swig/simulation_module.i b/tick/simulation/swig/simulation_module.i index 038c41ecd..90da63288 100644 --- a/tick/simulation/swig/simulation_module.i +++ b/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); }; diff --git a/tick/simulation/tests/hawkes_test.py b/tick/simulation/tests/hawkes_test.py index eff945632..c7c978728 100644 --- a/tick/simulation/tests/hawkes_test.py +++ b/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() diff --git a/tick/simulation/tests/src/hawkes_kernel_sumexp_gtest.cpp b/tick/simulation/tests/src/hawkes_kernel_sumexp_gtest.cpp index 15d7d4284..bee455e43 100644 --- a/tick/simulation/tests/src/hawkes_kernel_sumexp_gtest.cpp +++ b/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);