Merge pull request #63 from eurunuela/surrogate_histograms

Change how surrogate histograms are calculated

Author: eurunuela

connPFM/connectivity/connectivity_utils.py

@@ -2,7 +2,6 @@
 import logging
 
 import numpy as np
-from joblib import Parallel, delayed
 from scipy.stats import zscore
 
 LGR = logging.getLogger(__name__)
@@ -19,7 +18,7 @@ def calculate_ets(y, n):
     return ets, u, v
 
 
-def rss_surr(z_ts, u, v, surrprefix, sursufix, masker, irand):
+def rss_surr(z_ts, u, v, surrprefix, sursufix, masker, irand, nbins, hist_range=(0, 1)):
     """Calculate RSS on surrogate data."""
     [t, n] = z_ts.shape
 
@@ -42,7 +41,10 @@ def rss_surr(z_ts, u, v, surrprefix, sursufix, masker, irand):
     # calcuate rss
     rssr = np.sqrt(np.sum(np.square(etsr), axis=1))
 
-    return (rssr, etsr, np.min(etsr), np.max(etsr))
+    # Calculate histogram
+    ets_hist, bin_edges = np.histogram(etsr.flatten(), bins=nbins, range=hist_range)
+
+    return (rssr, etsr, ets_hist, bin_edges)
 
 
 def remove_neighboring_peaks(rss, idx):
@@ -136,35 +138,20 @@ def calculate_hist_threshold(hist, bins, percentile=95):
     return thr
 
 
-def surrogates_histogram(
-    surrprefix,
-    sursufix,
-    masker,
-    hist_range,
-    numrand=100,
-    nbins=500,
-    percentile=95,
+def sum_histograms(
+    hist_list,
 ):
     """
-    Read AUCs of surrogates, calculate histogram and sum of all histograms to
+    Get histograms of all surrogates and sum them to
     obtain a single histogram that summarizes the data.
     """
-    ets_hist = np.zeros((numrand, nbins))
-
-    # calculate histogram for each surrogate
-    hist = Parallel(n_jobs=-1, backend="multiprocessing")(
-        delayed(calculate_hist)(surrprefix, sursufix, irand, masker, hist_range, nbins)
-        for irand in range(numrand)
-    )
 
-    for irand in range(numrand):
-        ets_hist[irand, :] = hist[irand][0]
+    # Initialize matrix to store surrogate histograms
+    all_hists = np.zeros((len(hist_list), hist_list[0][3].shape[0] - 1))
 
-    bin_edges = hist[0][1]
+    for rand_idx in range(len(hist_list)):
+        all_hists[rand_idx, :] = hist_list[rand_idx][2]
 
-    ets_hist_sum = np.sum(ets_hist, axis=0)
+    ets_hist_sum = np.sum(all_hists, axis=0)
 
-    # calculate histogram threshold
-    thr = calculate_hist_threshold(ets_hist_sum, bin_edges, percentile)
-
-    return thr
+    return ets_hist_sum

connPFM/connectivity/ev.py

@@ -21,6 +21,7 @@ def event_detection(
     nsur=100,
     segments=True,
     peak_detection="rss",
+    nbins=1000,
 ):
     """Perform event detection on given data."""
     masker = NiftiLabelsMasker(
@@ -55,19 +56,14 @@ def event_detection(
     rssr = np.zeros([t, nsur])
 
     # calculate ets and rss of surrogate data
+    LGR.info("Calculating edge-time matrix, RSS and histograms for surrogates...")
     surrogate_events = Parallel(n_jobs=-1, backend="multiprocessing")(
-        delayed(connectivity_utils.rss_surr)(z_ts, u, v, surrprefix, sursufix, masker, irand)
+        delayed(connectivity_utils.rss_surr)(
+            z_ts, u, v, surrprefix, sursufix, masker, irand, nbins
+        )
         for irand in range(nsur)
     )
 
-    hist_ranges = np.zeros((2, nsur))
-    for irand in range(nsur):
-        hist_ranges[0, irand] = surrogate_events[irand][2]
-        hist_ranges[1, irand] = surrogate_events[irand][3]
-
-    hist_min = np.min(hist_ranges, axis=1)[0]
-    hist_max = np.max(hist_ranges, axis=1)[1]
-
     # Make selection of points with RSS
     if "rss" in peak_detection:
         LGR.info("Selecting points with RSS...")
@@ -117,13 +113,13 @@ def event_detection(
         LGR.info("Selecting points with edge time-series matrix...")
         if peak_detection == "ets":
             LGR.info("Reading AUC of surrogates to perform the thresholding step...")
-            thr = connectivity_utils.surrogates_histogram(
-                surrprefix,
-                sursufix,
-                masker,
-                hist_range=(hist_min, hist_max),
-                numrand=nsur,
+            hist_sum = connectivity_utils.sum_histograms(
+                surrogate_events,
             )
+            thr = connectivity_utils.calculate_hist_threshold(
+                hist_sum, surrogate_events[0][3][:-1], percentile=95
+            )
+
         elif peak_detection == "ets_time":
             # Initialize array for threshold
             thr = np.zeros(t)
@@ -137,9 +133,7 @@ def event_detection(
                     sur_ets_at_time[sur_idx, :] = surrogate_events[sur_idx][1][time_idx, :]
 
                 # calculate histogram of all surrogate ets at time point
-                hist, bins = np.histogram(
-                    sur_ets_at_time.flatten(), bins=500, range=(hist_min, hist_max)
-                )
+                hist, bins = np.histogram(sur_ets_at_time.flatten(), bins=nbins, range=(0, 1))
 
                 # calculate threshold for time point
                 thr[time_idx] = connectivity_utils.calculate_hist_threshold(

connPFM/tests/test_ev.py

@@ -26,11 +26,12 @@ def test_rss_surr(AUC_file, atlas_file, surr_dir, rssr_auc_file):
         standardize=False,
         strategy="mean",
     )
+    nbins = 1000
 
     AUC_img = masker.fit_transform(AUC_file)
     _, u, v = connectivity_utils.calculate_ets(AUC_img, AUC_img.shape[1])
     rssr, _, _, _ = connectivity_utils.rss_surr(
-        AUC_img, u, v, join(surr_dir, "surrogate_AUC_"), "", masker, 0
+        AUC_img, u, v, join(surr_dir, "surrogate_AUC_"), "", masker, 0, nbins
     )
     rssr_auc = np.loadtxt(rssr_auc_file)
     assert np.all(np.isclose(rssr, rssr_auc))