...
 
Commits (28)
# Change Log
## Unreleased
### Added
### Changed
### Fixed
- Make sure we can plot the PDF of equivalent widths. (Médéric Boquien)
- Fix a crash when generating a mock catalogue containing intensive properties. (Médéric Boquien)
- In the `sfhdelayed` and `sfhdelayedbq` modules, provide the correct description for the sfr_A parameter (Médéric Boquien & Laure Ciesla)
### Optimised
## 2018.0 (2018-11-06)
### Added
- It is now possible to optionally indicate the distance in Mpc in the input file. If present it will be used in lieu of the distance computed from the redshift. This is especially useful in the nearby universe where the redshift is a very poor indicator of the actual distance. (Médéric Boquien)
......
......@@ -66,6 +66,6 @@ class Counter:
def pprint(self, n):
dt = time.time() - self.t0
print("{}/{} computed in {:.1f} seconds ({:.1f}/s)".
print("{}/{} performed in {:.1f} seconds ({:.1f}/s)".
format(n, self.nmodels, dt, n / dt),
end="\n" if n == self.nmodels else "\r")
......@@ -235,7 +235,7 @@ class ObservationsManagerPassbands(object):
np.fabs(self.table[err]))
for prop in self.intprops:
err = prop + '_err'
self.table[name] = np.random.normal(fits.best.intprop[prop],
self.table[prop] = np.random.normal(fits.best.intprop[prop],
np.fabs(self.table[err]))
for prop in self.extprops:
err = prop + '_err'
......
......@@ -63,7 +63,8 @@ class SFHDelayed(SedModule):
)),
("sfr_A", (
"cigale_list(minvalue=0.)",
"Value of SFR at t = 0 in M_sun/yr.",
"Multiplicative factor controlling the SFR if normalise is False. "
"For instance without any burst: SFR(t)=sfr_A×t×exp(-t/τ)/τ²",
1.
)),
("normalise", (
......
......@@ -59,7 +59,8 @@ class SFHDelayedBQ(SedModule):
)),
("sfr_A", (
"cigale_list(minvalue=0.)",
"Value of SFR at t = 0 in M_sun/yr.",
"Multiplicative factor controlling the SFR if normalise is False. "
"For instance without any burst/quench: SFR(t)=sfr_A×t×exp(-t/τ)/τ²",
1.
)),
("normalise", (
......
......@@ -100,7 +100,7 @@ class Configuration(object):
["* sfhdelayed (delayed SFH with optional exponential burst)"] +
["* sfhdelayedbq (delayed SFH with optional constant burst/quench)"
] +
["* sfhfromfile (abitrary SFH read from an input file)"] +
["* sfhfromfile (arbitrary SFH read from an input file)"] +
["* sfhperiodic (periodic SFH, exponential, rectangle or delayed"
")"] +
["SSP:"] +
......
This diff is collapsed.
# -*- coding: utf-8 -*-
# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
# Copyright (C) 2013-2014 Yannick Roehlly
# Copyright (C) 2013 Institute of Astronomy
# Copyright (C) 2014 Laboratoire d'Astrophysique de Marseille
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Yannick Roehlly, Médéric Boquien & Denis Burgarella
import glob
from itertools import product
from os import path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import multiprocessing as mp
import numpy as np
from pcigale.utils import read_table
from pcigale.analysis_modules.utils import Counter, nothread
def pool_initializer(counter):
"""Initializer of the pool of processes to share variables between workers.
Parameters
----------
:param counter: Counter class object for the number of models plotted
"""
global gbl_counter
# Limit the number of threads to 1 if we use MKL in order to limit the
# oversubscription of the CPU/RAM.
nothread()
gbl_counter = counter
def chi2(config, format, outdir):
"""Plot the χ² values of analysed variables.
"""
file = path.join(path.dirname(outdir), config.configuration['data_file'])
input_data = read_table(file)
chi2_vars = config.configuration['analysis_params']['variables']
chi2_vars += [band for band in config.configuration['bands']
if band.endswith('_err') is False]
items = list(product(input_data['id'], chi2_vars, [format], [outdir]))
counter = Counter(len(items))
with mp.Pool(processes=config.configuration['cores'], initializer=pool_initializer,
initargs=(counter,)) as pool:
pool.starmap(_chi2_worker, items)
pool.close()
pool.join()
def _chi2_worker(obj_name, var_name, format, outdir):
"""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..
outdir: string
The absolute path to outdir
"""
gbl_counter.inc()
figure = plt.figure()
ax = figure.add_subplot(111)
var_name = var_name.replace('/', '_')
fnames = glob.glob("{}/{}_{}_chi2-block-*.npy".format(outdir, 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("{}/{}_{}_chi2.{}".format(outdir, obj_name, var_name, format))
plt.close(figure)
# -*- coding: utf-8 -*-
# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
# Copyright (C) 2013-2014 Yannick Roehlly
# Copyright (C) 2013 Institute of Astronomy
# Copyright (C) 2014 Laboratoire d'Astrophysique de Marseille
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Yannick Roehlly, Médéric Boquien & Denis Burgarella
from astropy.table import Table
import matplotlib
import sys
from os import path
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import multiprocessing as mp
import numpy as np
import pkg_resources
from scipy import stats
from pcigale.analysis_modules.utils import Counter, nothread
# Name of the file containing the best models information
BEST_RESULTS = "results.fits"
MOCK_RESULTS = "results_mock.fits"
def pool_initializer(counter):
"""Initializer of the pool of processes to share variables between workers.
Parameters
----------
:param counter: Counter class object for the number of models plotted
"""
global gbl_counter
# Limit the number of threads to 1 if we use MKL in order to limit the
# oversubscription of the CPU/RAM.
nothread()
gbl_counter = counter
def mock(config, nologo, outdir):
"""Plot the comparison of input/output values of analysed variables.
"""
best_results_file = path.abspath(path.join(outdir, BEST_RESULTS))
mock_results_file = path.abspath(path.join(outdir, MOCK_RESULTS))
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, outdir)
for param in params]
counter = Counter(len(arguments))
with mp.Pool(processes=config.configuration['cores'], initializer=pool_initializer,
initargs=(counter,)) as pool:
pool.starmap(_mock_worker, arguments)
pool.close()
pool.join()
def _mock_worker(exact, estimated, param, logo, outdir):
"""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.
outdir: string
The absolute path to outdir
"""
gbl_counter.inc()
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('{}/mock_{}.pdf'.format(outdir, param))
plt.close()
# -*- coding: utf-8 -*-
# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
# Copyright (C) 2013-2014 Yannick Roehlly
# Copyright (C) 2013 Institute of Astronomy
# Copyright (C) 2014 Laboratoire d'Astrophysique de Marseille
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Yannick Roehlly, Médéric Boquien & Denis Burgarella
import glob
from itertools import product
import matplotlib
from os import path
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import multiprocessing as mp
import numpy as np
from pcigale.utils import read_table
from pcigale.analysis_modules.utils import Counter, nothread
def pool_initializer(counter):
"""Initializer of the pool of processes to share variables between workers.
Parameters
----------
:param counter: Counter class object for the number of models plotted
"""
global gbl_counter
# Limit the number of threads to 1 if we use MKL in order to limit the
# oversubscription of the CPU/RAM.
nothread()
gbl_counter = counter
def pdf(config, format, outdir):
"""Plot the PDF of analysed variables.
"""
input_data = read_table(path.join(path.dirname(outdir), 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]
items = list(product(input_data['id'], pdf_vars, [format], [outdir]))
counter = Counter(len(items))
with mp.Pool(processes=config.configuration['cores'], initializer=pool_initializer,
initargs=(counter,)) as pool:
pool.starmap(_pdf_worker, items)
pool.close()
pool.join()
def _pdf_worker(obj_name, var_name, format, outdir):
"""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..
outdir: string
The absolute path to outdir
"""
gbl_counter.inc()
var_name = var_name.replace('/', '_')
if var_name.endswith('_log'):
fnames = glob.glob("{}/{}_{}_chi2-block-*.npy".format(outdir, obj_name,
var_name[:-4]))
log = True
else:
fnames = glob.glob("{}/{}_{}_chi2-block-*.npy".format(outdir, 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("{}/{}_{}_pdf.{}".format(outdir, obj_name, var_name, format))
plt.close(figure)
This diff is collapsed.
......@@ -48,7 +48,7 @@ setup(
cmdclass={"build": custom_build},
package_data={'pcigale': ['data/data.db'],
'pcigale_plots': ['data/CIGALE.png']},
'pcigale_plots': ['resources/CIGALE.png']},
author="The CIGALE team",
author_email="cigale@lam.fr",
......