Commit 444a153f authored by Médéric Boquien's avatar Médéric Boquien

As we loop over the bands to compute the χ², it is smarter that the bands are...

As we loop over the bands to compute the χ², it is smarter that the bands are the slower varying index and the models the faster varying index. That way it is more likely that models at the correct redshift are prefetched into the cache.
parent fb0570e1
......@@ -10,6 +10,7 @@
### Optimised
- A significant fraction of the total run time is spent computing integrals (e.g. fluxes in passbands). We can make the integration faster by rewriting the trapezoidal rule in terms of np.dot(). This allows to offload the computation to optimised libraries. The end result is that the integration is twice as fast, with a gain of ~10-15% on the total run time. (Médéric Boquien)
- The conversion from luminosity to flux is now a bit faster. (Médéric Boquien)
- The order the models are stored in memory has been changed to make the computation of the χ² faster. (Médéric Boquien)
## 0.9.0 (2016-04-04)
### Added
......
......@@ -161,9 +161,9 @@ class PdfAnalysis(AnalysisModule):
# We put the shape in a tuple along with the RawArray because workers
# need to know the shape to create the numpy array from the RawArray.
model_fluxes = (RawArray(ctypes.c_double, n_params * n_filters),
(n_params, n_filters))
(n_filters, n_params))
model_variables = (RawArray(ctypes.c_double, n_params * n_variables),
(n_params, n_variables))
(n_variables, n_params))
initargs = (params, filters, variables_nolog, model_fluxes,
model_variables, time.time(), mp.Value('i', 0))
......
......@@ -102,16 +102,16 @@ def save_pdf(obsid, names, mass_proportional, model_variables, scaling,
for i, name in enumerate(names):
if name.endswith('_log'):
if name[:-4] in mass_proportional:
model_variable = np.log10(model_variables[:, i][wlikely] *
model_variable = np.log10(model_variables[i, :][wlikely] *
scaling[wlikely])
else:
model_variable = np.log10(model_variables[:, i][wlikely])
model_variable = np.log10(model_variables[i, :][wlikely])
else:
if name in mass_proportional:
model_variable = (model_variables[:, i][wlikely] *
model_variable = (model_variables[i, :][wlikely] *
scaling[wlikely])
else:
model_variable = model_variables[:, i][wlikely]
model_variable = model_variables[i, :][wlikely]
_save_pdf(obsid, name, model_variable, likelihood)
......@@ -155,14 +155,14 @@ def save_chi2(obsid, names, mass_proportional, model_variables, scaling, chi2):
for i, name in enumerate(names):
if name.endswith('_log'):
if name[:-4] in mass_proportional:
model_variable = np.log10(model_variables[:, i] * scaling)
model_variable = np.log10(model_variables[i, :] * scaling)
else:
model_variable = np.log10(model_variables[:, i])
model_variable = np.log10(model_variables[i, :])
else:
if name in mass_proportional:
model_variable = model_variables[:, i] * scaling
model_variable = model_variables[i, :] * scaling
else:
model_variable = model_variables[:, i]
model_variable = model_variables[i, :]
_save_chi2(obsid, name, model_variable, chi2)
......@@ -344,13 +344,13 @@ def _compute_scaling(model_fluxes, obs_fluxes, obs_errors):
scaling: array
Scaling factors minimising the χ²
"""
num = np.zeros(model_fluxes.shape[0])
denom = np.zeros(model_fluxes.shape[0])
num = np.zeros(model_fluxes.shape[1])
denom = np.zeros(model_fluxes.shape[1])
for i in range(obs_fluxes.size):
if np.isfinite(obs_fluxes[i]):
num += model_fluxes[:, i] * (obs_fluxes[i] / (obs_errors[i] *
num += model_fluxes[i, :] * (obs_fluxes[i] / (obs_errors[i] *
obs_errors[i]))
denom += np.square(model_fluxes[:, i] / obs_errors[i])
denom += np.square(model_fluxes[i, :] * (1./obs_errors[i]))
return num/denom
......@@ -390,14 +390,14 @@ def compute_chi2(model_fluxes, obs_fluxes, obs_errors, lim_flag):
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
model_fluxes[:, imod])).x
# χ² of the comparison of each model to each observation.
chi2 = np.zeros(model_fluxes.shape[0])
chi2 = np.zeros(model_fluxes.shape[1])
for i in range(obs_fluxes.size):
if np.isfinite(obs_fluxes[i]) and obs_errors[i] > 0.:
chi2 += np.square(
(obs_fluxes[i] - model_fluxes[:, i] * scaling) / obs_errors[i])
(obs_fluxes[i] - model_fluxes[i, :] * scaling) * (1./obs_errors[i]))
# Some observations may not have flux values in some filter(s), but
# they can have upper limit(s).
......@@ -407,7 +407,7 @@ def compute_chi2(model_fluxes, obs_fluxes, obs_errors, lim_flag):
np.log(
np.sqrt(np.pi/2.)*(-obs_errors[mask_lim])*(
1.+erf(
(obs_fluxes[mask_lim]-model_fluxes[:, mask_lim] *
(obs_fluxes[mask_lim]-model_fluxes[mask_lim, :] *
scaling[:, np.newaxis]) /
(np.sqrt(2)*(-obs_errors[mask_lim]))))), axis=1)
......
......@@ -158,13 +158,13 @@ def sed(idx):
gbl_params.from_index(idx))
if 'sfh.age' in sed.info and sed.info['sfh.age'] > sed.info['universe.age']:
gbl_model_fluxes[idx, :] = np.full(len(gbl_filters), np.nan)
gbl_model_variables[idx, :] = np.full(len(gbl_analysed_variables),
gbl_model_fluxes[:, idx] = np.full(len(gbl_filters), np.nan)
gbl_model_variables[:, idx] = np.full(len(gbl_analysed_variables),
np.nan)
else:
gbl_model_fluxes[idx, :] = np.array([sed.compute_fnu(filter_) for
gbl_model_fluxes[:, idx] = np.array([sed.compute_fnu(filter_) for
filter_ in gbl_filters])
gbl_model_variables[idx, :] = np.array([sed.info[name]
gbl_model_variables[:, idx] = np.array([sed.info[name]
for name in
gbl_analysed_variables])
......@@ -225,7 +225,7 @@ def analysis(idx, obs):
wz = slice(0, None, 1)
corr_dz = 1.
chi2, scaling = compute_chi2(gbl_model_fluxes[wz, :], obs_fluxes,
chi2, scaling = compute_chi2(gbl_model_fluxes[:, wz], obs_fluxes,
obs_errors, gbl_lim_flag)
##################################################################
......@@ -285,17 +285,17 @@ def analysis(idx, obs):
maxstd = lambda mean, std: max(0.05 * mean, std)
if variable in sed.mass_proportional_info:
mean, std = weighted_param(_(gbl_model_variables[wz, i][wlikely]
mean, std = weighted_param(_(gbl_model_variables[i, wz][wlikely]
* scaling[wlikely] * corr_dz),
likelihood)
else:
mean, std = weighted_param(_(gbl_model_variables[wz, i][wlikely]),
mean, std = weighted_param(_(gbl_model_variables[i, wz][wlikely]),
likelihood)
gbl_analysed_averages[idx, i] = mean
gbl_analysed_std[idx, i] = maxstd(mean, std)
gbl_best_fluxes[idx, :] = gbl_model_fluxes[best_index, :] \
gbl_best_fluxes[idx, :] = gbl_model_fluxes[:, best_index] \
* scaling[best_index_z]
global gbl_keys
......@@ -310,11 +310,11 @@ def analysis(idx, obs):
save_best_sed(obs['id'], sed, scaling[best_index_z])
if gbl_save['chi2']:
save_chi2(obs['id'], gbl_analysed_variables,
sed.mass_proportional_info, gbl_model_variables[wz, :],
sed.mass_proportional_info, gbl_model_variables[:, wz],
scaling * corr_dz, chi2 / (nobs - 1))
if gbl_save['pdf']:
save_pdf(obs['id'], gbl_analysed_variables,
sed.mass_proportional_info, gbl_model_variables[wz, :],
sed.mass_proportional_info, gbl_model_variables[:, wz],
scaling * corr_dz, likelihood, wlikely)
with gbl_n_computed.get_lock():
......
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