Merge pull request #84 from CHIMEFRB/ss/bugfix-2024-08

ss/bugfix-2024-08
CHIMEFRB · Sep 26, 2024 · d85cafa · d85cafa
2 parents 57a86a8 + dfbcf3b
commit d85cafa
Show file tree

Hide file tree

Showing 7 changed files with 587 additions and 362 deletions.
diff --git a/fitburst/analysis/fitter.py b/fitburst/analysis/fitter.py
@@ -484,6 +484,10 @@ def _compute_fit_statistics(self, spectrum_observed: float, fit_result: object)
             covariance = np.linalg.inv(0.5 * hessian)
             uncertainties = [float(x) for x in np.sqrt(np.diag(covariance)).tolist()]
 
+            if np.any(np.isnan(uncertainties)):
+                print("WARNING: one or more NaNs in exact uncertainties, replacing with approximates...")
+                uncertainties = [float(x) for x in np.sqrt(np.diag(covariance_approx)).tolist()]
+
             self.covariance_approx = covariance_approx
             self.covariance = covariance
             self.covariance_labels = par_labels

diff --git a/fitburst/analysis/model.py b/fitburst/analysis/model.py
@@ -342,18 +342,15 @@ def compute_profile(self, times: float, arrival_time: float, sc_time_ref: float,
         # compute either Gaussian or pulse-broadening function, depending on inputs.
         profile = np.zeros(times_copy.shape, dtype=float)
         sc_time = sc_time_ref * (freqs / ref_freq) ** sc_index
+        normalize = general["flags"]["normalize_pbf"]
 
-        if np.any(sc_time < np.fabs(0.15 * width)):
-            profile = rt.profile.compute_profile_gaussian(times_copy, arrival_time, width)
-
-        else:
-            # the following times array manipulates the times array so that we avoid a
-            # floating-point overlow in the exp((-times - toa) / sc_time) term in the
-            # PBF call. TODO: use a better, more transparent method for avoiding this.
+        if np.any(sc_time > 0.0):
             times_copy[times_copy < -5 * width] = -5 * width
             profile = rt.profile.compute_profile_pbf(
-                times_copy, arrival_time, width, freqs, ref_freq, sc_time_ref, sc_index=sc_index
+                times_copy, arrival_time, width, freqs, ref_freq, sc_time_ref, sc_index=sc_index, normalize=normalize
             )
+        else:
+            profile = rt.profile.compute_profile_gaussian(times_copy, arrival_time, width)
 
         # if data are folded and time/profile data contain two realizations, then
         # average along the appropriate axis to obtain a single realization.

diff --git a/fitburst/backend/general.yaml b/fitburst/backend/general.yaml
@@ -7,3 +7,6 @@ constants:
     dispersion: 4149.377593360996   # in units of MHz**2 * s * cm**3 / pc
     index_dispersion: -2.0
     index_scattering: -4.0
+
+flags:
+    normalize_pbf: False
diff --git a/fitburst/pipelines/fitburst_pipeline.py b/fitburst/pipelines/fitburst_pipeline.py
@@ -540,6 +540,9 @@
                     current_uncertainties = bestfit_uncertainties[current_parameter_label]
                     print(f"    * {current_parameter_label}: {current_list} +/- {current_uncertainties}")        
 
+                print("INFO: ratio of hessian matrix (approximate / exact):")
+                print(fitter.hessian / fitter.hessian_approx)
+
             # now create plots.
             filename_elems = input_file.split(".")
             output_string = ".".join(filename_elems[:len(filename_elems)-1])
@@ -561,6 +564,9 @@
                         "initial_time": data.times_bin0,
                         "model_parameters": current_params,
                         "fit_statistics": fitter.fit_statistics,
+                        "fit_logistics" : {
+                            "weight_range" : weight_range,
+                        }
                     },
                     out,
                     indent=4