Commit 475a3c04 authored by Médéric Boquien's avatar Médéric Boquien

Do not access SharedArray.array directly anymore but rather use the fact that...

Do not access SharedArray.array directly anymore but rather use the fact that SharedArray behaves somewhat like a Numpy array for basical operations.
parent deb1940e
......@@ -144,15 +144,15 @@ def _compute_scaling(models, obs, wz):
"""
_ = list(models.flux.keys())[0]
num = np.zeros_like(models.flux[_].array[wz])
denom = np.zeros_like(models.flux[_].array[wz])
num = np.zeros_like(models.flux[_][wz])
denom = np.zeros_like(models.flux[_][wz])
for band, flux in obs.flux.items():
# Multiplications are faster than division, so we directly use the
# inverse error
inv_err2 = 1. / obs.flux_err[band] ** 2.
if np.isfinite(flux) and obs.flux_err[band] > 0.:
model = models.flux[band].array[wz]
model = models.flux[band][wz]
num += model * (flux * inv_err2)
denom += model**2. * inv_err2
......@@ -161,7 +161,7 @@ def _compute_scaling(models, obs, wz):
# inverse error
inv_err2 = 1. / obs.extprop_err[name] ** 2.
if np.isfinite(prop) and obs.extprop_err[name] > 0.:
model = models.extprop[name].array[wz]
model = models.extprop[name][wz]
num += model * (prop * inv_err2)
denom += model ** 2. * inv_err2
......@@ -223,13 +223,13 @@ def compute_chi2(models, obs, corr_dz, wz, lim_flag):
# inverse error
inv_flux_err = 1. / obs.flux_err[band]
if np.isfinite(flux) and inv_flux_err > 0.:
model = models.flux[band].array[wz]
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].array[wz]
model = models.intprop[name][wz]
chi2 += ((prop - model) * (1. / obs.intprops_err[name])) ** 2.
# Computation of the χ² from extensive properties
......@@ -238,7 +238,7 @@ def compute_chi2(models, obs, corr_dz, wz, lim_flag):
# inverse error
inv_prop_err = 1. / obs.extprop_err[name]
if np.isfinite(prop) and inv_prop_err > 0.:
model = models.extprop[name].array[wz]
model = models.extprop[name][wz]
chi2 += ((prop - (scaling * model) * corr_dz) * inv_prop_err) ** 2.
# Finally take the presence of upper limits into account
......@@ -247,7 +247,7 @@ def compute_chi2(models, obs, corr_dz, wz, lim_flag):
if obs_error < 0.:
chi2 -= 2. * np.log(.5 *
(1. + erf(((obs.fluxes[band] -
model_flux[band].array[wz] * scaling) /
model_flux[band][wz] * scaling) /
(-np.sqrt(2.)*obs_error)))))
return chi2, scaling
......
......@@ -124,19 +124,19 @@ def sed(idx, midx):
if 'sfh.age' in sed.info and sed.info['sfh.age'] > sed.info['universe.age']:
for band in gbl_models.flux:
gbl_models.flux[band].array[idx] = np.nan
gbl_models.flux[band][idx] = np.nan
for prop in gbl_models.extprop:
gbl_models.extprop[prop].array[idx] = np.nan
gbl_models.extprop[prop][idx] = np.nan
for prop in gbl_models.intprop:
gbl_models.intprop[prop].array[idx] = np.nan
gbl_models.intprop[prop][idx] = np.nan
else:
for band in gbl_models.flux.keys():
gbl_models.flux[band].array[idx] = sed.compute_fnu(band)
gbl_models.flux[band][idx] = sed.compute_fnu(band)
for prop in gbl_models.extprop.keys():
gbl_models.extprop[prop].array[idx] = sed.info[prop]
gbl_models.extprop[prop][idx] = sed.info[prop]
for prop in gbl_models.intprop.keys():
gbl_models.intprop[prop].array[idx] = sed.info[prop]
gbl_models.intprop[prop][idx] = sed.info[prop]
with gbl_ncomputed.get_lock():
gbl_ncomputed.value += 1
......@@ -186,7 +186,7 @@ def analysis(idx, obs):
# If all the models are valid, it is much more efficient to use a slice
if likelihood.size == wlikely[0].size:
wlikely = slice(None, None)
gbl_results.bayes.weight.array[idx] = np.nansum(likelihood)
gbl_results.bayes.weight[idx] = np.nansum(likelihood)
# We compute the weighted average and standard deviation using the
# likelihood as weight.
......@@ -196,10 +196,10 @@ def analysis(idx, obs):
_ = np.log10
else:
_ = lambda x: x
values = _(gbl_models.intprop[prop].array[wz])
values = _(gbl_models.intprop[prop][wz])
mean, std = weighted_param(values[wlikely], likelihood[wlikely])
gbl_results.bayes.intmean[prop].array[idx] = mean
gbl_results.bayes.interror[prop].array[idx] = std
gbl_results.bayes.intmean[prop][idx] = mean
gbl_results.bayes.interror[prop][idx] = std
if gbl_models.conf['analysis_params']['save_chi2'] is True:
save_chi2(obs, prop, gbl_models, chi2, values)
......@@ -209,11 +209,11 @@ def analysis(idx, obs):
_ = np.log10
else:
_ = lambda x: x
values = _(gbl_models.extprop[prop].array[wz])
values = _(gbl_models.extprop[prop][wz])
mean, std = weighted_param(values[wlikely] * scaling * corr_dz,
likelihood[wlikely])
gbl_results.bayes.extmean[prop].array[idx] = mean
gbl_results.bayes.exterror[prop].array[idx] = std
gbl_results.bayes.extmean[prop][idx] = mean
gbl_results.bayes.exterror[prop][idx] = std
if gbl_models.conf['analysis_params']['save_chi2'] is True:
save_chi2(obs, prop, gbl_models, chi2, values)
......@@ -269,13 +269,13 @@ def bestfit(oidx, obs):
scaling = gbl_results.best.scaling[oidx]
for band in gbl_results.best.flux:
gbl_results.best.flux[band].array[oidx] = sed.compute_fnu(band) * scaling
gbl_results.best.flux[band][oidx] = sed.compute_fnu(band) * scaling
for prop in gbl_results.best.intprop:
gbl_results.best.intprop[prop].array[oidx] = sed.info[prop]
gbl_results.best.intprop[prop][oidx] = sed.info[prop]
for prop in gbl_results.best.extprop:
gbl_results.best.extprop[prop].array[oidx] = sed.info[prop] * scaling \
gbl_results.best.extprop[prop][oidx] = sed.info[prop] * scaling \
* corr_dz
if gbl_conf['analysis_params']["save_best_sed"]:
......
......@@ -74,11 +74,11 @@ def fluxes(idx, midx):
gbl_models.intprop[prop][idx] = np.nan
else:
for band in gbl_models.flux.keys():
gbl_models.flux[band].array[idx] = sed.compute_fnu(band)
gbl_models.flux[band][idx] = sed.compute_fnu(band)
for prop in gbl_models.extprop.keys():
gbl_models.extprop[prop].array[idx] = sed.info[prop]
gbl_models.extprop[prop][idx] = sed.info[prop]
for prop in gbl_models.intprop.keys():
gbl_models.intprop[prop].array[idx] = sed.info[prop]
gbl_models.intprop[prop][idx] = sed.info[prop]
if gbl_save is True:
sed.to_fits("out/{}".format(midx))
......
......@@ -83,12 +83,12 @@ class ModelsManager(object):
table = Table()
table.add_column(Column(self.block, name='id'))
for band in sorted(self.flux.keys()):
table.add_column(Column(self.flux[band].array, name=band,
table.add_column(Column(self.flux[band], name=band,
unit='mJy'))
for prop in sorted(self.extprop.keys()):
table.add_column(Column(self.extprop[prop].array, name=prop))
table.add_column(Column(self.extprop[prop], name=prop))
for prop in sorted(self.intprop.keys()):
table.add_column(Column(self.intprop[prop].array, name=prop))
table.add_column(Column(self.intprop[prop], name=prop))
table.write("out/{}.fits".format(filename))
table.write("out/{}.txt".format(filename), format='ascii.fixed_width',
......
......@@ -86,24 +86,24 @@ class BayesResultsManager(object):
raise TypeError("A list of BayesResultsManager is required.")
merged = results[0]
intmean = {prop: np.array([result.intmean[prop].array
intmean = {prop: np.array([result.intmean[prop]
for result in results])
for prop in merged.intmean}
interror = {prop: np.array([result.interror[prop].array
interror = {prop: np.array([result.interror[prop]
for result in results])
for prop in merged.interror}
extmean = {prop: np.array([result.extmean[prop].array
extmean = {prop: np.array([result.extmean[prop]
for result in results])
for prop in merged.extmean}
exterror = {prop: np.array([result.exterror[prop].array
exterror = {prop: np.array([result.exterror[prop]
for result in results])
for prop in merged.exterror}
weight = np.array([result.weight.array for result in results])
weight = np.array([result.weight for result in results])
totweight = np.sum(weight, axis=0)
for prop in merged.intmean:
merged.intmean[prop].array[:] = np.sum(
merged.intmean[prop][:] = np.sum(
intmean[prop] * weight, axis=0) / totweight
# We compute the merged standard deviation by combining the
......@@ -111,11 +111,11 @@ class BayesResultsManager(object):
# http://stats.stackexchange.com/a/10445 where the number of
# datapoints has been substituted with the weights. In short we
# exploit the fact that Var(X) = E(Var(X)) + Var(E(X)).
merged.interror[prop].array[:] = np.sqrt(np.sum(
weight * (interror[prop]**2. + (intmean[prop]-merged.intmean[prop].array)**2), axis=0) / totweight)
merged.interror[prop][:] = np.sqrt(np.sum(
weight * (interror[prop]**2. + (intmean[prop]-merged.intmean[prop])**2), axis=0) / totweight)
for prop in merged.extmean:
merged.extmean[prop].array[:] = np.sum(
merged.extmean[prop][:] = np.sum(
extmean[prop] * weight, axis=0) / totweight
# We compute the merged standard deviation by combining the
......@@ -123,17 +123,17 @@ class BayesResultsManager(object):
# http://stats.stackexchange.com/a/10445 where the number of
# datapoints has been substituted with the weights. In short we
# exploit the fact that Var(X) = E(Var(X)) + Var(E(X)).
merged.exterror[prop].array[:] = np.sqrt(np.sum(
weight * (exterror[prop]**2. + (extmean[prop]-merged.extmean[prop].array)**2), axis=0) / totweight)
merged.exterror[prop][:] = np.sqrt(np.sum(
weight * (exterror[prop]**2. + (extmean[prop]-merged.extmean[prop])**2), axis=0) / totweight)
for prop in merged.extmean:
if prop.endswith('_log'):
merged.exterror[prop].array[:] = \
np.maximum(0.02, merged.exterror[prop].array)
merged.exterror[prop][:] = \
np.maximum(0.02, merged.exterror[prop])
else:
merged.exterror[prop].array[:] = \
np.maximum(0.05 * merged.extmean[prop].array,
merged.exterror[prop].array)
merged.exterror[prop][:] = \
np.maximum(0.05 * merged.extmean[prop],
merged.exterror[prop])
return merged
......@@ -210,7 +210,7 @@ class BestResultsManager(object):
each observation.
"""
return self._chi2.array
return self._chi2
@chi2.setter
def chi2(self, chi2):
......@@ -222,7 +222,7 @@ class BestResultsManager(object):
observation.
"""
return self._index.array
return self._index
@index.setter
def index(self, index):
......@@ -234,7 +234,7 @@ class BestResultsManager(object):
observation.
"""
return self._scaling.array
return self._scaling
@scaling.setter
def scaling(self, scaling):
......@@ -264,14 +264,14 @@ class BestResultsManager(object):
merged = results[0]
for iobs, bestidx in enumerate(best):
for band in merged.flux:
merged.flux[band].array[iobs] = \
results[bestidx].flux[band].array[iobs]
merged.flux[band][iobs] = \
results[bestidx].flux[band][iobs]
for prop in merged.intprop:
merged.intprop[prop].array[iobs] = \
results[bestidx].intprop[prop].array[iobs]
merged.intprop[prop][iobs] = \
results[bestidx].intprop[prop][iobs]
for prop in merged.extprop:
merged.extprop[prop].array[iobs] = \
results[bestidx].extprop[prop].array[iobs]
merged.extprop[prop][iobs] = \
results[bestidx].extprop[prop][iobs]
merged.chi2[iobs] = results[bestidx].chi2[iobs]
merged.scaling[iobs] = results[bestidx].scaling[iobs]
merged.index[iobs] = results[bestidx].index[iobs]
......@@ -343,14 +343,14 @@ class ResultsManager(object):
table.add_column(Column(self.obs.table['id'], name="id"))
for prop in sorted(self.bayes.intmean):
table.add_column(Column(self.bayes.intmean[prop].array,
table.add_column(Column(self.bayes.intmean[prop],
name="bayes."+prop))
table.add_column(Column(self.bayes.interror[prop].array,
table.add_column(Column(self.bayes.interror[prop],
name="bayes."+prop+"_err"))
for prop in sorted(self.bayes.extmean):
table.add_column(Column(self.bayes.extmean[prop].array,
table.add_column(Column(self.bayes.extmean[prop],
name="bayes."+prop))
table.add_column(Column(self.bayes.exterror[prop].array,
table.add_column(Column(self.bayes.exterror[prop],
name="bayes."+prop+"_err"))
table.add_column(Column(self.best.chi2, name="best.chi_square"))
......@@ -360,14 +360,14 @@ class ResultsManager(object):
name="best.reduced_chi_square"))
for prop in sorted(self.best.intprop):
table.add_column(Column(self.best.intprop[prop].array,
table.add_column(Column(self.best.intprop[prop],
name="best."+prop))
for prop in sorted(self.best.extprop):
table.add_column(Column(self.best.extprop[prop].array,
table.add_column(Column(self.best.extprop[prop],
name="best."+prop))
for band in self.obs.bands:
table.add_column(Column(self.best.flux[band].array,
table.add_column(Column(self.best.flux[band],
name="best."+band, unit='mJy'))
......
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