utils.py 9.74 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
51
52
    if model_z == 0.:
        return (obs_dist / (10. * parsec))**2.
    return (obs_dist / cosmo.luminosity_distance(model_z).value)**2.
53
54


55
def dchi2_over_ds2(s, obs_values, obs_errors, mod_values):
56
57
58
59
60
61
62
63
64
    """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
65
66
67
    obs_value: RawArray
        Contains observed fluxes for each filter and obseved extensive
        properties.
68
    obs_errors: RawArray
69
70
71
        Contains observed errors for each filter and extensive properties.
    model_values: RawArray
        Contains modeled fluxes for each filter and extensive properties.
72
73
74
75
76
77
78
79
80
81
82
83
84
85
    lim_flag: Boolean
        Tell whether we use upper limits (True) or not (False).

   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
86
    # i.e., when obs_errors < 0
87

88
    wlim = np.where(np.isfinite(obs_errors) & (obs_errors < 0.))
89
    wdata = np.where(obs_errors >= 0.)
90

Médéric Boquien's avatar
Médéric Boquien committed
91
92
    mod_values_data = mod_values[wdata]
    mod_values_lim = mod_values[wlim]
93

Médéric Boquien's avatar
Médéric Boquien committed
94
95
    obs_values_data = obs_values[wdata]
    obs_values_lim = obs_values[wlim]
96

97
98
    obs_errors_data = obs_errors[wdata]
    obs_errors_lim = -obs_errors[wlim]
99
100

    dchi2_over_ds_data = np.sum(
Médéric Boquien's avatar
Médéric Boquien committed
101
102
        (obs_values_data-s*mod_values_data) *
        mod_values_data/(obs_errors_data*obs_errors_data))
103
104

    dchi2_over_ds_lim = np.sqrt(2./np.pi)*np.sum(
Médéric Boquien's avatar
Médéric Boquien committed
105
        mod_values_lim*np.exp(
106
            -np.square(
Médéric Boquien's avatar
Médéric Boquien committed
107
                (obs_values_lim-s*mod_values_lim)/(np.sqrt(2)*obs_errors_lim)
108
109
110
111
                      )
                             )/(
            obs_errors_lim*(
                1.+erf(
Médéric Boquien's avatar
Médéric Boquien committed
112
                  (obs_values_lim-s*mod_values_lim)/(np.sqrt(2)*obs_errors_lim)
113
114
115
116
117
118
119
120
                      )
                           )
                               )
                                                )

    func = dchi2_over_ds_data - dchi2_over_ds_lim

    return func
121

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

123
def _compute_scaling(model_fluxes, model_propsmass, obs):
124
125
126
127
128
129
130
131
132
    """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
    ----------
    model_fluxes: array
        Fluxes of the models
133
134
    model_propsmass: array
        Extensive properties of the models to be fitted
135
136
137
    obs: Observation class instance
        Contains the fluxes, intensive properties, extensive properties and
        their errors, for a sigle observation.
138
139
140
141
142
    Returns
    -------
    scaling: array
        Scaling factors minimising the χ²
    """
143

144
145
    num = np.zeros(model_fluxes.shape[1])
    denom = np.zeros(model_fluxes.shape[1])
146
    for i in range(obs.fluxes.size):
147
148
149
        # Multiplications are faster than division, so we directly use the
        # inverse error
        inv_err2 = 1. / obs.fluxes_err[i] ** 2.
150
        if np.isfinite(obs.fluxes[i]) and obs.fluxes_err[i] > 0.:
151
            num += model_fluxes[i, :] * (obs.fluxes[i] * inv_err2)
152
153
            denom += np.square(model_fluxes[i, :] * (1./obs.fluxes_err[i]))
    for i in range(obs.extprops.size):
154
155
156
        # Multiplications are faster than division, so we directly use the
        # inverse error
        inv_err2 = 1. / obs.extprops_err[i] ** 2.
157
        if np.isfinite(obs.extprops[i]) and obs.extprops_err[i] > 0.:
158
159
            num += model_propsmass[i, :] * (obs.extprops[i] * inv_err2)
            denom += model_propsmass[i, :] ** 2. * inv_err2
160
161
162

    return num/denom

163

164
def compute_chi2(model_fluxes, model_props, model_propsmass, observation,
165
                 corr_dz, lim_flag):
166
167
168
169
170
    """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.
171
172
173
174
175

    Parameters
    ----------
    model_fluxes: array
        2D grid containing the fluxes of the models
176
177
178
179
180
181
182
    model_props: array
        2D grid containing the intensive properties of the models
    model_propsmass: array
        2D grid containing the extensive properties of the models
    observation: Class
        Class instance containing the fluxes, intensive properties, extensive
        properties and their errors, for a sigle observation.
183
184
    corr_dz: correction factor to scale the extensive properties to the right
        distance
185
    lim_flag: boolean
186
        Boolean indicating whether upper limits should be treated (True) or
187
188
189
190
191
192
193
194
195
        discarded (False)

    Returns
    -------
    chi2: array
        χ² for all the models in the grid
    scaling: array
        scaling of the models to obtain the minimum χ²
    """
196
197
198
    limits = lim_flag and np.any(observation.fluxes_err +
                                 observation.extprops_err <= 0.)
    scaling = _compute_scaling(model_fluxes, model_propsmass, observation)
199

200
201
202
203
204
205
    obs_fluxes = observation.fluxes
    obs_fluxes_err = observation.fluxes_err
    obs_props = observation.intprops
    obs_props_err = observation.intprops_err
    obs_propsmass = observation.extprops
    obs_propsmass_err = observation.extprops_err
206

207
208
209
    # Some observations may not have flux values in some filter(s), but
    # they can have upper limit(s).
    if limits == True:
210
211
        obs_values = np.concatenate((obs_fluxes, obs_propsmass))
        obs_values_err = np.concatenate((obs_fluxes_err, obs_propsmass_err))
212
        model_values = np.concatenate((model_fluxes, model_propsmass))
213
        for imod in range(scaling.size):
214
            scaling[imod] = optimize.root(dchi2_over_ds2, scaling[imod],
215
216
                                          args=(obs_values, obs_values_err,
                                                model_values[:, imod])).x
217
218

    # χ² of the comparison of each model to each observation.
219
    chi2 = np.zeros(model_fluxes.shape[1])
220
    for i in range(obs_fluxes.size):
221
222
223
        if np.isfinite(obs_fluxes[i]) and obs_fluxes_err[i] > 0.:
            chi2 += np.square((obs_fluxes[i] - model_fluxes[i, :] * scaling) *
                              (1./obs_fluxes_err[i]))
224

225
226
    for i in range(obs_propsmass.size):
        if np.isfinite(obs_propsmass[i]):
227
228
229
            chi2 += np.square((obs_propsmass[i] - corr_dz * (scaling *
                               model_propsmass[i, :])) *
                              (1./obs_propsmass_err[i]))
230
231
232
233
234

    for i in range(obs_props.size):
        if np.isfinite(obs_props[i]):
            chi2 += np.square((obs_props[i] - model_props[i, :]) *
                              (1./obs_props_err[i]))
235
    # they can have upper limit(s).
236
    if limits == True:
237
        for i, obs_error in enumerate(obs_fluxes_err):
238
            if obs_error < 0.:
239
240
                chi2 -= 2. * np.log(.5 *
                                    (1. + erf(((obs_fluxes[i] -
241
242
                                     model_fluxes[i, :] * scaling) /
                                     (-np.sqrt(2.)*obs_fluxes_err[i])))))
243
244

    return chi2, scaling
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269


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)