From 7b19cc5ef09b3dbfb1573f26cb9da09111e5db73 Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Wed, 3 Jul 2024 14:16:24 +0200
Subject: [PATCH] correct empricial_cdf for case with weights and include tests

---
 seismostats/analysis/estimate_mc.py            | 15 +++++++--------
 seismostats/analysis/tests/test_estimate_mc.py | 15 +++++++++++++++
 2 files changed, 22 insertions(+), 8 deletions(-)

diff --git a/seismostats/analysis/estimate_mc.py b/seismostats/analysis/estimate_mc.py
index 01d562b..6bc7105 100644
--- a/seismostats/analysis/estimate_mc.py
+++ b/seismostats/analysis/estimate_mc.py
@@ -95,8 +95,8 @@ def empirical_cdf(
     if mc is None:
         mc = np.min(sample)
 
-    idx = np.argsort(sample)
-    x = sample[idx]
+    idx1 = np.argsort(sample)
+    x = sample[idx1]
     x, y_count = np.unique(x, return_counts=True)
 
     # add empty bins
@@ -106,15 +106,15 @@ def empirical_cdf(
         if mag_bin not in x:
             x = np.append(x, mag_bin)
             y_count = np.append(y_count, 0)
-    idx = np.argsort(x)
-    x = x[idx]
-    y_count = y_count[idx]
+    idx2 = np.argsort(x)
+    x = x[idx2]
+    y_count = y_count[idx2]
 
     # estimate the CDF
     if weights is None:
         y = np.cumsum(y_count) / len(sample)
     else:
-        weights_sorted = weights[idx]
+        weights_sorted = weights[idx1]
         y = np.cumsum(weights_sorted) / weights_sorted.sum()
 
         # make sure that y is zero if there are no samples in the first bins
@@ -122,8 +122,7 @@ def empirical_cdf(
             if y_loop > 0:
                 break
         leading_zeros = np.zeros(ii)
-
-        y = leading_zeros.append(y[np.cumsum(y_count[ii:]) - 1])
+        y = np.append(leading_zeros, y[np.cumsum(y_count[ii:]) - 1])
 
     return x, y
 
diff --git a/seismostats/analysis/tests/test_estimate_mc.py b/seismostats/analysis/tests/test_estimate_mc.py
index 2eb66d1..abd018c 100644
--- a/seismostats/analysis/tests/test_estimate_mc.py
+++ b/seismostats/analysis/tests/test_estimate_mc.py
@@ -128,6 +128,21 @@ def test_empirical_cdf(setup_magnitudes, delta_m=0.1):
     assert_equal(len(x), len(y))
     assert_equal(y[0], 0.06)
 
+    # test that weights function the way that they should
+    # 1. with equal weights
+    magnitudes = np.array([0.5, 0.6, 0.7, 0.8, 0.9])
+    weights = np.array([1, 1, 1, 1, 1])
+    x, y = empirical_cdf(magnitudes, mc=0, delta_m=0.1, weights=weights)
+    assert_almost_equal(x, np.arange(0, 1, 0.1))
+    assert_almost_equal(y, [0, 0, 0, 0, 0, 0.2, 0.4, 0.6, 0.8, 1])
+
+    # 2. with different weights
+    magnitudes = np.array([0.5, 0.6, 0.7, 0.8, 0.9])
+    weights = np.array([0.5, 0.5, 0.5, 0.5, 3])
+    x, y = empirical_cdf(magnitudes, mc=0, delta_m=0.1, weights=weights)
+    assert_almost_equal(x, np.arange(0, 1, 0.1))
+    assert_almost_equal(y, [0, 0, 0, 0, 0, 0.1, 0.2, 0.3, 0.4, 1])
+
 
 @pytest.fixture
 def setup_ks_dists():