reformatted according to black

Author: Nabil Brice

stingray/pulse/pulsar.py

@@ -209,14 +209,17 @@ def phase_exposure(start_time, stop_time, period, nbin=16, gti=None):
     return expo / np.max(expo)
 
 
-def fold_events(times, *frequency_derivatives,
-                mode = "ef",
-                nbin = 16,
-                weights = 1,
-                gti = None,
-                ref_time = 0,
-                expocorr = False,
-                variances = None):
+def fold_events(
+    times,
+    *frequency_derivatives,
+    mode="ef",
+    nbin=16,
+    weights=1,
+    gti=None,
+    ref_time=0,
+    expocorr=False,
+    variances=None,
+):
     """Epoch folding with exposure correction.
 
     By default, the keyword `times` accepts a list of
@@ -321,7 +324,9 @@ def fold_events(times, *frequency_derivatives,
 
         # For weighted pdm, the measured uncertainties (variances) for each event are required
         # set the event_weights as 1/variances or else 1s if variances=None
-        event_weights = (1.0/np.asarray(variances)) if variances is not None else np.ones_like(weights)
+        event_weights = (
+            (1.0 / np.asarray(variances)) if variances is not None else np.ones_like(weights)
+        )
 
         bins = np.linspace(0, 1, nbin + 1)
         # phases are already fractional values in [0, 1]
@@ -331,19 +336,16 @@ def fold_events(times, *frequency_derivatives,
         bin_indices = np.clip(bin_indices, 0, nbin - 1)
 
         # for efficiency, the sum-of-squared deviations for each bin is split
-        # ss_dev = sum_in_bin( (values_in_bin - mean_in_bin)**2 )
+        # sosdev = sum_in_bin( (values_in_bin - mean_in_bin)**2 )
         # mean_in_bin = sum_in_bin(values_in_bin) / n_in_bin
-        # so need n_in_bin, sum_in_bin, sumofsquares_in_bin
+        # so need n_in_bin, sum_in_bin, sos_in_bin (sum-of-squares)
         # now compute the stats for each bin using bincount
         # minlength=nbin prevents array shortening when a bin index is not represented in bin_indices
-        n_in_bin = np.bincount(bin_indices, weights=event_weights,
-                                minlength=nbin)
-        sum_in_bin = np.bincount(bin_indices, weights=weights * event_weights,
-                                minlength=nbin)
-        sumofsquares_in_bin = np.bincount(bin_indices, weights=weights**2 * event_weights,
-                                minlength=nbin)
-        # put it together; avoid division by zero
-        raw_profile = sumofsquares_in_bin - sum_in_bin**2 / np.maximum(n_in_bin, 1)
+        n_in_bin = np.bincount(bin_indices, weights=event_weights, minlength=nbin)
+        sum_in_bin = np.bincount(bin_indices, weights=weights * event_weights, minlength=nbin)
+        sos_in_bin = np.bincount(bin_indices, weights=weights**2 * event_weights, minlength=nbin)
+        # put it together to make sum-of-squared deviations: sosdev; avoid division by zero
+        raw_profile = sos_in_bin - sum_in_bin**2 / np.maximum(n_in_bin, 1)
 
         # dummy array for the error, which we don't have for the variance
         raw_profile_err = np.zeros_like(raw_profile)

stingray/pulse/search.py

@@ -172,8 +172,15 @@ def stat_fun(t, f, fd=0, **kwargs):
 
 
 def phase_dispersion_search(
-    times, flux, frequencies, nbin=128, segment_size=5000, expocorr=False, gti=None, fdots=0,
-    variances=None
+    times,
+    flux,
+    frequencies,
+    nbin=128,
+    segment_size=5000,
+    expocorr=False,
+    gti=None,
+    fdots=0,
+    variances=None,
 ):
     """Performs folding at trial frequencies in time series data (i.e.~a light curve
     of flux or photon counts) and computes the Phase Dispersion Minimization statistic.
@@ -232,13 +239,17 @@ def stat_fun(t, f, fd=0, **kwargs):
         len_flux = len(flux)
         # Get the local (possibly modified copy of variances) from kwargs rather than as a closure
         variances_local = kwargs.get("variances", None)
-        event_weights = (1.0/np.asarray(variances_local)) if variances_local is not None else np.ones_like(flux)
+        event_weights = (
+            (1.0 / np.asarray(variances_local))
+            if variances_local is not None
+            else np.ones_like(flux)
+        )
         mean_grand = np.average(flux, weights=event_weights)
-        dev2_grand = (flux - mean_grand)**2
-        
-        sum_dev2_grand = np.dot( dev2_grand, event_weights )
+        dev2_grand = (flux - mean_grand) ** 2
+
+        sosdev_grand = np.dot(dev2_grand, event_weights)
 
-        return pdm_profile_stat(profile, sum_dev2_grand, len_flux)
+        return pdm_profile_stat(profile, sosdev_grand, len_flux)
 
     return _folding_search(
         stat_fun,