From eedcad3a45f75ccde6593fb5b8c22d04fd3cefbb Mon Sep 17 00:00:00 2001 From: = Date: Fri, 14 Dec 2018 13:46:52 -0300 Subject: [PATCH] Change from __init__.py to own module. --- pcigale_plots/parsers/chi2.py | 56 +++++++++++++++++++++ pcigale_plots/parsers/mock.py | 92 +++++++++++++++++++++++++++++++++++ pcigale_plots/parsers/pdf.py | 90 ++++++++++++++++++++++++++++++++++ 3 files changed, 238 insertions(+) create mode 100644 pcigale_plots/parsers/chi2.py create mode 100644 pcigale_plots/parsers/mock.py create mode 100644 pcigale_plots/parsers/pdf.py diff --git a/pcigale_plots/parsers/chi2.py b/pcigale_plots/parsers/chi2.py new file mode 100644 index 000000000..f787a2d39 --- /dev/null +++ b/pcigale_plots/parsers/chi2.py @@ -0,0 +1,56 @@ + +import glob +from itertools import product + +import matplotlib + +matplotlib.use('Agg') +import matplotlib.pyplot as plt +import multiprocessing as mp +import numpy as np +from pcigale.utils import read_table + + +def chi2(config): + """Plot the χ² values of analysed variables. + """ + input_data = read_table(config.configuration['data_file']) + chi2_vars = config.configuration['analysis_params']['variables'] + chi2_vars += [band for band in config.configuration['bands'] + if band.endswith('_err') is False] + + with mp.Pool(processes=config.configuration['cores']) as pool: + items = product(input_data['id'], chi2_vars) + pool.starmap(_chi2_worker, items) + pool.close() + pool.join() + + +def _chi2_worker(obj_name, var_name): + """Plot the reduced χ² associated with a given analysed variable + + Parameters + ---------- + obj_name: string + Name of the object. + var_name: string + Name of the analysed variable.. + + """ + figure = plt.figure() + ax = figure.add_subplot(111) + + var_name = var_name.replace('/', '_') + 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)) + 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) diff --git a/pcigale_plots/parsers/mock.py b/pcigale_plots/parsers/mock.py new file mode 100644 index 000000000..0649da008 --- /dev/null +++ b/pcigale_plots/parsers/mock.py @@ -0,0 +1,92 @@ +from astropy.table import Table +import matplotlib +import sys + +matplotlib.use('Agg') +import matplotlib.pyplot as plt +import multiprocessing as mp +import numpy as np +import pkg_resources +from scipy import stats + + +def mock(config, mock_results_file, best_results_file, nologo): + """Plot the comparison of input/output values of analysed variables. + """ + + try: + exact = Table.read(best_results_file) + except FileNotFoundError: + print("Best models file {} not found.".format(best_results_file)) + sys.exit(1) + + try: + estimated = Table.read(mock_results_file) + except FileNotFoundError: + print("Mock models file {} not found.".format(mock_results_file)) + sys.exit(1) + + params = config.configuration['analysis_params']['variables'] + + for param in params: + if param.endswith('_log'): + param = "best."+param + exact[param] = np.log10(exact[param[:-4]]) + + logo = False if nologo else plt.imread(pkg_resources.resource_filename(__name__, + "../resources/CIGALE.png")) + + arguments = ((exact["best."+param], estimated["bayes."+param], param, logo) + for param in params) + + with mp.Pool(processes=config.configuration['cores']) as pool: + pool.starmap(_mock_worker, arguments) + pool.close() + pool.join() + + +def _mock_worker(exact, estimated, param, logo): + """Plot the exact and estimated values of a parameter for the mock analysis + + Parameters + ---------- + exact: Table column + Exact values of the parameter. + estimated: Table column + Estimated values of the parameter. + param: string + Name of the parameter + nologo: boolean + Do not add the logo when set to true. + + """ + + range_exact = np.linspace(np.min(exact), np.max(exact), 100) + + # We compute the linear regression + if np.min(exact) < np.max(exact): + slope, intercept, r_value, p_value, std_err = stats.linregress(exact, + estimated) + else: + slope = 0.0 + intercept = 1.0 + r_value = 0.0 + + plt.errorbar(exact, estimated, marker='.', label=param, color='k', + linestyle='None', capsize=0.) + plt.plot(range_exact, range_exact, color='r', label='1-to-1') + plt.plot(range_exact, slope * range_exact + intercept, color='b', + label='exact-fit $r^2$ = {:.2f}'.format(r_value**2)) + plt.xlabel('Exact') + plt.ylabel('Estimated') + plt.title(param) + plt.legend(loc='best', fancybox=True, framealpha=0.5, numpoints=1) + plt.minorticks_on() + + if logo is not False: + plt.figimage(logo, 510, 55, origin='upper', zorder=10, alpha=1) + + plt.tight_layout() + plt.savefig('out/mock_{}.pdf'.format(param)) + + plt.close() diff --git a/pcigale_plots/parsers/pdf.py b/pcigale_plots/parsers/pdf.py new file mode 100644 index 000000000..6b5832a6c --- /dev/null +++ b/pcigale_plots/parsers/pdf.py @@ -0,0 +1,90 @@ + +import glob +from itertools import product +import matplotlib + +matplotlib.use('Agg') +import matplotlib.pyplot as plt +import multiprocessing as mp +import numpy as np +from pcigale.utils import read_table + + +def pdf(config): + """Plot the PDF of analysed variables. + """ + input_data = read_table(config.configuration['data_file']) + pdf_vars = config.configuration['analysis_params']['variables'] + pdf_vars += [band for band in config.configuration['bands'] + if band.endswith('_err') is False] + + with mp.Pool(processes=config.configuration['cores']) as pool: + items = product(input_data['id'], pdf_vars) + pool.starmap(_pdf_worker, items) + pool.close() + pool.join() + + +def _pdf_worker(obj_name, var_name): + """Plot the PDF associated with a given analysed variable + + Parameters + ---------- + obj_name: string + Name of the object. + var_name: string + Name of the analysed variable.. + + """ + var_name = var_name.replace('/', '_') + if var_name.endswith('_log'): + fnames = glob.glob("out/{}_{}_chi2-block-*.npy".format(obj_name, + var_name[:-4])) + log = True + else: + fnames = glob.glob("out/{}_{}_chi2-block-*.npy".format(obj_name, + var_name)) + log = False + 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.append(np.exp(-data[0, :] / 2.)) + model_variable.append(data[1, :]) + likelihood = np.concatenate(likelihood) + model_variable = np.concatenate(model_variable) + if log is True: + model_variable = np.log10(model_variable) + w = np.where(np.isfinite(likelihood) & np.isfinite(model_variable)) + likelihood = likelihood[w] + model_variable = model_variable[w] + + 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) -- 2.26.2