From c8ac1fa50ecbdc394a4b995a347b0e463e6cbb87 Mon Sep 17 00:00:00 2001
From: Eric Teunis de Boone <ericteunis@deboone.nl>
Date: Wed, 17 May 2023 20:23:03 +0200
Subject: [PATCH] Pulse: show residuals in  SNR vs accuracy plot: WIP

---
 simulations/11_pulsed_timing.py | 32 ++++++++++++++++++++------------
 1 file changed, 20 insertions(+), 12 deletions(-)

diff --git a/simulations/11_pulsed_timing.py b/simulations/11_pulsed_timing.py
index 7de08de..f762011 100755
--- a/simulations/11_pulsed_timing.py
+++ b/simulations/11_pulsed_timing.py
@@ -506,9 +506,11 @@ if __name__ == "__main__":
     #
     # Find time accuracies as a function of signal strength
     #
-    time_accuracies = np.zeros((len(template_dts), len(snr_factors)))
-    mask_counts = np.zeros_like(time_accuracies)
-    for l, template_dt in tqdm(enumerate(template_dts)):
+    time_residuals_data = []
+
+    for a, template_dt in tqdm(enumerate(template_dts)):
+
+        time_residuals_data.append(np.zeros( (len(snr_factors), 3, N_residuals)))# res, snr, masked
 
         # Create the template
         # This is sampled at a lower samplerate than the interpolation template
@@ -526,26 +528,29 @@ if __name__ == "__main__":
             print()# separating tqdm
             print()# separating tqdm
 
-            # Make a plot of the time residuals
-            if N_residuals > 1:
+            wrong_peak_condition = lambda t_res: abs(t_res) > antenna_dt*4
+            mask = wrong_peak_condition(time_residuals)
+
+            # Save directly to large data array
+            time_residuals_data[a][k] = time_residuals, snrs, ~mask
+
+            # Make a plot of the time residuals <<<
+            if True and N_residuals > 1:
                 for i in range(1 + cut_wrong_peak_matches):
                     mask_count = 0
 
                     if i==1: # if cut_wrong_peak_matches:
-                        wrong_peak_condition = lambda t_res: abs(t_res) > antenna_dt*4
-
-                        mask = wrong_peak_condition(time_residuals)
-
                         mask_count = np.count_nonzero(mask)
 
-                        print("Masking {} residuals".format(mask_count))
                         time_residuals = time_residuals[~mask]
 
+                        # None masked
                         if not mask_count:
                             continue
 
-                    time_accuracies[l, k] = np.std(time_residuals)
-                    mask_counts[l, k] = mask_count
+                        # All masked
+                        if not len(time_residuals):
+                            continue
 
                     hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step')
                     fig, ax = plt.subplots()
@@ -614,8 +619,11 @@ if __name__ == "__main__":
 
                     if True:
                         plt.close(fig)
+                # >>> End of plot
 
+    #
     # SNR time accuracy plot
+    #
     if True:
         threshold_markers = ['^', 'v', '8', 'o']
         mask_thresholds = [np.inf, N_residuals*0.5, N_residuals*0.1, 1, 0]