Commit c620bb2d authored by  Hector Salas's avatar Hector Salas

Modify ObservationsManagerPassbands to iterate over an instance of...

Modify ObservationsManagerPassbands to iterate over an instance of observationclass,  modify Workers accordingly and to include the properties and modify pdf_anaysis/utils to include properties in compute_chi2, _cimpute_scaling and dchi2_over_ds2
parent 793aa346
...@@ -53,7 +53,7 @@ def compute_corr_dz(model_z, obs_z): ...@@ -53,7 +53,7 @@ def compute_corr_dz(model_z, obs_z):
return (cosmo.luminosity_distance(obs_z).value * 1e5)**2. return (cosmo.luminosity_distance(obs_z).value * 1e5)**2.
def dchi2_over_ds2(s, obs_fluxes, obs_errors, mod_fluxes): def dchi2_over_ds2(s, obs_values, obs_errors, mod_values):
"""Function used to estimate the normalization factor in the SED fitting """Function used to estimate the normalization factor in the SED fitting
process when upper limits are included in the dataset to fit (from Eq. A11 process when upper limits are included in the dataset to fit (from Eq. A11
in Sawicki M. 2012, PASA, 124, 1008). in Sawicki M. 2012, PASA, 124, 1008).
...@@ -63,12 +63,13 @@ def dchi2_over_ds2(s, obs_fluxes, obs_errors, mod_fluxes): ...@@ -63,12 +63,13 @@ def dchi2_over_ds2(s, obs_fluxes, obs_errors, mod_fluxes):
s: Float s: Float
Contains value onto which we perform minimization = normalization Contains value onto which we perform minimization = normalization
factor factor
obs_fluxes: RawArray obs_value: RawArray
Contains observed fluxes for each filter. Contains observed fluxes for each filter and obseved extensive
properties.
obs_errors: RawArray obs_errors: RawArray
Contains observed errors for each filter. Contains observed errors for each filter and extensive properties.
model_fluxes: RawArray model_values: RawArray
Contains modeled fluxes for each filter. Contains modeled fluxes for each filter and extensive properties.
lim_flag: Boolean lim_flag: Boolean
Tell whether we use upper limits (True) or not (False). Tell whether we use upper limits (True) or not (False).
...@@ -88,28 +89,28 @@ def dchi2_over_ds2(s, obs_fluxes, obs_errors, mod_fluxes): ...@@ -88,28 +89,28 @@ def dchi2_over_ds2(s, obs_fluxes, obs_errors, mod_fluxes):
wlim = np.where(np.isfinite(obs_errors) & (obs_errors < 0.)) wlim = np.where(np.isfinite(obs_errors) & (obs_errors < 0.))
wdata = np.where(obs_errors >= 0.) wdata = np.where(obs_errors >= 0.)
mod_fluxes_data = mod_fluxes[wdata] mod_value_data = mod_value[wdata]
mod_fluxes_lim = mod_fluxes[wlim] mod_value_lim = mod_value[wlim]
obs_fluxes_data = obs_fluxes[wdata] obs_value_data = obs_value[wdata]
obs_fluxes_lim = obs_fluxes[wlim] obs_value_lim = obs_value[wlim]
obs_errors_data = obs_errors[wdata] obs_errors_data = obs_errors[wdata]
obs_errors_lim = -obs_errors[wlim] obs_errors_lim = -obs_errors[wlim]
dchi2_over_ds_data = np.sum( dchi2_over_ds_data = np.sum(
(obs_fluxes_data-s*mod_fluxes_data) * (obs_value_data-s*mod_value_data) *
mod_fluxes_data/(obs_errors_data*obs_errors_data)) mod_value_data/(obs_errors_data*obs_errors_data))
dchi2_over_ds_lim = np.sqrt(2./np.pi)*np.sum( dchi2_over_ds_lim = np.sqrt(2./np.pi)*np.sum(
mod_fluxes_lim*np.exp( mod_value_lim*np.exp(
-np.square( -np.square(
(obs_fluxes_lim-s*mod_fluxes_lim)/(np.sqrt(2)*obs_errors_lim) (obs_value_lim-s*mod_value_lim)/(np.sqrt(2)*obs_errors_lim)
) )
)/( )/(
obs_errors_lim*( obs_errors_lim*(
1.+erf( 1.+erf(
(obs_fluxes_lim-s*mod_fluxes_lim)/(np.sqrt(2)*obs_errors_lim) (obs_fluxes_lim-s*mod_value_lim)/(np.sqrt(2)*obs_errors_lim)
) )
) )
) )
...@@ -120,7 +121,7 @@ def dchi2_over_ds2(s, obs_fluxes, obs_errors, mod_fluxes): ...@@ -120,7 +121,7 @@ def dchi2_over_ds2(s, obs_fluxes, obs_errors, mod_fluxes):
return func return func
def _compute_scaling(model_fluxes, obs_fluxes, obs_errors): def _compute_scaling(model_fluxes, model_propsmass, observation):
"""Compute the scaling factor to be applied to the model fluxes to best fit """Compute the scaling factor to be applied to the model fluxes to best fit
the observations. Note that we look over the bands to avoid the creation of the observations. Note that we look over the bands to avoid the creation of
an array of the same size as the model_fluxes array. Because we loop on the an array of the same size as the model_fluxes array. Because we loop on the
...@@ -130,28 +131,42 @@ def _compute_scaling(model_fluxes, obs_fluxes, obs_errors): ...@@ -130,28 +131,42 @@ def _compute_scaling(model_fluxes, obs_fluxes, obs_errors):
---------- ----------
model_fluxes: array model_fluxes: array
Fluxes of the models Fluxes of the models
obs_fluxes: array model_propsmass: array
Observed fluxes Extensive properties of the models to be fitted
obs_errors: array observation: Class
Observed errors Class instance containing the fluxes, intensive properties, extensive
properties and their errors, for a sigle observation.
Returns Returns
------- -------
scaling: array scaling: array
Scaling factors minimising the χ² Scaling factors minimising the χ²
""" """
obs_fluxes = observation.fluxes
obs_fluxes_err = observation.fluxes_err
obs_propsmass = observation.extprops
obs_propsmass_err = observation.extprops_err
num = np.zeros(model_fluxes.shape[1]) num = np.zeros(model_fluxes.shape[1])
denom = np.zeros(model_fluxes.shape[1]) denom = np.zeros(model_fluxes.shape[1])
for i in range(obs_fluxes.size): for i in range(obs_fluxes.size):
if np.isfinite(obs_fluxes[i]) and obs_errors[i] > 0.: if np.isfinite(obs_fluxes[i]) and obs_fluxes_err[i] > 0.:
num += model_fluxes[i, :] * (obs_fluxes[i] / (obs_errors[i] * num += model_fluxes[i, :] * (obs_fluxes[i] / (obs_fluxes_err[i] *
obs_errors[i])) obs_fluxes_err[i]))
denom += np.square(model_fluxes[i, :] * (1./obs_errors[i])) denom += np.square(model_fluxes[i, :] * (1./obs_fluxes_err[i]))
for i in range(obs_propsmass.size):
if np.isfinite(obs_propsmass[i]) and obs_propsmass_err[i] > 0.:
num += model_propsmass[i, :] * (obs_propsmass[i] /
(obs_propsmass_err[i] *
obs_propsmass_err[i]))
denom += np.square(model_propsmass[i, :] *
(1./obs_propsmass_err[i]))
return num/denom return num/denom
def compute_chi2(model_fluxes, obs_fluxes, obs_errors, lim_flag): def compute_chi2(model_fluxes, model_props, model_propsmass, observation,
lim_flag):
"""Compute the χ² of observed fluxes with respect to the grid of models. We """Compute the χ² of observed fluxes with respect to the grid of models. We
take into account upper limits if need be. Note that we look over the bands take into account upper limits if need be. Note that we look over the bands
to avoid the creation of an array of the same size as the model_fluxes to avoid the creation of an array of the same size as the model_fluxes
...@@ -162,10 +177,13 @@ def compute_chi2(model_fluxes, obs_fluxes, obs_errors, lim_flag): ...@@ -162,10 +177,13 @@ def compute_chi2(model_fluxes, obs_fluxes, obs_errors, lim_flag):
---------- ----------
model_fluxes: array model_fluxes: array
2D grid containing the fluxes of the models 2D grid containing the fluxes of the models
obs_fluxes: array model_props: array
Fluxes of the observed object 2D grid containing the intensive properties of the models
obs_errors: array model_propsmass: array
Uncertainties on the fluxes of the observed object 2D grid containing the extensive properties of the models
observation: Class
Class instance containing the fluxes, intensive properties, extensive
properties and their errors, for a sigle observation.
lim_flag: boolean lim_flag: boolean
Boolean indicating whether upper limits should be treated (True) or Boolean indicating whether upper limits should be treated (True) or
discarded (False) discarded (False)
...@@ -177,34 +195,51 @@ def compute_chi2(model_fluxes, obs_fluxes, obs_errors, lim_flag): ...@@ -177,34 +195,51 @@ def compute_chi2(model_fluxes, obs_fluxes, obs_errors, lim_flag):
scaling: array scaling: array
scaling of the models to obtain the minimum χ² scaling of the models to obtain the minimum χ²
""" """
limits = lim_flag and np.any(obs_errors <= 0.) limits = lim_flag and np.any(observation.fluxes_err +
observation.extprops_err <= 0.)
scaling = _compute_scaling(model_fluxes, model_propsmass, observation)
scaling = _compute_scaling(model_fluxes, obs_fluxes, obs_errors) obs_fluxes = observation.fluxes
obs_fluxes_err = observation.fluxes_err
obs_props = observation.intprops
obs_props_err = observation.intprops_err
obs_propsmass = observation.extprops
obs_propsmass_err = observation.extprops_err
# Some observations may not have flux values in some filter(s), but # Some observations may not have flux values in some filter(s), but
# they can have upper limit(s). # they can have upper limit(s).
if limits == True: if limits == True:
obs_values = np.concatenate((obs_fluxes, obs_propsmass))
obs_values_err = np.concatenate((obs_fluxes_err, obs_propsmass_err))
model_values = np.concatenate((model_fluxes, model_propsmass))
for imod in range(scaling.size): for imod in range(scaling.size):
scaling[imod] = optimize.root(dchi2_over_ds2, scaling[imod], scaling[imod] = optimize.root(dchi2_over_ds2, scaling[imod],
args=(obs_fluxes, obs_errors, args=(obs_values, obs_values_err,
model_fluxes[:, imod])).x model_values[:, imod])).x
# χ² of the comparison of each model to each observation. # χ² of the comparison of each model to each observation.
chi2 = np.zeros(model_fluxes.shape[1]) chi2 = np.zeros(model_fluxes.shape[1])
for i in range(obs_fluxes.size): for i in range(obs_fluxes.size):
if np.isfinite(obs_fluxes[i]) and obs_errors[i] > 0.: if np.isfinite(obs_fluxes[i]) and obs_fluxes_err[i] > 0.:
chi2 += np.square( chi2 += np.square((obs_fluxes[i] - model_fluxes[i, :] * scaling) *
(obs_fluxes[i] - model_fluxes[i, :] * scaling) * (1./obs_errors[i])) (1./obs_fluxes_err[i]))
for i in range(obs_propsmass.size):
# Some observations may not have flux values in some filter(s), but if np.isfinite(obs_propsmass[i]):
chi2 += np.square((obs_propsmass[i] - model_propsmass[i, :] *
scaling) * (1./obs_propsmass_err[i]))
for i in range(obs_props.size):
if np.isfinite(obs_props[i]):
chi2 += np.square((obs_props[i] - model_props[i, :]) *
(1./obs_props_err[i]))
# they can have upper limit(s). # they can have upper limit(s).
if limits == True: if limits == True:
for i, obs_error in enumerate(obs_errors): for i, obs_error in enumerate(obs_fluxes_err):
if obs_error < 0.: if obs_error < 0.:
chi2 -= 2. * np.log(.5 * chi2 -= 2. * np.log(.5 *
(1. + erf(((obs_fluxes[i] - (1. + erf(((obs_fluxes[i] -
model_fluxes[i, :] * scaling) / model_fluxes[i, :] * scaling) /
(-np.sqrt(2.)*obs_errors[i]))))) (-np.sqrt(2.)*obs_fluxes_err[i])))))
return chi2, scaling return chi2, scaling
......
...@@ -126,12 +126,18 @@ def sed(idx, midx): ...@@ -126,12 +126,18 @@ def sed(idx, midx):
if 'sfh.age' in sed.info and sed.info['sfh.age'] > sed.info['universe.age']: if 'sfh.age' in sed.info and sed.info['sfh.age'] > sed.info['universe.age']:
gbl_models.fluxes[:, idx] = np.full(len(gbl_obs.bands), np.nan) gbl_models.fluxes[:, idx] = np.full(len(gbl_obs.bands), np.nan)
gbl_models.properties[:, idx] = np.full(len(gbl_properties), np.nan) gbl_models.properties[:, idx] = np.full(len(gbl_properties), np.nan)
gbl_models.intprops[:, idx] = np.full(len(gbl_obs.intprops), np.nan)
gbl_models.extprops[:, idx] = np.full(len(gbl_obs.extprops), np.nan)
else: else:
gbl_models.fluxes[:, idx] = [sed.compute_fnu(filter_) gbl_models.fluxes[:, idx] = [sed.compute_fnu(filter_)
for filter_ in gbl_obs.bands] for filter_ in gbl_obs.bands]
gbl_models.properties[:, idx] = [sed.info[name] gbl_models.properties[:, idx] = [sed.info[name]
for name in gbl_properties] for name in gbl_properties]
gbl_models.intprops[:, idx] = [sed.info[name]
for name in gbl_obs.intprops]
gbl_models.extprops[:, idx] = [sed.info[name]
for name in gbl_obs.extprops]
with gbl_ncomputed.get_lock(): with gbl_ncomputed.get_lock():
gbl_ncomputed.value += 1 gbl_ncomputed.value += 1
ncomputed = gbl_ncomputed.value ncomputed = gbl_ncomputed.value
...@@ -157,21 +163,22 @@ def analysis(idx, obs): ...@@ -157,21 +163,22 @@ def analysis(idx, obs):
""" """
np.seterr(invalid='ignore') np.seterr(invalid='ignore')
if obs['redshift'] >= 0.: if obs.redshift >= 0.:
# We pick the the models with the closest redshift using a slice to # We pick the the models with the closest redshift using a slice to
# work on views of the arrays and not on copies to save on RAM. # work on views of the arrays and not on copies to save on RAM.
z = np.array(gbl_models.conf['sed_modules_params']['redshifting']['redshift']) z = np.array(
wz = slice(np.abs(obs['redshift']-z).argmin(), None, z.size) gbl_models.conf['sed_modules_params']['redshifting']['redshift'])
corr_dz = compute_corr_dz(z[wz.start], obs['redshift']) wz = slice(np.abs(obs.redshift-z).argmin(), None, z.size)
corr_dz = compute_corr_dz(z[wz.start], obs.redshift)
else: # We do not know the redshift so we use the full grid else: # We do not know the redshift so we use the full grid
wz = slice(0, None, 1) wz = slice(0, None, 1)
corr_dz = 1. corr_dz = 1.
obs_fluxes = np.array([obs[name] for name in gbl_models.obs.bands]) observation = gbl_obs.observations[idx]
obs_errors = np.array([obs[name + "_err"] for name in chi2, scaling = compute_chi2(gbl_models.fluxes[:, wz],
gbl_models.obs.bands]) gbl_models.intprops[:, wz],
chi2, scaling = compute_chi2(gbl_models.fluxes[:, wz], obs_fluxes, gbl_models.extprops[:, wz], observation,
obs_errors, gbl_models.conf['analysis_params']['lim_flag']) gbl_models.conf['analysis_params']['lim_flag'])
if np.any(np.isfinite(chi2)): if np.any(np.isfinite(chi2)):
# We use the exponential probability associated with the χ² as # We use the exponential probability associated with the χ² as
...@@ -239,13 +246,15 @@ def bestfit(oidx, obs): ...@@ -239,13 +246,15 @@ def bestfit(oidx, obs):
# We compute the model at the exact redshift not to have to correct for the # We compute the model at the exact redshift not to have to correct for the
# difference between the object and the grid redshifts. # difference between the object and the grid redshifts.
params = deepcopy(gbl_params.from_index(best_index)) params = deepcopy(gbl_params.from_index(best_index))
params[-1]['redshift'] = obs['redshift'] params[-1]['redshift'] = obs.redshift
sed = gbl_warehouse.get_sed(gbl_params.modules, params) sed = gbl_warehouse.get_sed(gbl_params.modules, params)
fluxes = np.array([sed.compute_fnu(filt) for filt in gbl_obs.bands]) fluxes = np.array([sed.compute_fnu(filt) for filt in gbl_obs.bands])
obs_fluxes = np.array([obs[name] for name in gbl_obs.bands]) intprops = np.array([sed.info[prop] for prop in gbl_obs.intprops])
obs_errors = np.array([obs[name + '_err'] for name in gbl_obs.bands]) extprops = np.array([sed.info[prop] for prop in gbl_obs.extprops])
_, scaling = compute_chi2(fluxes[:, None], obs_fluxes, obs_errors,
_, scaling = compute_chi2(fluxes[:, None], intprops[:, None],
extprops[:, None], obs,
gbl_conf['analysis_params']['lim_flag']) gbl_conf['analysis_params']['lim_flag'])
gbl_results.best.properties[oidx, :] = [sed.info[k] for k in gbl_results.best.properties[oidx, :] = [sed.info[k] for k in
......
...@@ -81,18 +81,20 @@ class ObservationsManagerPassbands(object): ...@@ -81,18 +81,20 @@ class ObservationsManagerPassbands(object):
threshold) threshold)
self._add_model_error(modelerror) self._add_model_error(modelerror)
self.observations = list([Observation(row, self) for row in self.table])
def __len__(self): def __len__(self):
return len(self.table) return len(self.observations)
def __iter__(self): def __iter__(self):
self.idx = 0 self.idx = 0
self.max = len(self.table) self.max = len(self.observations)
return self return self
def __next__(self): def __next__(self):
if self.idx < self.max: if self.idx < self.max:
obs = self.table[self.idx] obs = self.observations[self.idx]
self.idx += 1 self.idx += 1
return obs return obs
raise StopIteration raise StopIteration
......
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