utils.py 11.4 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

257
    # Computation of the χ² from intensive properties
258
    for name, prop in obs.intprop.items():
259
        model = models.intprop[name][wz]
260
        chi2 += ((model - prop) * (1. / obs.intprop_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
        model = models.extprop[name][wz]
268
        chi2 += (((scaling * model) * corr_dz - prop) * inv_prop_err) ** 2.
269 270

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

    return chi2, scaling
284 285 286 287


def weighted_param(param, weights):
    """Compute the weighted mean and standard deviation of an array of data.
288 289
    Note that here we assume that the sum of the weights is normalised to 1.
    This simplifies and accelerates the computation.
290 291 292 293 294 295

    Parameters
    ----------
    param: array
        Values of the parameters for the entire grid of models
    weights: array
296
        Weights by which to weigh the parameter values
297 298 299 300 301 302 303 304 305 306

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

    """

307 308
    mean = np.sum(param * weights)
    std = np.sqrt(np.sum((param - mean)**2 * weights))
309 310

    return (mean, std)