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

8
9
from functools import lru_cache

10
from astropy import log
11
from ...utils.cosmology import luminosity_distance
12
import numpy as np
13
from scipy import optimize
14
from scipy.constants import parsec
15
from scipy.special import erf
16

Médéric Boquien's avatar
Médéric Boquien committed
17
18
log.setLevel('ERROR')

 Hector Salas's avatar
Hector Salas committed
19

20
21
def save_chi2(obs, variable, models, chi2, values):
    """Save the chi² and the associated physocal properties
22
23

    """
24
25
    fname = f"out/{obs.id}_{variable.replace('/', '_')}_chi2-block-" \
            f"{models.iblock}.npy"
26
27
28
29
    data = np.memmap(fname, dtype=np.float64, mode='w+',
                     shape=(2, chi2.size))
    data[0, :] = chi2
    data[1, :] = values
30
31
32


@lru_cache(maxsize=None)
33
def compute_corr_dz(model_z, obs):
34
35
36
    """The mass-dependent physical properties are computed assuming the
    redshift of the model. However because we round the observed redshifts to
    two decimals, there can be a difference of 0.005 in redshift between the
37
38
    models and the actual observation. This causes two issues. First there is a
    difference in the luminosity distance. At low redshift, this can cause a
39
    discrepancy in the mass-dependent physical properties: ~0.35 dex at z=0.010
40
41
    vs 0.015 for instance. In addition, the 1+z dimming will be different.
    We compute here the correction factor for these two effects.
42
43
44

    Parameters
    ----------
45
46
    model_z: float
        Redshift of the model.
47
48
    obs: instance of the Observation class
        Object containing the distance and redshift of an object
49
50

    """
51
    return (obs.distance / luminosity_distance(model_z)) ** 2. * \
52
           (1. + model_z) / (1. + obs.redshift)
53
54


55
56
def dchi2_over_ds2(s, obsdata, obsdata_err, obslim, obslim_err, moddata,
                   modlim):
57
58
59
60
61
62
63
64
65
    """Function used to estimate the normalization factor in the SED fitting
    process when upper limits are included in the dataset to fit (from Eq. A11
    in Sawicki M. 2012, PASA, 124, 1008).

    Parameters
    ----------
    s: Float
        Contains value onto which we perform minimization = normalization
        factor
66
67
68
69
70
71
72
73
74
75
76
77
    obsdata: array
        Fluxes and extensive properties
    obsdata_err: array
        Errors on the fluxes and extensive properties
    obslim: array
        Fluxes and extensive properties upper limits
    obslim_err: array
        Errors on the fluxes and extensive properties upper limits
    moddata: array
        Model fluxes and extensive properties
    modlim: array
        Model fluxes and extensive properties for upper limits
78
79
80
81
82
83
84
85
86
87
88
89

   Returns
    -------
    func: Float
        Eq. A11 in Sawicki M. 2012, PASA, 124, 1008).

    """
    # We enter into this function if lim_flag = True.

    # The mask "data" selects the filter(s) for which measured fluxes are given
    # i.e., when obs_fluxes is >=0. and obs_errors >=0.
    # The mask "lim" selects the filter(s) for which upper limits are given
90
    # i.e., when obs_errors < 0
91
92
93
94
95
96
97
98
99
    dchi2_over_ds_data = np.sum((obsdata-s*moddata) * moddata/obsdata_err**2.)

    dchi2_over_ds_lim = np.sqrt(2./np.pi) * np.sum(
        modlim*np.exp(
            -np.square((obslim - s*modlim)/(np.sqrt(2)*obslim_err))
                     )/(
            obslim_err*(1. + erf((obslim - s*modlim)/(np.sqrt(2)*obslim_err)))
                       )
                                                  )
100
101
102
    func = dchi2_over_ds_data - dchi2_over_ds_lim

    return func
103

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

105
def _compute_scaling(models, obs, corr_dz, wz):
106
107
108
109
110
111
112
    """Compute the scaling factor to be applied to the model fluxes to best fit
    the observations. Note that we look over the bands to avoid the creation of
    an array of the same size as the model_fluxes array. Because we loop on the
    bands and not on the models, the impact on the performance should be small.

    Parameters
    ----------
113
114
    models: ModelsManagers class instance
        Contains the models (fluxes, intensive, and extensive properties).
115
116
117
    obs: Observation class instance
        Contains the fluxes, intensive properties, extensive properties and
        their errors, for a sigle observation.
118
119
120
    corr_dz: float
        Correction factor to scale the extensive properties to the right
        distance
121
122
123
124
    wz: slice
        Selection of the models at the redshift of the observation or all the
        redshifts in photometric-redshift mode.

125
126
127
128
129
    Returns
    -------
    scaling: array
        Scaling factors minimising the χ²
    """
