pdf.py 3.96 KB
Newer Older
1 2 3 4 5 6 7
# -*- 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
8 9 10 11

import glob
from itertools import product
import matplotlib
12
from os import path
13 14 15 16 17 18

matplotlib.use('Agg')
import matplotlib.pyplot as plt
import multiprocessing as mp
import numpy as np
from pcigale.utils import read_table
19
from pcigale.analysis_modules.utils import Counter, nothreading
20 21 22 23 24 25 26 27 28 29 30


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.
31
    nothreading()
32
    gbl_counter = counter
33 34


35
def pdf(config, format, outdir):
36 37
    """Plot the PDF of analysed variables.
    """
38
    input_data = read_table(path.join(path.dirname(outdir), config.configuration['data_file']))
39 40 41 42
    pdf_vars = config.configuration['analysis_params']['variables']
    pdf_vars += [band for band in config.configuration['bands']
                 if band.endswith('_err') is False]

43
    items = list(product(input_data['id'], pdf_vars, [format], [outdir]))
44 45 46
    counter = Counter(len(items))
    with mp.Pool(processes=config.configuration['cores'], initializer=pool_initializer,
                 initargs=(counter,)) as pool:
47 48 49 50 51
        pool.starmap(_pdf_worker, items)
        pool.close()
        pool.join()


52
def _pdf_worker(obj_name, var_name, format, outdir):
53 54 55 56 57 58 59 60
    """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..
61 62
    outdir: string
        The absolute path to outdir
63 64

    """
65
    gbl_counter.inc()
66 67
    var_name = var_name.replace('/', '_')
    if var_name.endswith('_log'):
68 69
        fnames = glob.glob(f"{outdir}/{obj_name}_{var_name[:-4]}_chi2-block-"
                           f"*.npy")
70 71
        log = True
    else:
72
        fnames = glob.glob(f"{outdir}/{obj_name}_{var_name}_chi2-block-*.npy")
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
        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()
113 114 115
    figure.suptitle(f"Probability distribution function of {var_name} for "
                    f"{obj_name}")
    figure.savefig(f"{outdir}/{obj_name}_{var_name}_pdf.{format}")
116
    plt.close(figure)