From a444d8deb02f268e84dbd4cb9ad54047f2e7386b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A9d=C3=A9ric=20Boquien?= Date: Mon, 11 Dec 2017 12:42:19 -0300 Subject: [PATCH] =?UTF-8?q?For=20more=20reliability=20we=20now=20save=20th?= =?UTF-8?q?e=20chi=C2=B2=20by=20blocks.=20Otherwise=20it=20could=20lead=20?= =?UTF-8?q?to=20some=20crashes=20when=20there=20was=20more=20than=201=20re?= =?UTF-8?q?dshift.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../analysis_modules/pdf_analysis/utils.py | 16 ++-- pcigale/managers/models.py | 1 + pcigale_plots/__init__.py | 93 ++++++++++--------- 3 files changed, 55 insertions(+), 55 deletions(-) diff --git a/pcigale/analysis_modules/pdf_analysis/utils.py b/pcigale/analysis_modules/pdf_analysis/utils.py index fc78c32e..d39fe5ba 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 ec9e4383..228820cf 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 b2695ff0..8254d413 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): -- GitLab