130

131
    _ = list(models.flux.keys())[0]
132
133
    num = np.zeros_like(models.flux[_][wz])
    denom = np.zeros_like(models.flux[_][wz])
134
135

    for band, flux in obs.flux.items():
136
        # Multiplications are faster than divisions, so we directly use the
137
        # inverse error
138
        inv_err2 = 1. / obs.flux_err[band] ** 2.
139
140
141
        model = models.flux[band][wz]
        num += model * (flux * inv_err2)
        denom += model ** 2. * inv_err2
142
143

    for name, prop in obs.extprop.items():
144
        # Multiplications are faster than divisions, so we directly use the
145
        # inverse error
146
        inv_err2 = 1. / obs.extprop_err[name] ** 2.
147
148
149
        model = models.extprop[name][wz]
        num += model * (prop * inv_err2 * corr_dz)
        denom += model ** 2. * (inv_err2 * corr_dz ** 2.)
150
151
152

    return num/denom

153

154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
def _correct_scaling_ul(scaling, mod, obs, wz):
    """Correct the scaling factor when one or more fluxes and/or properties are
    upper limits.

    Parameters
    ----------
    scaling: array
        Contains the scaling factors of all the models
    mod: ModelsManager
        Contains the models
    obs: ObservationsManager
        Contains the observations
    wz: slice
        Selection of the models at the redshift of the observation or all the
        redshifts in photometric-redshift mode.
    """
    # Store the keys so we always read them in the same order
    fluxkeys = obs.flux.keys()
    fluxulkeys = obs.flux_ul.keys()
    extpropkeys = obs.extprop.keys()
    extpropulkeys = obs.extprop_ul.keys()

    # Put the fluxes and extensive properties in the same array for simplicity
    obsdata = [obs.flux[k] for k in fluxkeys]
    obsdata += [obs.extprop[k] for k in extpropkeys]
    obsdata_err = [obs.flux_err[k] for k in fluxkeys]
    obsdata_err += [obs.extprop_err[k] for k in extpropkeys]
    obslim = [obs.flux_ul[k] for k in fluxulkeys]
    obslim += [obs.extprop_ul[k] for k in extpropulkeys]
    obslim_err = [obs.flux_ul_err[k] for k in fluxulkeys]
    obslim_err += [obs.extprop_ul_err[k] for k in extpropulkeys]
    obsdata = np.array(obsdata)
    obsdata_err = np.array(obsdata_err)
    obslim = np.array(obslim)
    obslim_err = np.array(obslim_err)

    # We store views models at the right redshifts to avoid having SharedArray
    # recreate a numpy array for each model
    modflux = {k: mod.flux[k][wz] for k in mod.flux.keys()}
    modextprop = {k: mod.extprop[k][wz] for k in mod.extprop.keys()}
    for imod in range(scaling.size):
        moddata = [modflux[k][imod] for k in fluxkeys]
        moddata += [modextprop[k][imod] for k in extpropkeys]
        modlim = [modflux[k][imod] for k in fluxulkeys]
        modlim += [modextprop[k][imod] for k in extpropulkeys]
        moddata = np.array(moddata)
        modlim = np.array(modlim)
        scaling[imod] = optimize.root(dchi2_over_ds2, scaling[imod],
                                      args=(obsdata, obsdata_err,
                                            obslim, obslim_err,
                                            moddata, modlim)).x


207
def compute_chi2(models, obs, corr_dz, wz, lim_flag):
208
209
210
211
212
    """Compute the χ² of observed fluxes with respect to the grid of models. We
    take into account upper limits if need be. Note that we look over the bands
    to avoid the creation of an array of the same size as the model_fluxes
    array. Because we loop on the bands and not on the models, the impact on
    the performance should be small.
213
214
215

    Parameters
    ----------
216
217
    models: ModelsManagers class instance
        Contains the models (fluxes, intensive, and extensive properties).
218
219
220
    obs: Observation class instance
        Contains the fluxes, intensive properties, extensive properties and
        their errors, for a sigle observation.
221
222
    corr_dz: float
        Correction factor to scale the extensive properties to the right
223
        distance
224
225
226
    wz: slice
        Selection of the models at the redshift of the observation or all the
        redshifts in photometric-redshift mode.
227
    lim_flag: boolean
228
        Boolean indicating whether upper limits should be treated (True) or
229
230
231
232
233
234
235
236
237
        discarded (False)

    Returns
    -------
    chi2: array
        χ² for all the models in the grid
    scaling: array
        scaling of the models to obtain the minimum χ²
    """
