diff --git a/pcigale/analysis_modules/pdf_analysis/utils.py b/pcigale/analysis_modules/pdf_analysis/utils.py index fc78c32ecaba188701e1bb0644032a5db8a47a03..d39fe5ba90699b4f4aca140ff5cf076c9e1edb9c 100644 --- a/pcigale/analysis_modules/pdf_analysis/utils.py +++ b/pcigale/analysis_modules/pdf_analysis/utils.py @@ -6,7 +6,6 @@ # Author: Yannick Roehlly & Médéric Boquien from functools import lru_cache -import os from astropy import log from astropy.cosmology import WMAP7 as cosmo @@ -21,15 +20,12 @@ def save_chi2(obs, variable, models, chi2, values): """Save the chi² and the associated physocal properties """ - fname = 'out/{}_{}_chi2.npy'.format(obs['id'], variable) - if os.path.exists(fname): - data = np.memmap(fname, dtype=np.float64, mode='r+', - shape=(2, models.params.size)) - else: - data = np.memmap(fname, dtype=np.float64, mode='w+', - shape=(2, models.params.size)) - data[0, models.block] = chi2 - data[1, models.block] = values + fname = 'out/{}_{}_chi2-block-{}.npy'.format(obs['id'], variable, + models.iblock) + data = np.memmap(fname, dtype=np.float64, mode='w+', + shape=(2, chi2.size)) + data[0, :] = chi2 + data[1, :] = values @lru_cache(maxsize=None) diff --git a/pcigale/managers/models.py b/pcigale/managers/models.py index ec9e43831bd30bb6dc828b8f7a15355b3b8ef52b..228820cfa8ca7bc7fc319c660bd46ebe4e70ec41 100644 --- a/pcigale/managers/models.py +++ b/pcigale/managers/models.py @@ -30,6 +30,7 @@ class ModelsManager(object): self.conf = conf self.obs = obs self.params = params + self.iblock = iblock self.block = params.blocks[iblock] self.propertiesnames = conf['analysis_params']['variables'] diff --git a/pcigale_plots/__init__.py b/pcigale_plots/__init__.py index b2695ff0c4440fb9b80ebe822071070bbf80cee4..8254d413c7127d5465b489cac835da139e6b8dc7 100644 --- a/pcigale_plots/__init__.py +++ b/pcigale_plots/__init__.py @@ -7,6 +7,7 @@ # Author: Yannick Roehlly, Médéric Boquien & Denis Burgarella import argparse +import glob from itertools import product, repeat from collections import OrderedDict import sys @@ -49,23 +50,22 @@ def _chi2_worker(obj_name, var_name): Name of the analysed variable.. """ - fname = "out/{}_{}_chi2.npy".format(obj_name, var_name) - if os.path.isfile(fname): + figure = plt.figure() + ax = figure.add_subplot(111) + + fnames = glob.glob("out/{}_{}_chi2-block-*.npy".format(obj_name, var_name)) + for fname in fnames: data = np.memmap(fname, dtype=np.float64) data = np.memmap(fname, dtype=np.float64, shape=(2, data.size // 2)) - figure = plt.figure() - ax = figure.add_subplot(111) ax.scatter(data[1, :], data[0, :], color='k', s=.1) - ax.set_xlabel(var_name) - ax.set_ylabel("Reduced $\chi^2$") - ax.set_ylim(0., ) - ax.minorticks_on() - figure.suptitle("Reduced $\chi^2$ distribution of {} for {}." - .format(var_name, obj_name)) - figure.savefig("out/{}_{}_chi2.pdf".format(obj_name, var_name)) - plt.close(figure) - else: - print("No chi² found for {}. No plot created.".format(obj_name)) + ax.set_xlabel(var_name) + ax.set_ylabel("Reduced $\chi^2$") + ax.set_ylim(0., ) + ax.minorticks_on() + figure.suptitle("Reduced $\chi^2$ distribution of {} for {}." + .format(var_name, obj_name)) + figure.savefig("out/{}_{}_chi2.pdf".format(obj_name, var_name)) + plt.close(figure) def _pdf_worker(obj_name, var_name): @@ -79,43 +79,46 @@ def _pdf_worker(obj_name, var_name): Name of the analysed variable.. """ - fname = "out/{}_{}_chi2.npy".format(obj_name, var_name) - if os.path.isfile(fname): + + fnames = glob.glob("out/{}_{}_chi2-block-*.npy".format(obj_name, var_name)) + likelihood = [] + model_variable = [] + for fname in fnames: data = np.memmap(fname, dtype=np.float64) data = np.memmap(fname, dtype=np.float64, shape=(2, data.size // 2)) - likelihood = np.exp(-data[0, :] / 2.) - model_variable = data[1, :] + likelihood.append(np.exp(-data[0, :] / 2.)) + model_variable.append(data[1, :]) + likelihood = np.concatenate(likelihood) + model_variable = np.concatenate(model_variable) - Npdf = 100 - min_hist = np.min(model_variable) - max_hist = np.max(model_variable) - Nhist = min(Npdf, len(np.unique(model_variable))) + Npdf = 100 + min_hist = np.min(model_variable) + max_hist = np.max(model_variable) + Nhist = min(Npdf, len(np.unique(model_variable))) - if min_hist == max_hist: - pdf_grid = np.array([min_hist, max_hist]) - pdf_prob = np.array([1., 1.]) - else: - pdf_prob, pdf_grid = np.histogram(model_variable, Nhist, - (min_hist, max_hist), - weights=likelihood, density=True) - pdf_x = (pdf_grid[1:]+pdf_grid[:-1]) / 2. - - pdf_grid = np.linspace(min_hist, max_hist, Npdf) - pdf_prob = np.interp(pdf_grid, pdf_x, pdf_prob) - - figure = plt.figure() - ax = figure.add_subplot(111) - ax.plot(pdf_grid, pdf_prob, color='k') - ax.set_xlabel(var_name) - ax.set_ylabel("Probability density") - ax.minorticks_on() - figure.suptitle("Probability distribution function of {} for {}" - .format(var_name, obj_name)) - figure.savefig("out/{}_{}_pdf.pdf".format(obj_name, var_name)) - plt.close(figure) + if min_hist == max_hist: + pdf_grid = np.array([min_hist, max_hist]) + pdf_prob = np.array([1., 1.]) else: - print("No PDF found for {}. No plot created.".format(obj_name)) + pdf_prob, pdf_grid = np.histogram(model_variable, Nhist, + (min_hist, max_hist), + weights=likelihood, density=True) + pdf_x = (pdf_grid[1:]+pdf_grid[:-1]) / 2. + + pdf_grid = np.linspace(min_hist, max_hist, Npdf) + pdf_prob = np.interp(pdf_grid, pdf_x, pdf_prob) + + figure = plt.figure() + ax = figure.add_subplot(111) + ax.plot(pdf_grid, pdf_prob, color='k') + ax.set_xlabel(var_name) + ax.set_ylabel("Probability density") + ax.minorticks_on() + figure.suptitle("Probability distribution function of {} for {}" + .format(var_name, obj_name)) + figure.savefig("out/{}_{}_pdf.pdf".format(obj_name, var_name)) + plt.close(figure) def _sed_worker(obs, mod, filters, sed_type, nologo):