Change from __init__.py to own module.

pcigale_plots/parsers/chi2.py

+
+import glob
+from itertools import product
+
+import matplotlib
+
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import multiprocessing as mp
+import numpy as np
+from pcigale.utils import read_table
+
+
+def chi2(config):
+    """Plot the χ² values of analysed variables.
+    """
+    input_data = read_table(config.configuration['data_file'])
+    chi2_vars = config.configuration['analysis_params']['variables']
+    chi2_vars += [band for band in config.configuration['bands']
+                  if band.endswith('_err') is False]
+
+    with mp.Pool(processes=config.configuration['cores']) as pool:
+        items = product(input_data['id'], chi2_vars)
+        pool.starmap(_chi2_worker, items)
+        pool.close()
+        pool.join()
+
+
+def _chi2_worker(obj_name, var_name):
+    """Plot the reduced χ² associated with a given analysed variable
+
+    Parameters
+    ----------
+    obj_name: string
+        Name of the object.
+    var_name: string
+        Name of the analysed variable..
+
+    """
+    figure = plt.figure()
+    ax = figure.add_subplot(111)
+
+    var_name = var_name.replace('/', '_')
+    fnames = glob.glob("out/{}_{}_chi2-block-*.npy".format(obj_name, var_name))
+    for fname in fnames:
+        data = np.memmap(fname, dtype=np.float64)
+        data = np.memmap(fname, dtype=np.float64, shape=(2, data.size // 2))
+        ax.scatter(data[1, :], data[0, :], color='k', s=.1)
+    ax.set_xlabel(var_name)
+    ax.set_ylabel("Reduced $\chi^2$")
+    ax.set_ylim(0., )
+    ax.minorticks_on()
+    figure.suptitle("Reduced $\chi^2$ distribution of {} for {}."
+                    .format(var_name, obj_name))
+    figure.savefig("out/{}_{}_chi2.pdf".format(obj_name, var_name))
+    plt.close(figure)

pcigale_plots/parsers/mock.py

+from astropy.table import Table
+import matplotlib
+import sys
+
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import multiprocessing as mp
+import numpy as np
+import pkg_resources
+from scipy import stats
+
+
+def mock(config, mock_results_file, best_results_file, nologo):
+    """Plot the comparison of input/output values of analysed variables.
+    """
+
+    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"))
+
+    arguments = ((exact["best."+param], estimated["bayes."+param], param, logo)
+                 for param in params)
+
+    with mp.Pool(processes=config.configuration['cores']) as pool:
+        pool.starmap(_mock_worker, arguments)
+        pool.close()
+        pool.join()
+
+
+def _mock_worker(exact, estimated, param, logo):
+    """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.
+
+    """
+
+    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()
+    plt.savefig('out/mock_{}.pdf'.format(param))
+
+    plt.close()

pcigale_plots/parsers/pdf.py

+
+import glob
+from itertools import product
+import matplotlib
+
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import multiprocessing as mp
+import numpy as np
+from pcigale.utils import read_table
+
+
+def pdf(config):
+    """Plot the PDF of analysed variables.
+    """
+    input_data = read_table(config.configuration['data_file'])
+    pdf_vars = config.configuration['analysis_params']['variables']
+    pdf_vars += [band for band in config.configuration['bands']
+                 if band.endswith('_err') is False]
+
+    with mp.Pool(processes=config.configuration['cores']) as pool:
+        items = product(input_data['id'], pdf_vars)
+        pool.starmap(_pdf_worker, items)
+        pool.close()
+        pool.join()
+
+
+def _pdf_worker(obj_name, var_name):
+    """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..
+
+    """
+    var_name = var_name.replace('/', '_')
+    if var_name.endswith('_log'):
+        fnames = glob.glob("out/{}_{}_chi2-block-*.npy".format(obj_name,
+                                                               var_name[:-4]))
+        log = True
+    else:
+        fnames = glob.glob("out/{}_{}_chi2-block-*.npy".format(obj_name,
+                                                               var_name))
+        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()
+    figure.suptitle("Probability distribution function of {} for {}"
+                    .format(var_name, obj_name))
+    figure.savefig("out/{}_{}_pdf.pdf".format(obj_name, var_name))
+    plt.close(figure)