CLEANUP: code in plotting.fw now modular

multi_locus_analysis/finite_window/__init__.py

@@ -64,3 +64,4 @@
 from .simulation import *
 from .munging import *
 from .stats import *
+from .stats import _ccdf_int

multi_locus_analysis/finite_window/simulation.py

@@ -5,7 +5,7 @@
 
 import bruno_util
 import bruno_util.random
-from .. import stats
+from .. import stats as _mla_stats
 
 
 @bruno_util.random.strong_default_seed
@@ -253,7 +253,7 @@ def _boot_int_i(N_var_T):
         ].copy()
         window_sizes = obs.groupby('replicate')['window_size'].first().values
         # now sorted
-        window_sizes, window_cdf = stats.ecdf(window_sizes)
+        window_sizes, window_cdf = _mla_stats.ecdf(window_sizes)
         window_sf = 1 - window_cdf
         all_times, cdf_int, cdf_ext, Z_X, F_T = fw.ecdf_ext_int(
             exterior.wait_time.values,

multi_locus_analysis/finite_window/stats.py

@@ -5,13 +5,20 @@
 
 from ..stats import ecdf
 
+
+def window_sf(obs, traj_cols):
+    window_sizes = obs.groupby(traj_cols)['window_size'].first().values
+    # now sorted
+    window_sizes, window_cdf = ecdf(window_sizes, pad_left_at_x=0)
+    return 1 - window_cdf
+
+
 def ecdf_combined(exterior, interior, window_sizes, ext_bins='auto', **kwargs):
 
     all_times, cdf_int, cdf_ext, Z_X, F_T = ecdf_ext_int(
         exterior, interior, window_sizes, **kwargs
     )
 
-    N_ext = len(exterior)
     y, t_bins = np.histogram(exterior, bins=ext_bins)
     bin_centers = (t_bins[:-1] + t_bins[1:]) / 2
     Z_hist = np.sum(y*np.diff(t_bins))
@@ -37,6 +44,15 @@ def ecdf_combined(exterior, interior, window_sizes, ext_bins='auto', **kwargs):
     )
     return bin_centers, final_est
 
+
+def _ccdf_int(x_int, cdf_int):
+    if not np.isclose(cdf_int[0], 0):
+        raise ValueError("Need to integrate interior CDF, but values start >0")
+    ccdf_int = np.zeros_like(cdf_int)
+    ccdf_int[1:] = np.cumsum(cdf_int[1:] * np.diff(x_int))
+    return ccdf_int
+
+
 def ecdf_ext_int(exterior, interior, window_sizes, window_sf=None,
                  times_allowed=None, pad_left_at_x=0):
     if len(exterior) < 5 or len(interior) < 5:
@@ -48,9 +64,8 @@ def ecdf_ext_int(exterior, interior, window_sizes, window_sf=None,
         interior, window_sizes, window_sf, times_allowed,
         pad_left_at_x=pad_left_at_x,
     )
-    # now compute integral of CDF w.r.t. t
-    ccdf_int = np.zeros_like(cdf_int)
-    ccdf_int[1:] = np.cumsum(cdf_int[1:] * np.diff(x_int))
+    # integral of CDF w.r.t. t
+    ccdf_int = _ccdf_int(x_int, cdf_int)
 
     if times_allowed is None:
         times_allowed = np.sort(np.concatenate((x_int, x_ext)))
@@ -144,15 +159,6 @@ def ecdf_windowed(
     else:
         x = np.unique(y)
     x.sort()
-    if auto_pad_left:
-        dx = np.mean(np.diff(x))
-        x = np.insert(x, 0, x[0] - dx)
-    elif pad_left_at_x is not None:
-        if x[0] <= pad_left_at_x:
-            warnings.warn('pad_left_at_x not left of x in ecdf_windowed! '
-                          'Ignoring...')
-        else:
-            x = np.insert(x, 0, pad_left_at_x)
 
     num_obs = len(y)
     cdf = np.zeros(x.shape, dtype=np.dtype('float'))
@@ -179,6 +185,17 @@ def ecdf_windowed(
         cdf[xi] = full_cdf[i]
     if normalize:
         cdf = cdf/full_cdf[-1]
+    if auto_pad_left or pad_left_at_x is not None:
+        cdf = np.insert(cdf, 0, 0)
+    if auto_pad_left:
+        dx = np.mean(np.diff(x))
+        x = np.insert(x, 0, x[0] - dx)
+    elif pad_left_at_x is not None:
+        if x[0] <= pad_left_at_x:
+            warnings.warn('pad_left_at_x not left of x in ecdf_windowed! '
+                          'Ignoring...')
+        else:
+            x = np.insert(x, 0, pad_left_at_x)
     return x, cdf