workers.py 9.07 KB
Newer Older
1 2 3
# -*- coding: utf-8 -*-
# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
# Copyright (C) 2013-2014 Institute of Astronomy
4
# Copyright (C) 2014 Laboratoire d'Astrophysique de Marseille, AMU
5 6
# Copyright (C) 2014 Yannick Roehlly <yannick@iaora.eu>
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
7
# Author: Yannick Roehlly, Médéric Boquien & Denis Burgarella
8

9
from copy import deepcopy
10

11
import numpy as np
12

13
from ..utils import nothread
14
from .utils import save_chi2, compute_corr_dz, compute_chi2, weighted_param
15
from ...warehouse import SedWarehouse
16

Médéric Boquien's avatar
Médéric Boquien committed
17

18
def init_sed(models, counter):
19
    """Initializer of the pool of processes to share variables between workers.
20 21 22

    Parameters
    ----------
23 24
    models: ModelsManagers
        Manages the storage of the computed models (fluxes and properties).
25 26
    counter: Counter class object
        Counter for the number of models computed
27 28

    """
29
    global gbl_warehouse, gbl_models, gbl_counter
30

31 32 33 34
    # Limit the number of threads to 1 if we use MKL in order to limit the
    # oversubscription of the CPU/RAM.
    nothread()

35
    gbl_warehouse = SedWarehouse()
36

37
    gbl_models = models
38
    gbl_counter = counter
Médéric Boquien's avatar
Médéric Boquien committed
39

 Hector Salas's avatar
Hector Salas committed
40

41
def init_analysis(models, results, counter):
42
    """Initializer of the pool of processes to share variables between workers.
43 44 45

    Parameters
    ----------
46 47 48 49
    models: ModelsManagers
        Manages the storage of the computed models (fluxes and properties).
    results: ResultsManager
        Contains the estimates and errors on the properties.
50 51
    counter: Counter class object
        Counter for the number of objects analysed
52 53

    """
54
    global gbl_models, gbl_obs, gbl_results, gbl_counter
55

56
    gbl_models = models
57
    gbl_obs = models.obs
58
    gbl_results = results
59
    gbl_counter = counter
60 61


62
def init_bestfit(conf, params, observations, results, counter):
63
    """Initializer of the pool of processes to share variables between workers.
64

65 66 67 68 69 70 71 72 73 74
    Parameters
    ----------
    conf: dict
        Contains the pcigale.ini configuration.
    params: ParametersManager
        Manages the parameters from a 1D index.
    observations: astropy.Table
        Contains the observations including the filter names.
    results: ResultsManager
        Contains the estimates and errors on the properties.
75 76
    counter: Counter class object
        Counter for the number of objects analysed
77

78
    """
79
    global gbl_warehouse, gbl_conf, gbl_params, gbl_obs
80
    global gbl_results, gbl_counter
81

82
    gbl_warehouse = SedWarehouse()
83

84 85 86 87
    gbl_conf = conf
    gbl_params = params
    gbl_obs = observations
    gbl_results = results
88
    gbl_counter = counter
Médéric Boquien's avatar
Médéric Boquien committed
89

 Hector Salas's avatar
Hector Salas committed
90

91 92 93
def sed(idx, midx):
    """Worker process to retrieve a SED and affect the relevant data to an
    instance of ModelsManager.
94 95 96

    Parameters
    ----------
97
    idx: int
98 99 100
        Index of the model within the current block of models.
    midx: int
        Global index of the model.
101 102

    """
103 104
    sed = gbl_warehouse.get_sed(gbl_models.params.modules,
                                gbl_models.params.from_index(midx))
105

106
    if 'sfh.age' in sed.info and sed.info['sfh.age'] > sed.info['universe.age']:
107
        for band in gbl_models.flux:
108
            gbl_models.flux[band][idx] = np.nan
109
        for prop in gbl_models.extprop:
110
            gbl_models.extprop[prop][idx] = np.nan
111
        for prop in gbl_models.intprop:
112
            gbl_models.intprop[prop][idx] = np.nan
113
    else:
114
        for band in gbl_models.flux:
115
            gbl_models.flux[band][idx] = sed.compute_fnu(band)
116
        for prop in gbl_models.extprop:
117
            gbl_models.extprop[prop][idx] = sed.info[prop]
118
        for prop in gbl_models.intprop:
119
            gbl_models.intprop[prop][idx] = sed.info[prop]
120

121
    gbl_counter.inc()
122

 Hector Salas's avatar
Hector Salas committed
123

124
def analysis(idx, obs):
125 126
    """Worker process to analyse the PDF and estimate parameters values and
    store them in an instance of ResultsManager.
127 128 129

    Parameters
    ----------
130 131
    idx: int
        Index of the observation. This is necessary to put the computed values
132
        at the right location in the ResultsManager.
133 134 135 136
    obs: row
        Input data for an individual object

    """
137 138
    np.seterr(invalid='ignore')

139
    if obs.redshift >= 0.:
140 141
        # We pick the the models with the closest redshift using a slice to
        # work on views of the arrays and not on copies to save on RAM.
142 143 144
        z = np.array(
            gbl_models.conf['sed_modules_params']['redshifting']['redshift'])
        wz = slice(np.abs(obs.redshift-z).argmin(), None, z.size)
145
        corr_dz = compute_corr_dz(z[wz.start], obs.distance)
