utils.py 11.3 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 astropy.cosmology import WMAP7 as cosmo
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 = 'out/{}_{}_chi2-block-{}.npy'.format(obs.id, variable.replace('/',
                                                 '\/'), models.iblock)
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_dist):
34
35
36
37
38
39
40
    """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
    models and the actual observation. At low redshift, this can cause a
    discrepancy in the mass-dependent physical properties: ~0.35 dex at z=0.010
    vs 0.015 for instance. Therefore we correct these physical quantities by
    multiplying them by corr_dz.
41
42
43

    Parameters
    ----------
44
45
    model_z: float
        Redshift of the model.
46
47
    obs_dist: float
        Luminosity distance of the observed object.
48
49

    """
50
    if model_z == 0.:
51
52
        return (obs_dist / (10. * parsec)) ** 2.
    return (obs_dist / cosmo.luminosity_distance(model_z).value) ** 2.
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
255
        model = models.flux[band][wz]
        chi2 += ((flux - model * scaling) * inv_flux_err) ** 2.
256

257
    # Computation of the χ² from intensive properties
258
    for name, prop in obs.intprop.items():
259
260
        model = models.intprop[name][wz]
        chi2 += ((prop - model) * (1. / obs.intprops_err[name])) ** 2.
261

262
    # Computation of the χ² from extensive properties
263
    for name, prop in obs.extprop.items():
264
        # Multiplications are faster than divisions, so we directly use the
265
266
        # inverse error
        inv_prop_err = 1. / obs.extprop_err[name]
267
268
        model = models.extprop[name][wz]
        chi2 += ((prop - (scaling * model) * corr_dz) * inv_prop_err) ** 2.
269
270

    # Finally take the presence of upper limits into account
271
    if limits is True:
272
273
274
275
276
277
        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] -
                                 model * scaling) /
                                 (-np.sqrt(2.)*obs_error)))))
278
279
280
281
282
283
        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] -
                                 model * scaling) /
                                 (-np.sqrt(2.)*obs_error)))))
284
285

    return chi2, scaling
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310


def weighted_param(param, weights):
    """Compute the weighted mean and standard deviation of an array of data.

    Parameters
    ----------
    param: array
        Values of the parameters for the entire grid of models
    weights: array
        Weights by which to weight the parameter values

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

    """

    mean = np.average(param, weights=weights)
    std = np.sqrt(np.average((param-mean)**2, weights=weights))

    return (mean, std)