Commit 7a90054b authored by Médéric Boquien's avatar Médéric Boquien

Store the upper limits and the regular fluxes/properties in separate...

Store the upper limits and the regular fluxes/properties in separate dictionaries so that we do not have to test when fitting whether they are upper limits of not.
parent ad67241e
......@@ -222,33 +222,30 @@ def compute_chi2(models, obs, corr_dz, wz, lim_flag):
# Multiplications are faster than divisions, so we directly use the
# inverse error
inv_flux_err = 1. / obs.flux_err[band]
if np.isfinite(flux) and inv_flux_err > 0.:
model = models.flux[band][wz]
chi2 += ((flux - model * scaling) * inv_flux_err) ** 2.
model = models.flux[band][wz]
chi2 += ((flux - model * scaling) * inv_flux_err) ** 2.
# Computation of the χ² from intensive properties
for name, prop in obs.intprop.items():
if np.isfinite(prop):
model = models.intprop[name][wz]
chi2 += ((prop - model) * (1. / obs.intprops_err[name])) ** 2.
model = models.intprop[name][wz]
chi2 += ((prop - model) * (1. / obs.intprops_err[name])) ** 2.
# Computation of the χ² from extensive properties
for name, prop in obs.extprop.items():
# Multiplications are faster than divisions, so we directly use the
# inverse error
inv_prop_err = 1. / obs.extprop_err[name]
if np.isfinite(prop) and inv_prop_err > 0.:
model = models.extprop[name][wz]
chi2 += ((prop - (scaling * model) * corr_dz) * inv_prop_err) ** 2.
model = models.extprop[name][wz]
chi2 += ((prop - (scaling * model) * corr_dz) * inv_prop_err) ** 2.
# Finally take the presence of upper limits into account
if limits == True:
for band, obs_error in obs.flux_err.items():
if obs_error < 0.:
chi2 -= 2. * np.log(.5 *
(1. + erf(((obs.fluxes[band] -
model_flux[band][wz] * scaling) /
(-np.sqrt(2.)*obs_error)))))
for band, obs_error in obs.flux_ul_err.items():
model = models.flux[band][wz]
chi2 -= 2. * np.log(.5 *
(1. + erf(((obs.flux_ul[band] -
model * scaling) /
(-np.sqrt(2.)*obs_error)))))
return chi2, scaling
......
......@@ -306,9 +306,20 @@ class Observation(object):
self.distance = cosmo.luminosity_distance(self.redshift).value
else:
self.distance = np.nan
self.flux = {band: row[band] for band in cls.bands}
self.flux_err = {band: row[band + '_err'] for band in cls.bands}
self.intprop = {prop: row[prop] for prop in cls.intprops}
self.intprop_err = {prop: row[prop + '_err'] for prop in cls.intprops}
self.extprop = {prop: row[prop] for prop in cls.extprops}
self.extprop_err = {prop: row[prop + '_err'] for prop in cls.extprops}
self.flux = {k: row[k] for k in cls.bands
if np.isfinite(row[k]) and row[k + '_err'] > 0.}
self.flux_ul = {k: row[k] for k in cls.bands
if np.isfinite(row[k]) and row[k + '_err'] <= 0.}
self.flux_err = {k: row[k + '_err'] for k in self.flux.keys()}
self.flux_ul_err = {k: row[k + '_err'] for k in self.flux_ul.keys()}
self.extprop = {k: row[k] for k in cls.extprops
if np.isfinite(row[k]) and row[k + '_err'] > 0.}
self.extprop_ul = {k: row[k] for k in cls.extprops
if np.isfinite(row[k]) and row[k + '_err'] <= 0.}
self.extprop_err = {k: row[k + '_err'] for k in self.extprop.keys()}
self.extprop_ul_err = {k: row[k + '_err']
for prop in self.extprop_ul.keys()}
self.intprop = {k: row[k] for k in cls.intprops}
self.intprop_err = {k: row[k + '_err'] for k in cls.intprops}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment