Fix a crash in the computation of the χ² because previous optimisations had...

Fix a crash in the computation of the χ² because previous optimisations had been been not fully propagated to the case of upper limits.

CHANGELOG.md

@@ -22,6 +22,7 @@
 - Remove a Numpy warning in the computation of the IGM absorption (Médéric Boquien)
 - The Herschel passbands only included the filter. Now they include the full throughput of the instrument. Flux differences are expected to be no more than 1%. (Médéric Boquien)
 - The `dustatt_powerlaw` module indicated that the Milky Way bump has an amplitude of 3. This is only valid for the `dustatt_calzleit` module. As `dustatt_powerlaw` is normalised to A(V) rather than E(B-V) for `dustatt_calzleit`, the bump is a factor Rv larger. A more reasonable value is now given. (Médéric Boquien)
+- The correction of the χ² for the upper limits now properly takes into account earlier changes intended to reduce memory usage and speed up the analysis. This prevents a crash. (Médéric Boquien)
 
 ### Optimised
 - By default the MKL library created many threads for each for the parallel processes. Not only was this not necessary as a high-level parallelisation already exists, but it generated a strong oversubscription on the CPU and on the RAM. The slowdown was over a factor of ~2 in some cases. Now we mandate KML to use only 1 thread fo each process. (Médéric Boquien & Yannick Roehlly)

pcigale/analysis_modules/pdf_analysis/utils.py

@@ -402,14 +402,11 @@ def compute_chi2(model_fluxes, obs_fluxes, obs_errors, lim_flag):
     # Some observations may not have flux values in some filter(s), but
     # they can have upper limit(s).
     if limits == True:
-        mask_lim = (obs_errors <= 0.)
-        chi2 += -2. * np.sum(
-            np.log(
-                np.sqrt(np.pi/2.)*(-obs_errors[mask_lim])*(
-                    1.+erf(
-                        (obs_fluxes[mask_lim]-model_fluxes[mask_lim, :] *
-                         scaling[:, np.newaxis]) /
-                        (np.sqrt(2)*(-obs_errors[mask_lim]))))), axis=1)
+        for i, obs_error in enumerate(obs_errors):
+            if obs_error < 0.:
+                chi2 += -2. * np.log(np.sqrt(np.pi/2.) * (-obs_errors[i]) * (
+                        1.+erf((obs_fluxes[i] - model_fluxes[i, :]*scaling) /
+                            (np.sqrt(2)*(-obs_errors[i])))))
 
     return chi2, scaling