238
    limits = lim_flag and (len(obs.flux_ul) > 0 or len(obs.extprop_ul) > 0)
239
    scaling = _compute_scaling(models, obs, corr_dz, wz)
240

241
242
    # Some observations may not have flux values in some filter(s), but
    # they can have upper limit(s).
243
244
    if limits is True:
        _correct_scaling_ul(scaling, models, obs, wz)
245

246
247
248
    # χ² of the comparison of each model to each observation.
    chi2 = np.zeros_like(scaling)

249
    # Computation of the χ² from fluxes
250
    for band, flux in obs.flux.items():
251
        # Multiplications are faster than divisions, so we directly use the
252
253
        # inverse error
        inv_flux_err = 1. / obs.flux_err[band]
254
        model = models.flux[band][wz]
255
        chi2 += ((model * scaling - flux) * inv_flux_err) ** 2.
256

Guang's avatar
Guang committed
257
    # Penalize det_alpha_ox which lie out of the user-set range
258
    if 'xray' in models.params.modules:
Guang's avatar
Guang committed
259
260
        # Get the model indice that have valid AGN component
        agn_idxs = np.where( models.extprop['agn.agn_intrin_Lnu_2500A'][wz]>0 )[0]
261
        # Calculate expected alpha_ox from Lnu_2500 (Just et al. 2007)
Guang's avatar
Guang committed
262
        exp_alpha_ox = -0.137*np.log10(models.extprop['agn.agn_intrin_Lnu_2500A'][wz][agn_idxs]*1e7*scaling[agn_idxs])+2.638
263
        # Calculate det_alpha_ox = alpha_ox - alpha_ox(Lnu_2500)
Guang's avatar
Guang committed
264
        det_alpha_ox = models.intprop['xray.alpha_ox'][wz][agn_idxs] - exp_alpha_ox
265
        # If det_alpha_ox out of range, set corresponding chi2 to nan
Guang's avatar
Guang committed
266
267
        nan_idxs = agn_idxs[ np.abs(det_alpha_ox) > models.intprop['xray.max_dev_alpha_ox'][wz][agn_idxs] ]
        chi2[nan_idxs] = np.nan
268

269
    # Computation of the χ² from intensive properties
270
    for name, prop in obs.intprop.items():
271
        model = models.intprop[name][wz]
272
        chi2 += ((model - prop) * (1. / obs.intprop_err[name])) ** 2.
273

274
    # Computation of the χ² from extensive properties
275
    for name, prop in obs.extprop.items():
276
        # Multiplications are faster than divisions, so we directly use the
277
278
        # inverse error
        inv_prop_err = 1. / obs.extprop_err[name]
279
        model = models.extprop[name][wz]
280
        chi2 += (((scaling * model) * corr_dz - prop) * inv_prop_err) ** 2.
281
282

    # Finally take the presence of upper limits into account
283
    if limits is True:
284
285
286
287
        for band, obs_error in obs.flux_ul_err.items():
            model = models.flux[band][wz]
            chi2 -= 2. * np.log(.5 *
                                (1. + erf(((obs.flux_ul[band] -
288
                                 model * scaling) / (np.sqrt(2.)*obs_error)))))
289
290
291
292
        for band, obs_error in obs.extprop_ul_err.items():
            model = models.extprop[band][wz]
            chi2 -= 2. * np.log(.5 *
                                (1. + erf(((obs.extprop_ul[band] -
293
                                 model * scaling) / (np.sqrt(2.)*obs_error)))))
294
295

    return chi2, scaling
296
297
298
299


def weighted_param(param, weights):
    """Compute the weighted mean and standard deviation of an array of data.
300
301
    Note that here we assume that the sum of the weights is normalised to 1.
    This simplifies and accelerates the computation.
302
303
304
305
306
307

    Parameters
    ----------
    param: array
        Values of the parameters for the entire grid of models
    weights: array
308
        Weights by which to weigh the parameter values
309
310
311
312
313
314
315
316
317
318

    Returns
    -------
    mean: float
        Weighted mean of the parameter values
    std: float
        Weighted standard deviation of the parameter values

    """

319
320
    mean = np.sum(param * weights)
    std = np.sqrt(np.sum((param - mean)**2 * weights))
321
322

    return (mean, std)