pdf.py 3.82 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
20
21
22
23
24
25
26
27
28


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
29

30
    gbl_counter = counter
31
32


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

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


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

    """
63
    gbl_counter.inc()
64
65
    var_name = var_name.replace('/', '_')
    if var_name.endswith('_log'):
66
67
        fnames = glob.glob(f"{outdir}/{obj_name}_{var_name[:-4]}_chi2-block-"
                           f"*.npy")
68
69
        log = True
    else:
70
        fnames = glob.glob(f"{outdir}/{obj_name}_{var_name}_chi2-block-*.npy")
71
72
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
        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()
111
112
113
    figure.suptitle(f"Probability distribution function of {var_name} for "
                    f"{obj_name}")
    figure.savefig(f"{outdir}/{obj_name}_{var_name}_pdf.{format}")
114
    plt.close(figure)