Commit df30a6e3 authored by Médéric Boquien's avatar Médéric Boquien

Merge pcigale-mock into pcigale-plots.

parent 69338e2b
......@@ -9,11 +9,13 @@
- To homogenize input and output files, the "observation_id" has been changed to "id" in the output files. (Médéric Boquien)
- The output files providing estimates of the physical properties are now generated both in the form of text and FITS files. (Médéric Boquien)
- When using the dustatt_calzleit module, choosing ẟ≠0 leads to an effective E(B-V) different from the one set by the user. Now the E(B-V) will always correspond to the one specified by the user. This means that at fixed E(B-V), A(V) depends on ẟ. (Médéric Boquien)
- The pcigale-mock tool has been merged into pcigale-plots; the mock plots can be obtained with the "mock" command.
### Fixed
- To estimate parameters in log, pcigale determines which variables end with the "_log" string and removed it to compute the models. However in some circumstances, it was overzealous. This has been fixed. (Médéric Boquien)
- When estimating a parameter in log, these were not scaled appropriately and taken in log when saving the related χ² and PDF. (Médéric Boquien)
- In the presence of upper limits, correct the scaling factor of the models to the observations before computing the χ², not after. (Médéric Boquien)
- When called without arguments, pcigale-plots would crash and display the backtrace. Now it displays the a short help showing how to use it. (Médéric Boquien)
### Optimised
- Prior to version 0.7.0, we needed to maintain the list of redshifts for all the computed models. Past 0.7.0 we just infer the redshift from a list unique redshifts. This means that we can now discard the list of redshifts for all the models and only keep the list of unique redshifts. This saves ~8 MB of memory for every 10⁶ models. the models should be computed slightly faster but it is in the measurement noise. (Médéric Boquien)
......
# -*- 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: Médéric Boquien & Denis Burgarella
import argparse
import multiprocessing as mp
import os
import sys
from astropy.table import Table
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
import pkg_resources
from scipy import stats
from pcigale.session.configuration import Configuration
__version__ = "0.1-alpha"
# Name of the file containing the best models information
BEST_MODEL_FILE = "results.fits"
MOCK_OUTPUT_FILE = "results_mock.fits"
# Directory where the output files are stored
OUT_DIR = "out/"
def worker(exact, estimated, param, nologo):
"""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 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.
Parameters
----------
config: configuration object
Contains the pcigale.ini configuration data
nologo: boolean
True when the CIGALE logo should not be drawn
"""
try:
exact = Table.read(OUT_DIR + BEST_MODEL_FILE)
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)
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()
def main():
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('--nologo', action="store_true",
help='Do not plot the CIGALE logo')
args = parser.parse_args()
config = Configuration()
mock(config, args.nologo)
......@@ -21,6 +21,7 @@ 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
......@@ -30,7 +31,8 @@ __version__ = "0.1-alpha"
# Name of the file containing the best models information
BEST_MODEL_FILE = "results.fits"
BEST_RESULTS = "results.fits"
MOCK_RESULTS = "results_mock.fits"
# Directory where the output files are stored
OUT_DIR = "out/"
# Wavelength limits (restframe) when plotting the best SED.
......@@ -316,6 +318,54 @@ def _sed_worker(obs, mod, filters, sed_type, nologo):
print("No SED found for {}. No plot created.".format(obs['id']))
def _mock_worker(exact, estimated, param, nologo):
"""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 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))
plt.close()
def chi2(config):
"""Plot the χ² values of analysed variables.
"""
......@@ -348,7 +398,7 @@ def sed(config, sed_type, nologo):
"""Plot the best SED with associated observed and modelled fluxes.
"""
obs = read_table(config.configuration['data_file'])
mod = Table.read(OUT_DIR + BEST_MODEL_FILE)
mod = Table.read(OUT_DIR + BEST_RESULTS)
with Database() as base:
filters = OrderedDict([(name, base.get_filter(name))
......@@ -362,6 +412,38 @@ def sed(config, sed_type, nologo):
pool.join()
def mock(config, nologo):
"""Plot the comparison of input/output values of analysed variables.
"""
try:
exact = Table.read(OUT_DIR + BEST_RESULTS)
except FileNotFoundError:
print("Best models file {} not found.".format(OUT_DIR + BEST_RESULTS))
sys.exit(1)
try:
estimated = Table.read(OUT_DIR + MOCK_RESULTS)
except FileNotFoundError:
print("Mock models file {} not found.".format(OUT_DIR + MOCK_RESULTS))
sys.exit(1)
params = config.configuration['analysis_method_params']['analysed_variables']
for param in params:
if param.endswith('_log'):
param = "best."+param
exact[param] = np.log10(exact[param[:-4]])
arguments = ((exact["best."+param], estimated["bayes."+param], param, nologo)
for param in params)
with mp.Pool(processes=config.configuration['cores']) as pool:
pool.starmap(_mock_worker, arguments)
pool.close()
pool.join()
def main():
if sys.version_info[:2] >= (3, 4):
......@@ -388,6 +470,10 @@ def main():
sed_parser.add_argument('--nologo', action="store_true")
sed_parser.set_defaults(parser='sed')
sed_parser = subparsers.add_parser('mock', help=mock.__doc__)
sed_parser.add_argument('--nologo', action="store_true")
sed_parser.set_defaults(parser='mock')
args = parser.parse_args()
if args.config_file:
......@@ -395,9 +481,14 @@ def main():
else:
config = Configuration()
if args.parser == 'chi2':
chi2(config)
elif args.parser == 'pdf':
pdf(config)
elif args.parser == 'sed':
sed(config, args.type, args.nologo)
if len(sys.argv) == 1:
parser.print_usage()
else:
if args.parser == 'chi2':
chi2(config)
elif args.parser == 'pdf':
pdf(config)
elif args.parser == 'sed':
sed(config, args.type, args.nologo)
elif args.parser == 'mock':
mock(config, args.nologo)
......@@ -20,8 +20,7 @@ class custom_build(build):
entry_points = {
'console_scripts': ['pcigale = pcigale:main',
'pcigale-plots = pcigale_plots:main',
'pcigale-filters = pcigale_filters:main',
'pcigale-mock = pcigale_mock:main']
'pcigale-filters = pcigale_filters:main']
}
setup(
......@@ -36,8 +35,7 @@ setup(
cmdclass={"build": custom_build},
package_data={'pcigale': ['data/data.db'],
'pcigale_plots': ['data/CIGALE.png'],
'pcigale_mock': ['data/CIGALE.png']},
'pcigale_plots': ['data/CIGALE.png']},
author="The CIGALE team",
author_email="cigale@lam.fr",
......
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