146 147
    else:  # We do not know the redshift so we use the full grid
        wz = slice(0, None, 1)
148
        corr_dz = 1.
149

150
    chi2, scaling = compute_chi2(gbl_models, obs, corr_dz, wz,
151
                                 gbl_models.conf['analysis_params']['lim_flag'])
152

153
    if np.any(chi2 < -np.log(np.finfo(np.float64).tiny) * 2.):
154 155
        # We use the exponential probability associated with the χ² as
        # likelihood function.
156
        likelihood = np.exp(-chi2 / 2.)
157 158 159 160
        wlikely = np.where(np.isfinite(likelihood))
        # If all the models are valid, it is much more efficient to use a slice
        if likelihood.size == wlikely[0].size:
            wlikely = slice(None, None)
161 162 163
        likelihood = likelihood[wlikely]
        scaling_l = scaling[wlikely]

164
        gbl_results.bayes.weight[idx] = np.nansum(likelihood)
165

Médéric Boquien's avatar
Médéric Boquien committed
166 167
        # We compute the weighted average and standard deviation using the
        # likelihood as weight.
168 169
        for prop in gbl_results.bayes.intmean:
            if prop.endswith('_log'):
170
                values = gbl_models.intprop[prop[:-4]][wz]
171 172
                _ = np.log10
            else:
173
                values = gbl_models.intprop[prop][wz]
174
                _ = lambda x: x
175
            mean, std = weighted_param(_(values[wlikely]), likelihood)
176 177
            gbl_results.bayes.intmean[prop][idx] = mean
            gbl_results.bayes.interror[prop][idx] = std
178 179
            if gbl_models.conf['analysis_params']['save_chi2'] is True:
                save_chi2(obs, prop, gbl_models, chi2, values)
180

181 182
        for prop in gbl_results.bayes.extmean:
            if prop.endswith('_log'):
183
                values = gbl_models.extprop[prop[:-4]][wz]
184
                _ = np.log10
185
            else:
186
                values = gbl_models.extprop[prop][wz]
187
                _ = lambda x: x
188 189
            mean, std = weighted_param(_(values[wlikely] * scaling_l * corr_dz),
                                       likelihood)
190 191
            gbl_results.bayes.extmean[prop][idx] = mean
            gbl_results.bayes.exterror[prop][idx] = std
192
            if gbl_models.conf['analysis_params']['save_chi2'] is True:
193 194
                save_chi2(obs, prop, gbl_models, chi2,
                          values * scaling * corr_dz)
195

196 197 198 199 200 201 202 203 204
        for band in gbl_results.bayes.fluxmean:
            values = gbl_models.flux[band][wz]
            mean, std = weighted_param(values[wlikely] * scaling_l,
                                       likelihood)
            gbl_results.bayes.fluxmean[band][idx] = mean
            gbl_results.bayes.fluxerror[band][idx] = std
            if gbl_models.conf['analysis_params']['save_chi2'] is True:
                save_chi2(obs, band, gbl_models, chi2, values * scaling)

205 206
        best_idx_z = np.nanargmin(chi2)
        gbl_results.best.chi2[idx] = chi2[best_idx_z]
207
        gbl_results.best.scaling[idx] = scaling[best_idx_z]
208 209 210 211
        gbl_results.best.index[idx] = (wz.start + best_idx_z*wz.step +
                                       gbl_models.block.start)
    else:
        # It sometimes happens because models are older than the Universe's age
212 213 214
        print("No suitable model found for the object {}. It may be that "
              "models are older than the Universe or that your chi² are very "
              "large.".format(obs.id))
215

216
    gbl_counter.inc()
 Hector Salas's avatar
Hector Salas committed
217

218 219 220 221 222 223 224 225 226 227 228 229 230 231 232
def bestfit(oidx, obs):
    """Worker process to compute and save the best fit.

    Parameters
    ----------
    oidx: int
        Index of the observation. This is necessary to put the computed values
        at the right location in the ResultsManager.
    obs: row
        Input data for an individual object

    """
    np.seterr(invalid='ignore')

    best_index = int(gbl_results.best.index[oidx])
233 234 235 236

    # We compute the model at the exact redshift not to have to correct for the
    # difference between the object and the grid redshifts.
    params = deepcopy(gbl_params.from_index(best_index))
237
    if obs.redshift >= 0.:
238
        params[gbl_params.modules.index('redshifting')]['redshift'] = obs.redshift
239
    sed = gbl_warehouse.get_sed(gbl_params.modules, params)
240

241
    corr_dz = compute_corr_dz(obs.redshift, obs.distance)
242
    scaling = gbl_results.best.scaling[oidx]
243 244

    for band in gbl_results.best.flux:
245
        gbl_results.best.flux[band][oidx] = sed.compute_fnu(band) * scaling
246 247

    for prop in gbl_results.best.intprop:
248
        gbl_results.best.intprop[prop][oidx] = sed.info[prop]
249 250

    for prop in gbl_results.best.extprop:
251
        gbl_results.best.extprop[prop][oidx] = sed.info[prop] * scaling \
252
                                                   * corr_dz
253 254

    if gbl_conf['analysis_params']["save_best_sed"]:
255
        sed.to_fits('out/{}'.format(obs.id), scaling)
256

257
    gbl_counter.inc()