Commit 58e65f0e authored by Médéric Boquien's avatar Médéric Boquien

Merge branch 'release/v0.9.0'

parents b8a59999 c502664a
This diff is collapsed.
......@@ -19,9 +19,10 @@ import itertools
import numpy as np
from scipy import interpolate
import scipy.constants as cst
from astropy.table import Table
from pcigale.data import (Database, Filter, M2005, BC03, Fritz2006,
Dale2014, DL2007, DL2014, NebularLines,
NebularContinuum)
NebularContinuum, Schreiber2016)
def read_bc03_ssp(filename):
......@@ -664,6 +665,34 @@ def build_nebular(base):
base.add_nebular_continuum(models_cont)
base.add_nebular_lines(models_lines)
def build_schreiber2016(base):
models = []
schreiber2016_dir = os.path.join(os.path.dirname(__file__),
'schreiber2016/')
print("Importing {}...".format(schreiber2016_dir + 'g15_pah.fits'))
pah = Table.read(schreiber2016_dir + 'g15_pah.fits')
print("Importing {}...".format(schreiber2016_dir + 'g15_dust.fits'))
dust = Table.read(schreiber2016_dir + 'g15_dust.fits')
# Getting the lambda grid for the templates and convert from μm to nm.
wave = dust['LAM'][0, 0, :].data * 1e3
for td in np.arange(15., 100.):
# Find the closest temperature in the model list of tdust
tsed = np.argmin(np.absolute(dust['TDUST'][0].data-td))
# The models are in νFν. We convert this to W/nm.
lumin_dust = dust['SED'][0, tsed, :].data / wave
lumin_pah = pah['SED'][0, tsed, :].data / wave
models.append(Schreiber2016(0, td, wave, lumin_dust))
models.append(Schreiber2016(1, td, wave, lumin_pah))
base.add_schreiber2016(models)
def build_base():
base = Database(writable=True)
base.upgrade_base()
......@@ -709,6 +738,11 @@ def build_base():
print("\nDONE\n")
print('#' * 78)
print("9- Importing Schreiber et al (2016) models\n")
build_schreiber2016(base)
print("\nDONE\n")
print('#' * 78)
base.session.close_all()
......
......@@ -2,10 +2,10 @@
# energy
# H beta pseudo filter
4827.87 0.0
4827.875 -0.2857959417043478
4847.87 -0.2857959417043478
4847.875 0.34782608695652173
4876.625 0.34782608695652173
4876.63 -0.2857959417043478
4891.625 -0.2857959417043478
4827.875 0.2857959417043478
4847.87 0.2857959417043478
4847.875 -0.34782608695652173
4876.625 -0.34782608695652173
4876.63 0.2857959417043478
4891.625 0.2857959417043478
4891.63 0.0
......@@ -2,14 +2,14 @@
# energy
# H delta pseudo filter
4041.5 0.0
4041.6000000000004 -0.1415428166967742
4079.75 -0.1415428166967742
4041.6000000000004 0.1415428166967742
4079.75 0.1415428166967742
4079.8 0.0
4083.3999999999996 0.0
4083.5 0.25806451612903225
4122.25 0.25806451612903225
4083.5 -0.25806451612903225
4122.25 -0.25806451612903225
4122.3 0.0
4128.4 0.0
4128.5 -0.1415428166967742
4161.0 -0.1415428166967742
4128.5 0.1415428166967742
4161.0 0.1415428166967742
4161.1 0.0
......@@ -2,12 +2,12 @@
# energy
# H gamma pseudo filter
4283.4 0.0
4283.5 -0.11268875366857142
4319.74 -0.11268875366857142
4319.75 0.22857142857142856
4363.5 0.22857142857142856
4283.5 0.11268875366857142
4319.74 0.11268875366857142
4319.75 -0.22857142857142856
4363.5 -0.22857142857142856
4363.6 0.0
4367.200000000001 0.0
4367.25 -0.11268875366857142
4419.75 -0.11268875366857142
4367.25 0.11268875366857142
4419.75 0.11268875366857142
4419.8 0.0
......@@ -2,14 +2,14 @@
# energy
# Mg2 pseudo filter
4895.12 0.0
4895.125 -0.07843137254117646
4957.625 -0.07843137254117646
4895.125 0.07843137254117646
4957.625 0.07843137254117646
4957.63 0.0
5154.120000000001 0.0
5154.125 0.23529411764705882
5196.625 0.23529411764705882
5154.125 -0.23529411764705882
5196.625 -0.23529411764705882
5196.63 0.0
5301.12 0.0
5301.125 -0.07843137254117646
5366.125 -0.07843137254117646
5301.125 0.07843137254117646
5366.125 0.07843137254117646
5366.130000000001 0.0
......@@ -8,8 +8,8 @@ import multiprocessing as mp
import sys
from .session.configuration import Configuration
from .analysis_modules import get_module as get_analysis_module
from .analysis_modules.utils import ParametersHandler
from .analysis_modules import get_module
from .handlers.parameters_handler import ParametersHandler
__version__ = "0.1-alpha"
......@@ -36,28 +36,21 @@ def check(config):
"""
# TODO: Check if all the parameters that don't have default values are
# given for each module.
print("With this configuration, pcigale must compute {} "
"SEDs.".format(ParametersHandler(
config.configuration['creation_modules'],
config.configuration['creation_modules_params']
).size))
configuration = config.configuration
if configuration:
print("With this configuration cigale will compute {} "
"models.".format(ParametersHandler(configuration).size))
def run(config):
"""Run the analysis.
"""
data_file = config.configuration['data_file']
column_list = config.configuration['column_list']
creation_modules = config.configuration['creation_modules']
creation_modules_params = config.configuration['creation_modules_params']
analysis_module = get_analysis_module(config.configuration[
'analysis_method'])
analysis_module_params = config.configuration['analysis_method_params']
cores = config.configuration['cores']
analysis_module.process(data_file, column_list, creation_modules,
creation_modules_params, analysis_module_params,
cores)
configuration = config.configuration
if configuration:
analysis_module = get_module(configuration['analysis_method'])
analysis_module.process(configuration)
def main():
......
......@@ -31,8 +31,7 @@ class AnalysisModule(object):
# module parameter.
self.parameters = kwargs
def _process(self, data_file, column_list, creation_modules,
creation_modules_params, parameters):
def _process(self, configuration):
"""Do the actual analysis
This method is responsible for the fitting / analysis process
......@@ -40,19 +39,8 @@ class AnalysisModule(object):
Parameters
----------
data_file: string
Name of the file containing the observations to be fitted.
column_list: array of strings
Names of the columns from the data file to use in the analysis.
creation_modules: array of strings
Names (in the right order) of the modules to use to build the SED.
creation_modules_params: array of array of dictionaries
Array containing all the possible combinations of configurations
for the creation_modules. Each 'inner' array has the same length as
the creation_modules array and contains the configuration
dictionary for the corresponding module.
parameters: dictionary
Configuration for the module.
configuration: dictionary
Configuration file
Returns
-------
......@@ -61,8 +49,7 @@ class AnalysisModule(object):
"""
raise NotImplementedError()
def process(self, data_file, column_list, creation_modules,
creation_modules_params, parameters):
def process(self, configuration):
"""Process with the analysis
This method is responsible for checking the module parameters before
......@@ -72,19 +59,8 @@ class AnalysisModule(object):
Parameters
----------
data_file: string
Name of the file containing the observations to be fitted.
column_list: array of strings
Names of the columns from the data file to use in the analysis.
creation_modules: array of strings
Names (in the right order) of the modules to use to build the SED.
creation_modules_params: array of array of dictionaries
Array containing all the possible combinations of configurations
for the creation_modules. Each 'inner' array has the same length as
the creation_modules array and contains the configuration
dictionary for the corresponding module.
parameters: dictionary
Configuration for the module.
configuration: dictionary
Contents of pcigale.ini in the form of a dictionary
Returns
-------
......@@ -95,6 +71,7 @@ class AnalysisModule(object):
KeyError: when not all the needed parameters are given.
"""
parameters = configuration['analysis_params']
# For parameters that are present on the parameter_list with a default
# value and that are not in the parameters dictionary, we add them
# with their default value.
......@@ -124,8 +101,7 @@ class AnalysisModule(object):
"expected one." + message)
# We do the actual processing
self._process(data_file, column_list, creation_modules,
creation_modules_params, parameters)
self._process(configuration)
def get_module(module_name):
......
......@@ -29,10 +29,10 @@ def save_best_sed(obsid, sed, norm):
Normalisation factor to scale the scale to the observations
"""
sed.to_votable(OUT_DIR + "{}_best_model.xml".format(obsid), mass=norm)
sed.to_fits(OUT_DIR + "{}".format(obsid), mass=norm)
def save_pdf(obsid, name, model_variable, likelihood):
def _save_pdf(obsid, name, model_variable, likelihood):
"""Compute and save the PDF to a FITS file
We estimate the probability density functions (PDF) of the parameter from
......@@ -83,7 +83,40 @@ def save_pdf(obsid, name, model_variable, likelihood):
table.write(OUT_DIR + "{}_{}_pdf.fits".format(obsid, name))
def save_chi2(obsid, name, model_variable, chi2):
def save_pdf(obsid, names, mass_proportional, model_variables, scaling,
likelihood, wlikely):
"""Compute and save the PDF of analysed variables
Parameters
----------
obsid: string
Name of the object. Used to prepend the output file name
names: list of strings
Analysed variables names
model_variables: array
Values of the model variables
likelihood: 1D array
Likelihood of the "likely" models
"""
for i, name in enumerate(names):
if name.endswith('_log'):
if name[:-4] in mass_proportional:
model_variable = np.log10(model_variables[:, i][wlikely] *
scaling[wlikely])
else:
model_variable = np.log10(model_variables[:, i][wlikely])
else:
if name in mass_proportional:
model_variable = (model_variables[:, i][wlikely] *
scaling[wlikely])
else:
model_variable = model_variables[:, i][wlikely]
_save_pdf(obsid, name, model_variable, likelihood)
def _save_chi2(obsid, name, model_variable, chi2):
"""Save the best reduced χ² versus an analysed variable
Parameters
......@@ -103,94 +136,107 @@ def save_chi2(obsid, name, model_variable, chi2):
table.write(OUT_DIR + "{}_{}_chi2.fits".format(obsid, name))
def save_table_analysis(filename, obsid, analysed_variables, analysed_averages,
analysed_std):
"""Save the estimated values derived from the analysis of the PDF
def save_chi2(obsid, names, mass_proportional, model_variables, scaling, chi2):
"""Save the best reduced χ² versus analysed variables
Parameters
----------
filename: name of the file to save
Name of the output file
obsid: table column
Names of the objects
analysed_variables: list
Analysed variable names
analysed_averages: RawArray
Analysed variables values estimates
analysed_std: RawArray
Analysed variables errors estimates
obsid: string
Name of the object. Used to prepend the output file name
name: list of strings
Analysed variables names
model_variables: array
Values of the model variables
scaling: array
Scaling factors of the models
chi2:
Reduced χ²
"""
np_analysed_averages = np.ctypeslib.as_array(analysed_averages[0])
np_analysed_averages = np_analysed_averages.reshape(analysed_averages[1])
np_analysed_std = np.ctypeslib.as_array(analysed_std[0])
np_analysed_std = np_analysed_std.reshape(analysed_std[1])
result_table = Table()
result_table.add_column(Column(obsid.data, name="observation_id"))
for index, variable in enumerate(analysed_variables):
result_table.add_column(Column(
np_analysed_averages[:, index],
name=variable
))
result_table.add_column(Column(
np_analysed_std[:, index],
name=variable+"_err"
))
result_table.write(OUT_DIR + filename, format='ascii.fixed_width',
delimiter=None)
def save_table_best(filename, obsid, chi2, chi2_red, variables, fluxes,
filters, info_keys):
"""Save the values corresponding to the best fit
for i, name in enumerate(names):
if name.endswith('_log'):
if name[:-4] in mass_proportional:
model_variable = np.log10(model_variables[:, i] * scaling)
else:
model_variable = np.log10(model_variables[:, i])
else:
if name in mass_proportional:
model_variable = model_variables[:, i] * scaling
else:
model_variable = model_variables[:, i]
_save_chi2(obsid, name, model_variable, chi2)
def save_results(filename, obsid, bayes_variables, bayes_mean, bayes_std, chi2,
chi2_red, best, fluxes, filters, info_keys):
"""Save the estimated values derived from the analysis of the PDF and the
parameters associated with the best fit. An simple text file and a FITS
file are generated.
Parameters
----------
filename: name of the file to save
Name of the output file
filename: string
Name of the output file without the extension
obsid: table column
Names of the objects
bayes_variables: list
Analysed variable names
bayes_mean: RawArray
Analysed variables values estimates
bayes_std: RawArray
Analysed variables errors estimates
chi2: RawArray
Best χ² for each object
chi2_red: RawArray
Best reduced χ² for each object
variables: RawArray
All variables corresponding to a SED
best: RawArray
All variables corresponding to the best SED
fluxes: RawArray
Fluxes in all bands for each object
Fluxes in all bands for the best SED
filters: list
Filters used to compute the fluxes
info_keys: list
Parameters names
"""
np_bayes_mean = np.ctypeslib.as_array(bayes_mean[0])
np_bayes_mean = np_bayes_mean.reshape(bayes_mean[1])
np_bayes_std = np.ctypeslib.as_array(bayes_std[0])
np_bayes_std = np_bayes_std.reshape(bayes_std[1])
np_fluxes = np.ctypeslib.as_array(fluxes[0])
np_fluxes = np_fluxes.reshape(fluxes[1])
np_variables = np.ctypeslib.as_array(variables[0])
np_variables = np_variables.reshape(variables[1])
np_best = np.ctypeslib.as_array(best[0])
np_best = np_best.reshape(best[1])
np_chi2 = np.ctypeslib.as_array(chi2[0])
np_chi2_red = np.ctypeslib.as_array(chi2_red[0])
best_model_table = Table()
best_model_table.add_column(Column(obsid.data, name="observation_id"))
best_model_table.add_column(Column(np_chi2, name="chi_square"))
best_model_table.add_column(Column(np_chi2_red, name="reduced_chi_square"))
table = Table()
table.add_column(Column(obsid.data, name="id"))
for idx, name in enumerate(bayes_variables):
table.add_column(Column(np_bayes_mean[:, idx], name="bayes."+name))
table.add_column(Column(np_bayes_std[:, idx],
name="bayes."+name+"_err"))
for index, name in enumerate(info_keys):
column = Column(np_variables[:, index], name=name)
best_model_table.add_column(column)
table.add_column(Column(np_chi2, name="best.chi_square"))
table.add_column(Column(np_chi2_red, name="best.reduced_chi_square"))
for index, name in enumerate(filters):
column = Column(np_fluxes[:, index], name=name, unit='mJy')
best_model_table.add_column(column)
for idx, name in enumerate(info_keys):
table.add_column(Column(np_best[:, idx], name="best."+name))
best_model_table.write(OUT_DIR + filename, format='ascii.fixed_width',
delimiter=None)
for idx, name in enumerate(filters):
table.add_column(Column(np_fluxes[:, idx], name="best."+name,
unit='mJy'))
table.write(OUT_DIR+filename+".txt", format='ascii.fixed_width',
delimiter=None)
table.write(OUT_DIR+filename+".fits", format='fits')
def dchi2_over_ds2(s, obs_fluxes, obs_errors, mod_fluxes):
......@@ -335,7 +381,16 @@ def compute_chi2(model_fluxes, obs_fluxes, obs_errors, lim_flag):
scaling: array
scaling of the models to obtain the minimum χ²
"""
limits = lim_flag and np.any(obs_errors <= 0.)
scaling = _compute_scaling(model_fluxes, obs_fluxes, obs_errors)
# Some observations may not have flux values in some filter(s), but
# they can have upper limit(s).
if limits == True:
for imod in range(len(model_fluxes)):
scaling[imod] = optimize.root(dchi2_over_ds2, scaling[imod],
args=(obs_fluxes, obs_errors,
model_fluxes[imod, :])).x
# χ² of the comparison of each model to each observation.
chi2 = np.zeros(model_fluxes.shape[0])
......@@ -346,11 +401,7 @@ def compute_chi2(model_fluxes, obs_fluxes, obs_errors, lim_flag):
# Some observations may not have flux values in some filter(s), but
# they can have upper limit(s).
if (lim_flag and np.any(obs_errors <= 0.)) == True:
for imod in range(len(model_fluxes)):
scaling[imod] = optimize.root(dchi2_over_ds2, scaling[imod],
args=(obs_fluxes, obs_errors,
model_fluxes[imod, :])).x
if limits == True:
mask_lim = (obs_errors <= 0.)
chi2 += -2. * np.sum(
np.log(
......
......@@ -8,6 +8,7 @@
import time
from astropy.cosmology import WMAP7 as cosmology
import numpy as np
import scipy.stats as st
......@@ -16,8 +17,7 @@ from .utils import (save_best_sed, save_pdf, save_chi2, compute_chi2,
from ...warehouse import SedWarehouse
def init_sed(params, filters, analysed, redshifts, fluxes, variables,
t_begin, n_computed):
def init_sed(params, filters, analysed, fluxes, variables, t_begin, n_computed):
"""Initializer of the pool of processes. It is mostly used to convert
RawArrays into numpy arrays. The latter are defined as global variables to
be accessible from the workers.
......@@ -30,8 +30,6 @@ def init_sed(params, filters, analysed, redshifts, fluxes, variables,
Contains the names of the filters to compute the fluxes.
analysed: list
Variable names to be analysed.
redshifts: RawArray and tuple containing the shape
Redshifts of individual models. Shared among workers.
fluxes: RawArray and tuple containing the shape
Fluxes of individual models. Shared among workers.
variables: RawArray and tuple containing the shape
......@@ -42,11 +40,9 @@ def init_sed(params, filters, analysed, redshifts, fluxes, variables,
Time of the beginning of the computation.
"""
global gbl_model_redshifts, gbl_model_fluxes, gbl_model_variables
global gbl_n_computed, gbl_t_begin, gbl_params, gbl_previous_idx
global gbl_filters, gbl_analysed_variables, gbl_warehouse
gbl_model_redshifts = np.ctypeslib.as_array(redshifts[0])
global gbl_model_fluxes, gbl_model_variables, gbl_n_computed, gbl_t_begin
global gbl_params, gbl_previous_idx, gbl_filters, gbl_analysed_variables
global gbl_warehouse
gbl_model_fluxes = np.ctypeslib.as_array(fluxes[0])
gbl_model_fluxes = gbl_model_fluxes.reshape(fluxes[1])
......@@ -67,7 +63,7 @@ def init_sed(params, filters, analysed, redshifts, fluxes, variables,
gbl_warehouse = SedWarehouse()
def init_analysis(params, filters, analysed, redshifts, fluxes, variables,
def init_analysis(params, filters, analysed, z, fluxes, variables,
t_begin, n_computed, analysed_averages, analysed_std,
best_fluxes, best_parameters, best_chi2, best_chi2_red, save,
lim_flag, n_obs):
......@@ -83,7 +79,7 @@ def init_analysis(params, filters, analysed, redshifts, fluxes, variables,
Contains filters to compute the fluxes.
analysed: list
Variable names to be analysed
redshifts: RawArray and tuple containing the shape.
z: RawArray and tuple containing the shape.
Redshifts of individual models. Shared among workers.
fluxes: RawArray
Fluxes of individual models. Shared among workers.
......@@ -112,9 +108,8 @@ def init_analysis(params, filters, analysed, redshifts, fluxes, variables,
Number of observations.
"""
init_sed(params, filters, analysed, redshifts, fluxes, variables,
t_begin, n_computed)
global gbl_redshifts, gbl_analysed_averages, gbl_analysed_std
init_sed(params, filters, analysed, fluxes, variables, t_begin, n_computed)
global gbl_z, gbl_analysed_averages, gbl_analysed_std
global gbl_best_fluxes, gbl_best_parameters, gbl_best_chi2
global gbl_best_chi2_red, gbl_save, gbl_n_obs, gbl_lim_flag, gbl_keys
......@@ -134,8 +129,7 @@ def init_analysis(params, filters, analysed, redshifts, fluxes, variables,
gbl_best_chi2_red = np.ctypeslib.as_array(best_chi2_red[0])
gbl_redshifts = gbl_model_redshifts[np.unique(gbl_model_redshifts,
return_index=True)[1]]
gbl_z = z
gbl_save = save
gbl_lim_flag = lim_flag
......@@ -174,8 +168,6 @@ def sed(idx):
for name in
gbl_analysed_variables])
gbl_model_redshifts[idx] = sed.info['universe.redshift']
with gbl_n_computed.get_lock():
gbl_n_computed.value += 1
n_computed = gbl_n_computed.value
......@@ -204,15 +196,34 @@ def analysis(idx, obs):
obs_fluxes = np.array([obs[name] for name in gbl_filters])
obs_errors = np.array([obs[name + "_err"] for name in gbl_filters])
obs_z = obs['redshift']
nobs = np.where(np.isfinite(obs_fluxes))[0].size
if obs['redshift'] >= 0.:
if obs_z >= 0.:
# We pick the the models with the closest redshift using a slice to
# work on views of the arrays and not on copies to save on RAM.
wz = slice(np.abs(obs['redshift'] - gbl_redshifts).argmin(), None,
gbl_redshifts.size)
idx_z = np.abs(obs_z - gbl_z).argmin()
model_z = gbl_z[idx_z]
wz = slice(idx_z, None, gbl_z.size)
# The mass-dependent physical properties are computed assuming the
# redshift of the model. However because we round the observed redshifts
# to two decimals, there can be a difference of 0.005 in redshift
# between the models and the actual observation. At low redshift, this
# can cause a discrepancy in the mass-dependent physical properties:
# ~0.35 dex at z=0.010 vs 0.015 for instance. Therefore we correct these
# physical quantities by multiplying them by corr_dz.
if model_z == obs_z:
corr_dz = 1.
else:
if model_z > 0.:
corr_dz = (cosmology.luminosity_distance(obs_z).value /
cosmology.luminosity_distance(model_z).value)**2.
else:
corr_dz = (cosmology.luminosity_distance(obs_z).value * 1e5)**2.
else: # We do not know the redshift so we use the full grid
wz = slice(0, None, 1)
corr_dz = 1.
chi2, scaling = compute_chi2(gbl_model_fluxes[wz, :], obs_fluxes,
obs_errors, gbl_lim_flag)
......@@ -260,13 +271,13 @@ def analysis(idx, obs):
# We correct the mass-dependent parameters
for key in sed.mass_proportional_info:
sed.info[key] *= scaling[best_index_z]
sed.info[key] *= scaling[best_index_z] * corr_dz
# We compute the weighted average and standard deviation using the
# likelihood as weight.
for i, variable in enumerate(gbl_analysed_variables):
if variable.endswith('_log'):
variable = variable.strip('_log')
variable = variable[:-4]
_ = np.log10
maxstd = lambda mean, std: max(0.02, std)
else:
......@@ -275,7 +286,8 @@ def analysis(idx, obs):
if variable in sed.mass_proportional_info:
mean, std = weighted_param(_(gbl_model_variables[wz, i][wlikely]
* scaling[wlikely]), likelihood)
* scaling[wlikely] * corr_dz),
likelihood)
else:
mean, std = weighted_param(_(gbl_model_variables[wz, i][wlikely]),
likelihood)
......@@ -297,22 +309,13 @@ def analysis(idx, obs):
if gbl_save['best_sed']:
save_best_sed(obs['id'], sed, scaling[best_index_z])
if gbl_save['chi2']:
for i, variable in enumerate(gbl_analysed_variables):
if variable in sed.mass_proportional_info:
save_chi2(obs['id'], variable, gbl_model_variables[wz, i] *
scaling, chi2 / (nobs - 1))
else:
save_chi2(obs['id'], variable, gbl_model_variables[wz, i],
chi2 / (nobs - 1))
save_chi2(obs['id'], gbl_analysed_variables,
sed.mass_proportional_info, gbl_model_variables[wz, :],
scaling * corr_dz, chi2 / (nobs - 1))
if gbl_save['pdf']:
for i, variable in enumerate(gbl_analysed_variables):
if variable in sed.mass_proportional_info:
save_pdf(obs['id'], variable,
gbl_model_variables[wz, i][wlikely] *
scaling[wlikely], likelihood)
else:
save_pdf(obs['id'], variable,
gbl_model_variables[wz, i][wlikely], likelihood)
save_pdf(obs['id'], gbl_analysed_variables,
sed.mass_proportional_info, gbl_model_variables[wz, :],
scaling * corr_dz, likelihood, wlikely)