In the presence of upper limits, correct the scaling factor of the models to...

In the presence of upper limits, correct the scaling factor of the models to the observations before computing the χ², not after.

CHANGELOG.md

@@ -8,6 +8,7 @@
 ### Fixed
 - To estimate parameters in log, pcigale determines which variables end with the "_log" string and removed it to compute the models. However in some circumstances, it was overzealous. This has been fixed. (Médéric Boquien)
 - When estimating a parameter in log, these were not scaled appropriately and taken in log when saving the related χ² and PDF. (Médéric Boquien)
+- In the presence of upper limits, correct the scaling factor of the models to the observations before computing the χ², not after. (Médéric Boquien)
 
 ### Optimised
 - Prior to version 0.7.0, we needed to maintain the list of redshifts for all the computed models. Past 0.7.0 we just infer the redshift from a list unique redshifts. This means that we can now discard the list of redshifts for all the models and only keep the list of unique redshifts. This saves ~8 MB of memory for every 10⁶ models. the models should be computed slightly faster but it is in the measurement noise. (Médéric Boquien)

pcigale/analysis_modules/pdf_analysis/utils.py

@@ -399,7 +399,17 @@ def compute_chi2(model_fluxes, obs_fluxes, obs_errors, lim_flag):
     scaling: array
         scaling of the models to obtain the minimum χ²
     """
+    limits = lim_flag and np.any(obs_errors <= 0.)
+
     scaling = _compute_scaling(model_fluxes, obs_fluxes, obs_errors)
+    # Some observations may not have flux values in some filter(s), but
+    # they can have upper limit(s).
+    if limits == True:
+        scaling_orig = scaling.copy()
+        for imod in range(len(model_fluxes)):
+            scaling[imod] = optimize.root(dchi2_over_ds2, scaling[imod],
+                                          args=(obs_fluxes, obs_errors,
+                                                model_fluxes[imod, :])).x
 
     # χ² of the comparison of each model to each observation.
     chi2 = np.zeros(model_fluxes.shape[0])
@@ -410,11 +420,7 @@ 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 (lim_flag and np.any(obs_errors <= 0.)) == True:
-        for imod in range(len(model_fluxes)):
-            scaling[imod] = optimize.root(dchi2_over_ds2, scaling[imod],
-                                          args=(obs_fluxes, obs_errors,
-                                                model_fluxes[imod, :])).x
+    if limits == True:
         mask_lim = (obs_errors <= 0.)
         chi2 += -2. * np.sum(
             np.log(