mock.py 4.08 KB
Newer Older
1 2 3 4 5 6 7 8
# -*- 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

9 10 11
from astropy.table import Table
import matplotlib
import sys
12
from os import path
13 14 15 16 17 18 19

matplotlib.use('Agg')
import matplotlib.pyplot as plt
import multiprocessing as mp
import numpy as np
import pkg_resources
from scipy import stats
20
from pcigale.analysis_modules.utils import Counter, nothread
21

22 23 24
# Name of the file containing the best models information
BEST_RESULTS = "results.fits"
MOCK_RESULTS = "results_mock.fits"
25

26 27 28 29 30 31 32 33 34 35 36 37 38 39

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


40
def mock(config, nologo, outdir):
41 42
    """Plot the comparison of input/output values of analysed variables.
    """
43 44
    best_results_file = path.abspath(path.join(outdir, BEST_RESULTS))
    mock_results_file = path.abspath(path.join(outdir, MOCK_RESULTS))
45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67

    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"))

68 69
    arguments = [(exact["best."+param], estimated["bayes."+param], param, logo, outdir)
                 for param in params]
70

71 72 73
    counter = Counter(len(arguments))
    with mp.Pool(processes=config.configuration['cores'], initializer=pool_initializer,
                 initargs=(counter,)) as pool:
74 75 76 77 78
        pool.starmap(_mock_worker, arguments)
        pool.close()
        pool.join()


79
def _mock_worker(exact, estimated, param, logo, outdir):
80 81 82 83 84 85 86 87 88 89 90 91
    """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.
92 93
    outdir: string
        The absolute path to outdir
94 95

    """
96
    gbl_counter.inc()
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
    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()
123
    plt.savefig('{}/mock_{}.pdf'.format(outdir, param))
124 125

    plt.close()