Skip to content

GitLab

Commit 11a1d3aa authored by Médéric Boquien's avatar Médéric Boquien
Browse files

Simplify and streamline the production of mock plots. Now just output one plot... 

Simplify and streamline the production of mock plots. Now just output one plot per file and do so in parallel. Also output in log when it seems to make sense.
parent c29ab4a9
Hide whitespace changes
Inline Side-by-side
pcigale_mock/__init__.py
View file @ 11a1d3aa
# -*- coding: utf-8 -*-
# Copyright (C) 2015 Laboratoire d'Astrophysique de Marseille
# Copyright (C) 2015 Universidad de Antofagasta
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Denis Burgarella
# Author: Médéric Boquien & Denis Burgarella
import argparse
from itertools import product, repeat
from collections import OrderedDict
import multiprocessing as mp
import os
import sys
from astropy.table import Table, join
from astropy.table import Table
import matplotlib
from matplotlib.backends.backend_pdf import PdfPages
#matplotlib.use('Agg')
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import multiprocessing as mp
import numpy as np
import os
import pkg_resources
from scipy.constants import c
from scipy import stats
from pcigale.data import Database
from pcigale.utils import read_table
from pcigale.session.configuration import Configuration
import matplotlib.gridspec as gridspec
__version__ = "0.1-alpha"
......@@ -32,157 +27,109 @@ MOCK_OUTPUT_FILE = "analysis_mock_results.txt"
# Directory where the output files are stored
OUT_DIR = "out/"
def worker(best_model, mock_output, mock_params, NRows, NCols, nologo):
"""Plot the diagnostic diagrammes from the mock analysis
def worker(exact, estimated, param, nologo):
"""Plot the exact and estimated values of a parameter for the mock analysis
Parameters
----------
best_model: Table row
Data from the input best_model file.
mock_output: Table row
Data from the out analysis of the mock catalogue
mock_params: Table row
Name of the parameter analysed that we want to plot
NRows: Integer
How many rows per page?
NCols: Integer
How many columns per page?
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.
"""
id_best = best_model['observation_id']
id_mock = mock_output['observation_id']
range_exact = np.linspace(np.min(exact), np.max(exact), 100)
if (len(id_best)>0 and len(id_mock)>0):
joined_table = join(left=best_model,right=mock_output,
keys='observation_id',
table_names=['Best', 'Mock'],
uniq_col_name='{table_name}_{col_name}')
else:
print("*** WARNING ***: No object found in {} or {}."
.format(OUT_DIR + BEST_MODEL_FILE, OUT_DIR + MOCK_OUTPUT_FILE))
NParams = len(mock_params)
if NParams <=0:
print("*** WARNING ***: No analysed parameter found ")
NPages = int(np.ceil(NParams / float(NRows*NCols)))
#print(NRows, NCols, NPages, NParams)
with PdfPages(OUT_DIR + 'mock.pdf') as pdf:
for Page in range(NPages):
fig = plt.figure()
plt.suptitle('Page '+str(Page))
gs = gridspec.GridSpec(NRows, NCols)
minParam = Page*NRows*NCols
maxParam = Page*NRows*NCols+NRows*NCols
for Param in enumerate(mock_params[minParam: maxParam]):
# We exclude outliers
mean_x = np.mean(best_model[Param[1]])
std_x = np.std(best_model[Param[1]])
if std_x > 0.:
mask = np.logical_and(mock_output[Param[1]] > mean_x-5.*std_x,
mock_output[Param[1]] < mean_x+5.*std_x
)
PercentExcluded = 100*(1.-np.sum(mask)/len(best_model[Param[1]]))
else:
mask = True
PercentExcluded = 0.0
# We compute the linear regression (excludind outliers)
if (np.min(best_model[Param[1]]) < np.max(best_model[Param[1]])):
slope, intercept, r_value, p_value, std_err = stats.linregress(
best_model[Param[1]][mask], mock_output[Param[1]][mask])
if p_value < 0.001:
p_value = 0.001
else:
slope=0.0
intercept=1.0
r_value=0.0
p_value=0.0
ax = fig.add_subplot(gs[Param[0]])
ax.set_xmargin(0.2)
ax.set_ymargin(0.2)
ax.errorbar(best_model[Param[1]][mask], mock_output[Param[1]][mask],
yerr=mock_output[Param[1]+'_err'][mask],
marker='.',
label=Param[1]+'\n'+ str("%.1f" % PercentExcluded+'% of outliers ($> 5\sigma$)'),
color='k', linestyle='None', capsize=0.)
ax.errorbar(best_model[Param[1]], best_model[Param[1]],
color='r', linestyle='-', label='1-to-1')
ax.errorbar(best_model[Param[1]], slope*best_model[Param[1]]+intercept,
color='b', linestyle='-',
label='best-fit'+' $r^2$ = '+ str("%.2f" % np.square(r_value))+', outliers not used')
ax.set_xlabel('Input', fontsize=10)
ax.set_ylabel('Output', fontsize=10)
ax.yaxis.labelpad = 0
if NRows > 2:
font_size = 4
else:
font_size = 6
ax.legend(fontsize=font_size, loc='best', fancybox=True, framealpha=0.5, numpoints=1)
ax.tick_params(axis='x', labelsize=8)
ax.tick_params(axis='y', labelsize=8)
if nologo is False:
image = plt.imread(pkg_resources.resource_filename(__name__,
"data/CIGALE.png"))
# Where do we plot CIGALE's logo?
if NRows == 1:
x0 = 225
y0 = 380
elif NRows == 2:
x0 = 300
y0 = 525
else:
x0 = 320
y0 = 525
fig.figimage(image, x0, y0, origin='upper', zorder=10,
alpha=1)
pdf.savefig() # saves the current figure into a pdf page
fig.tight_layout()
plt.close(fig)
def mock(config, NRowsNCols, nologo):
"""Plot the comparison of input/output values of analysed variables.
"""
if os.path.isfile(OUT_DIR + BEST_MODEL_FILE):
input = Table.read(OUT_DIR + BEST_MODEL_FILE, format='ascii')
# 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:
print("*** WARNING ***: No best model file found {}: error."
.format(OUT_DIR + BEST_MODEL_FILE))
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 nologo is False:
image = plt.imread(pkg_resources.resource_filename(__name__,
"data/CIGALE.png"))
plt.figimage(image, 510, 55, origin='upper', zorder=10, alpha=1)
plt.tight_layout()
plt.savefig(OUT_DIR + 'mock_{}.pdf'.format(param))
if np.all(exact > 0.) and range_exact[-1]/range_exact[0] > 20.:
plt.loglog()
plt.savefig(OUT_DIR + 'mock_log_{}.pdf'.format(param))
plt.close()
def mock(config, nologo):
"""Plot the comparison of input/output values of analysed variables in
parallel.
if os.path.isfile(OUT_DIR + MOCK_OUTPUT_FILE):
output = Table.read(OUT_DIR + MOCK_OUTPUT_FILE, format='ascii')
else:
print("*** WARNING ***: No best model file found {}."
.format(OUT_DIR + MOCK_OUTPUT_FILE))
Parameters
----------
config: configuration object
Contains the pcigale.ini configuration data
nologo: boolean
True when the CIGALE logo should not be drawn
analysed_params = config.configuration['analysis_method_params']
mock_params = (analysed_params['analysed_variables'])
"""
try:
exact = Table.read(OUT_DIR + BEST_MODEL_FILE, format='ascii')
except FileNotFoundError:
print("Best models file {} not found.".format(OUT_DIR +
BEST_MODEL_FILE))
sys.exit(1)
try:
estimated = Table.read(OUT_DIR + MOCK_OUTPUT_FILE, format='ascii')
except FileNotFoundError:
print("Mock models file {} not found.".format(OUT_DIR +
MOCK_OUTPUT_FILE))
sys.exit(1)
params = config.configuration['analysis_method_params']['analysed_variables']
arguments = ((exact[param], estimated[param], param, nologo) for param in
params)
with mp.Pool(processes=config.configuration['cores']) as pool:
pool.starmap(worker, arguments)
pool.close()
pool.join()
NRows = NRowsNCols
NCols = NRowsNCols
worker(input, output, mock_params, NRows, NCols, nologo)
def main():
parser = argparse.ArgumentParser(description='Create diagnotic plots for CIGALE')
if sys.version_info[:2] >= (3, 4):
mp.set_start_method('spawn')
else:
print("Could not set the multiprocessing start method to spawn. If "
"you encounter a deadlock, please upgrade to Python≥3.4.")
parser = argparse.ArgumentParser(description="Create diagnostic plots for "
"CIGALE")
parser.add_argument('--NRowsNCols', type=int, default='2',
help='How many columns/rows (1/1 or 2/2) per page in the mosaic?')
parser.add_argument('--nologo', action="store_true",
help='if you do not want the CIGALE logo in the figure')
help='Do not plot the CIGALE logo')
args = parser.parse_args()
config = Configuration()
mock(config, args.NRowsNCols, args.nologo)
mock(config, args.nologo)
Supports Markdown
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