Commit 79f1425e authored by BURGARELLA Denis's avatar BURGARELLA Denis

create mock catalogue from best model and analyse it

parent 8db7d37d
......@@ -44,6 +44,8 @@ from .workers import init_analysis as init_worker_analysis
from .workers import analysis as worker_analysis
from ..utils import ParametersHandler, backup_dir
from ..utils import OUT_DIR
# Tolerance threshold under which any flux or error is considered as 0.
TOLERANCE = 1e-12
# Limit the redshift to this number of decimals
......@@ -82,6 +84,12 @@ class PdfAnalysis(AnalysisModule):
"If true, for each object check whether upper limits are present "
"and analyse them.",
False
)),
("mock_flag", (
"boolean",
"If true, for each object we create a mock object "
"and analyse them.",
False
))
])
......@@ -124,6 +132,7 @@ class PdfAnalysis(AnalysisModule):
save = {key: config["save_{}".format(key)].lower() == "true"
for key in ["best_sed", "chi2", "pdf"]}
lim_flag = config["lim_flag"].lower() == "true"
mock_flag = config["mock_flag"].lower() == "true"
# Get the needed filters in the pcigale database. We use an ordered
# dictionary because we need the keys to always be returned in the
......@@ -193,7 +202,7 @@ class PdfAnalysis(AnalysisModule):
initargs=initargs) as pool:
pool.map(worker_sed, range(n_params))
print('\nAnalysing models...')
print("\nAnalysing models...")
# We use RawArrays for the same reason as previously
analysed_averages = (RawArray(ctypes.c_double, n_obs * n_variables),
......@@ -221,12 +230,96 @@ class PdfAnalysis(AnalysisModule):
initargs=initargs) as pool:
pool.starmap(worker_analysis, enumerate(obs_table))
val_best_chi2_red = np.ctypeslib.as_array(best_chi2_red[0])
verylow_redchi2 = sum(chi2_red < TOLERANCE
for chi2_red in val_best_chi2_red)/n_obs*100.
low_redchi2 = sum(chi2_red < 0.5
for chi2_red in val_best_chi2_red)/n_obs*100.
# If low values of reduced chi^2, it means that the data are overfitted
# Errors might be under-estimated or not enough valid data.
if verylow_redchi2 > 10. or low_redchi2 > 20.:
print("\nWarning: "
"{} % of the objects have chi^2_red ~ 0. and {} % chi^2_red < 0.5"
.format(np.round(verylow_redchi2,1), np.round(low_redchi2,1)))
save_table_analysis('analysis_results.txt', obs_table['id'],
analysed_variables, analysed_averages, analysed_std)
save_table_best('best_models.txt', obs_table['id'], best_chi2,
best_chi2_red, best_parameters, best_fluxes, filters, info)
if mock_flag is True:
print("\nMock analysis...")
# We use RawArrays for the same reason as previously
analysed_averages_mock = (RawArray(ctypes.c_double, n_obs * n_variables),
(n_obs, n_variables))
analysed_std_mock = (RawArray(ctypes.c_double, n_obs * n_variables),
(n_obs, n_variables))
best_fluxes_mock = (RawArray(ctypes.c_double, n_obs * n_filters),
(n_obs, n_filters))
best_parameters_mock = (RawArray(ctypes.c_double, n_obs * n_info),
(n_obs, n_info))
best_chi2_mock = (RawArray(ctypes.c_double, n_obs), (n_obs))
best_chi2_red_mock = (RawArray(ctypes.c_double, n_obs), (n_obs))
obs_fluxes = np.array([obs_table[name] for name in filters]).T
obs_errors = np.array([obs_table[name + "_err"] for name in filters]).T
mock_fluxes = np.empty_like(obs_fluxes)
mock_errors = np.empty_like(obs_errors)
bestmod_fluxes = np.ctypeslib.as_array(best_fluxes[0])
bestmod_fluxes = bestmod_fluxes.reshape(best_fluxes[1])
wdata = np.where((obs_fluxes > 0.)&(obs_errors > 0.))
wlim = np.where((obs_fluxes > 0.)&
(obs_errors >= -9990.)&(obs_errors < 0.))
wnodata = np.where((obs_fluxes <= -9990.)&(obs_errors <= -9990.))
mock_fluxes[wdata] = np.random.normal(
bestmod_fluxes[wdata],
obs_errors[wdata],
len(bestmod_fluxes[wdata]))
mock_errors[wdata] = obs_errors[wdata]
mock_fluxes[wlim] = obs_fluxes[wlim]
mock_errors[wlim] = obs_errors[wlim]
mock_fluxes[wnodata] = obs_fluxes[wnodata]
mock_errors[wnodata] = obs_errors[wnodata]
mock_table = obs_table.copy()
mock_table['id'] = obs_table['id']
mock_table['redshift'] = obs_table['redshift']
indx = 0
for name in filters:
mock_table[name] = mock_fluxes[:, indx]
mock_table[name + "_err"] = mock_errors[:, indx]
indx += 1
initargs = (params, filters, analysed_variables, model_redshifts,
model_fluxes, model_variables, time.time(),
mp.Value('i', 0), analysed_averages_mock, analysed_std_mock,
best_fluxes_mock, best_parameters_mock, best_chi2_mock,
best_chi2_red_mock, save, lim_flag, n_obs)
if cores == 1: # Do not create a new process
init_worker_analysis(*initargs)
for idx, mock in enumerate(mock_table):
worker_analysis(idx, mock)
else: # Analyse observations in parallel
with mp.Pool(processes=cores, initializer=init_worker_analysis,
initargs=initargs) as pool:
pool.starmap(worker_analysis, enumerate(mock_table))
print("\nSaving results...")
save_table_analysis(obs_table['id'], analysed_variables,
analysed_averages, analysed_std)
save_table_best(obs_table['id'], best_chi2, best_chi2_red,
best_parameters, best_fluxes, filters, info)
if mock_flag is True:
save_table_analysis('analysis_mock_results.txt', mock_table['id'],
analysed_variables, analysed_averages_mock, analysed_std_mock)
save_table_best('best_mock_models.txt', mock_table['id'],
best_chi2_mock, best_chi2_red_mock,
best_parameters_mock, best_fluxes_mock, filters, info)
print("Run completed!")
......
......@@ -15,11 +15,6 @@ from ..utils import OUT_DIR
# Number of points in the PDF
PDF_NB_POINTS = 1000
# Name of the file containing the analysis results
RESULT_FILE = "analysis_results.txt"
# Name of the file containing the best models information
BEST_MODEL_FILE = "best_models.txt"
def save_best_sed(obsid, sed, norm):
"""Save the best SED to a VO table.
......@@ -95,12 +90,14 @@ def save_chi2(obsid, analysed_variables, model_variables, reduced_chi2):
table.write(OUT_DIR + "{}_{}_chi2.fits".format(obsid, var_name))
def save_table_analysis(obsid, analysed_variables, analysed_averages,
def save_table_analysis(filename, obsid, analysed_variables, analysed_averages,
analysed_std):
"""Save the estimated values derived from the analysis of the PDF
Parameters
----------
filename: name of the file to save
Name of the output file
obsid: table column
Names of the objects
analysed_variables: list
......@@ -128,15 +125,17 @@ def save_table_analysis(obsid, analysed_variables, analysed_averages,
np_analysed_std[:, index],
name=variable+"_err"
))
result_table.write(OUT_DIR + RESULT_FILE, format='ascii.commented_header')
result_table.write(OUT_DIR + filename, format='ascii.commented_header')
def save_table_best(obsid, chi2, chi2_red, variables, fluxes, filters,
def save_table_best(filename, obsid, chi2, chi2_red, variables, fluxes, filters,
info_keys):
"""Save the values corresponding to the best fit
Parameters
----------
filename: name of the file to save
Name of the output file
obsid: table column
Names of the objects
chi2: RawArray
......@@ -176,11 +175,11 @@ def save_table_best(obsid, chi2, chi2_red, variables, fluxes, filters,
column = Column(np_fluxes[:, index], name=name, unit='mJy')
best_model_table.add_column(column)
best_model_table.write(OUT_DIR + BEST_MODEL_FILE,
best_model_table.write(OUT_DIR + filename,
format='ascii.commented_header')
def dchi2_over_ds2(s, obs_fluxes, obs_errors, mod_fluxes):
def dchi2_over_ds2(s):
"""Function used to estimate the normalization factor in the SED fitting
process when upper limits are included in the dataset to fit (from Eq. A11
in Sawicki M. 2012, PASA, 124, 1008).
......@@ -212,17 +211,17 @@ def dchi2_over_ds2(s, obs_fluxes, obs_errors, mod_fluxes):
# The mask "lim" selects the filter(s) for which upper limits are given
# i.e., when obs_fluxes is >=0. and obs_errors = 9990 <= obs_errors < 0.
wlim = np.where((obs_errors >= -9990.)&(obs_errors < 0.))
wdata = np.where(obs_errors>=0.)
wlim = np.where((gbl_obs_errors >= -9990.)&(gbl_obs_errors < 0.))
wdata = np.where(gbl_obs_errors>=0.)
mod_fluxes_data = mod_fluxes[wdata]
mod_fluxes_lim = mod_fluxes[wlim]
mod_fluxes_data = gbl_mod_fluxes[wdata]
mod_fluxes_lim = gbl_mod_fluxes[wlim]
obs_fluxes_data = obs_fluxes[wdata]
obs_fluxes_lim = bs_fluxes[wlim]
obs_fluxes_data = gbl_obs_fluxes[wdata]
obs_fluxes_lim = gbl_obs_fluxes[wlim]
obs_errors_data = obs_errors[wdata]
obs_errors_lim = -obs_errors[wlim]
obs_errors_data = gbl_obs_errors[wdata]
obs_errors_lim = -gbl_obs_errors[wlim]
dchi2_over_ds_data = np.sum(
(obs_fluxes_data-s*mod_fluxes_data) *
......
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