Repository: SCECcode/pycsep (language: python)

Issue:
T-test does not work for single magnitude bin
Hi,
I noticed that T-test does not work for single magnitude bins. E.g. mbins = [5.95].

The error is traced down to "target_event_rates --> get_rates --> get_magnitude_index --> bin1d_vec" to line 75: h = bins[1] - bins[0].

Since there is only one bin, so bins[1] gives error. For now, I manually put h=0.1 in my copy of repo and its working.

Perhaps, as a solution, a check can be put here, that if there size(bin)==1, then default h=0.1. It is already taking last bins as open. So it should work. 

Cheers,
Thanks


T-test does not work for single magnitude bin
Hi,
I noticed that T-test does not work for single magnitude bins. E.g. mbins = [5.95].

The error is traced down to "target_event_rates --> get_rates --> get_magnitude_index --> bin1d_vec" to line 75: h = bins[1] - bins[0].

Since there is only one bin, so bins[1] gives error. For now, I manually put h=0.1 in my copy of repo and its working.

Perhaps, as a solution, a check can be put here, that if there size(bin)==1, then default h=0.1. It is already taking last bins as open. So it should work. 

Cheers,
Thanks




Discussion / hints:
Hi @khawajasim! 

Good find on this issue. I'll work on a fix for this tomorrow. Given the implementation right noow, we'd expect that bin to be [5.95, inf) instead of [5.95, 6.05]. Do you want to be able to specify a finite bin-width if there is only a single bin?

In other words, It seems if we are treating that bin as open, we'd want to just check for that case and return index=0 if tthe magnitude >= 5.96. Do you agree?


Hi Bill,

Thank you for your response. Since I am just playing around with actual or dummy data, so I put all the magnitudes in just one bin, say 5.95 etc. And I had only one bin for all the data, so I just putting dh = 0.1 as a quick fix, which worked for the time being. But it may not be a robust solution for someone else with single bin. So yes better solution is [5.95, inf). 


Sounds good @khawajasim. I'm working on a fix for this right now. I'll ask you to review the code before I merge into master.
Sure.
Hi @khawajasim! 

Good find on this issue. I'll work on a fix for this tomorrow. Given the implementation right noow, we'd expect that bin to be [5.95, inf) instead of [5.95, 6.05]. Do you want to be able to specify a finite bin-width if there is only a single bin?

In other words, It seems if we are treating that bin as open, we'd want to just check for that case and return index=0 if tthe magnitude >= 5.96. Do you agree?


Hi Bill,

Thank you for your response. Since I am just playing around with actual or dummy data, so I put all the magnitudes in just one bin, say 5.95 etc. And I had only one bin for all the data, so I just putting dh = 0.1 as a quick fix, which worked for the time being. But it may not be a robust solution for someone else with single bin. So yes better solution is [5.95, inf). 


Sounds good @khawajasim. I'm working on a fix for this right now. I'll ask you to review the code before I merge into master.
Sure.

Fix patch:
diff --git a/csep/core/poisson_evaluations.py b/csep/core/poisson_evaluations.py
index 7f3a6564..0901d1f9 100644
--- a/csep/core/poisson_evaluations.py
+++ b/csep/core/poisson_evaluations.py
@@ -3,6 +3,7 @@
 
 import numpy
 import scipy.stats
+import scipy.spatial
 
 from csep.models import EvaluationResult
 from csep.utils.stats import poisson_joint_log_likelihood_ndarray
@@ -226,7 +227,7 @@ def poisson_spatial_likelihood(forecast, catalog):
     scale = catalog.event_count / forecast.event_count
     
     first_term = -forecast.spatial_counts() * scale
-    second_term = catalog.spatial_counts() * np.log(forecast.spatial_counts() * scale)
+    second_term = catalog.spatial_counts() * numpy.log(forecast.spatial_counts() * scale)
     third_term = -scipy.special.loggamma(catalog.spatial_counts() + 1)
     
     poll = first_term + second_term + third_term
diff --git a/csep/utils/calc.py b/csep/utils/calc.py
index 3ac07461..e0071f47 100644
--- a/csep/utils/calc.py
+++ b/csep/utils/calc.py
@@ -72,7 +72,12 @@ def bin1d_vec(p, bins, tol=None, right_continuous=False):
     bins = numpy.array(bins)
     p = numpy.array(p)
     a0 = numpy.min(bins)
-    h = bins[1] - bins[0]
+    # if user supplies only a single bin, do 2 things: 1) fix right continuous to true, and use of h is arbitrary
+    if bins.size == 1:
+        right_continuous = True
+        h = 1
+    else:
+        h = bins[1] - bins[0]
 
     a0_tol = numpy.abs(a0) * numpy.finfo(numpy.float).eps
     h_tol = numpy.abs(h) * numpy.finfo(numpy.float).eps


Test patch:
diff --git a/tests/test_calc.py b/tests/test_calc.py
index 0dea9703..618443d1 100644
--- a/tests/test_calc.py
+++ b/tests/test_calc.py
@@ -98,6 +98,14 @@ def test_bin1d_vec_int(self):
         expected = [0, 0, 0, 1, 2]
         self.assertListEqual(test.tolist(), expected)
 
+    def test_bin1d_single_bin1(self):
+        data = [-1, 0, 2, 3, 1, 1.5, 1.0, 0.999999999999999]
+        bin_edges = [1]
+        # purposely leaving right_continous flag=False bc it should be forced in the bin1d_vec function
+        test = bin1d_vec(data, bin_edges)
+        expected = [-1, -1, 0, 0, 0, 0, 0, -1]
+        self.assertListEqual(test.tolist(), expected)
+
     def test_upper_limit_right_continuous(self):
         data = [40, 40, 40]
         bin_edges = [0, 10, 20, 30]