results.py 14 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# -*- coding: utf-8 -*-
# Copyright (C) 2017 Universidad de Antofagasta
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Médéric Boquien

"""These classes manage the results from the analysis. The main class
ResultsManager contains instances of BayesResultsManager and BestResultsManager
that contain the bayesian and best-fit estimates of the physical properties
along with the names of the parameters, which are proportional to the mass,
etc. Each of these classes contain a merge() method that allows to combine
results of the analysis with different blocks of models.
"""

from astropy.table import Table, Column
import numpy as np

17
from .utils import SharedArray
18
19
20
21
22
23
24
25
26
27
28
29
30


class BayesResultsManager(object):
    """This class contains the results of the bayesian estimates of the
    physical properties of the analysed objects. It is constructed from a
    ModelsManager instance, which provides the required information on the
    shape of the arrays. Because it can contain the partial results for only
    one block of models, we also store the sum of the model weights (that is
    the likelihood) so we can merge different instances to compute the combined
    estimates of the physical properties.

    """
    def __init__(self, models):
31
        nobs = len(models.obs)
32
        self.propertiesnames = models.allpropnames
33
34
35
36
37
38
        extpropnames = [prop for prop in models.obs.conf['analysis_params']['variables']
                        if (prop in models.allextpropnames or
                            prop[:-4] in models.allextpropnames)]
        intpropnames = [prop for prop in models.obs.conf['analysis_params']['variables']
                        if (prop in models.allintpropnames or
                            prop[:-4] in models.allintpropnames)]
39
        self.nproperties = len(intpropnames) + len(extpropnames)
40
41
42
43
44
45

        # Arrays where we store the data related to the models. For memory
        # efficiency reasons, we use RawArrays that will be passed in argument
        # to the pool. Each worker will fill a part of the RawArrays. It is
        # important that there is no conflict and that two different workers do
        # not write on the same section.
46
47
48
49
        self.intmean = {prop: SharedArray(nobs) for prop in intpropnames}
        self.interror = {prop: SharedArray(nobs) for prop in intpropnames}
        self.extmean = {prop: SharedArray(nobs) for prop in extpropnames}
        self.exterror = {prop: SharedArray(nobs) for prop in extpropnames}
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
        self.weight = SharedArray(nobs)

    @property
    def mean(self):
        return self._mean

    @mean.setter
    def mean(self, mean):
        self._mean = mean

    @property
    def error(self):
        return self._error

    @error.setter
    def error(self, error):
        self._error = error

    @property
    def weight(self):
        return self._weight

    @weight.setter
    def weight(self, weight):
        self._weight = weight
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92

    @staticmethod
    def merge(results):
        """Merge a list of partial results computed on individual blocks of
        models.

        Parameters
        ----------
        results: list of BayesResultsManager instances
            List of the partial results to be merged

        """
        if not (isinstance(results, list) and
                all(isinstance(result, BayesResultsManager)
                    for result in results)):
            raise TypeError("A list of BayesResultsManager is required.")

        merged = results[0]
93
        intmean = {prop: np.array([result.intmean[prop]
94
95
                                   for result in results])
                   for prop in merged.intmean}
96
        interror = {prop: np.array([result.interror[prop]
97
98
                                    for result in results])
                    for prop in merged.interror}
99
        extmean = {prop: np.array([result.extmean[prop]
100
101
                                   for result in results])
                   for prop in merged.extmean}
102
        exterror = {prop: np.array([result.exterror[prop]
103
104
                                    for result in results])
                    for prop in merged.exterror}
105
        weight = np.array([result.weight for result in results])
106
107
108
109

        totweight = np.sum(weight, axis=0)

        for prop in merged.intmean:
110
            merged.intmean[prop][:] = np.sum(
111
112
113
114
115
116
117
                intmean[prop] * weight, axis=0) / totweight

            # We compute the merged standard deviation by combining the
            # standard deviations for each block. See
            # http://stats.stackexchange.com/a/10445 where the number of
            # datapoints has been substituted with the weights. In short we
            # exploit the fact that Var(X) = E(Var(X)) + Var(E(X)).
118
119
            merged.interror[prop][:] = np.sqrt(np.sum(
                weight * (interror[prop]**2. + (intmean[prop]-merged.intmean[prop])**2), axis=0) / totweight)
120
121

        for prop in merged.extmean:
122
            merged.extmean[prop][:] = np.sum(
123
124
125
126
127
128
129
                extmean[prop] * weight, axis=0) / totweight

            # We compute the merged standard deviation by combining the
            # standard deviations for each block. See
            # http://stats.stackexchange.com/a/10445 where the number of
            # datapoints has been substituted with the weights. In short we
            # exploit the fact that Var(X) = E(Var(X)) + Var(E(X)).
130
131
            merged.exterror[prop][:] = np.sqrt(np.sum(
                weight * (exterror[prop]**2. + (extmean[prop]-merged.extmean[prop])**2), axis=0) / totweight)
132
133
134

        for prop in merged.extmean:
            if prop.endswith('_log'):
135
136
                merged.exterror[prop][:] = \
                    np.maximum(0.02, merged.exterror[prop])
137
            else:
138
139
140
                merged.exterror[prop][:] = \
                    np.maximum(0.05 * merged.extmean[prop],
                               merged.exterror[prop])
141

142
143
144
145
146
147
148
149
150
151
152
153
154
155
        return merged


class BestResultsManager(object):
    """This class contains the physical properties of the best fit of the
    analysed objects. It is constructed from a ModelsManager instance, which
    provides the required information on the shape of the arrays. Because it
    can contain the partial results for only one block of models, we also store
    the index so we can merge different instances to compute the best fit.

    """

    def __init__(self, models):
        self.obs = models.obs
156
        nobs = len(models.obs)
157
158
159
160
161
        # Arrays where we store the data related to the models. For memory
        # efficiency reasons, we use RawArrays that will be passed in argument
        # to the pool. Each worker will fill a part of the RawArrays. It is
        # important that there is no conflict and that two different workers do
        # not write on the same section.
162
        self.flux = {band: SharedArray(nobs) for band in models.obs.bands}
163
164
        allintpropnames = models.allintpropnames
        allextpropnames = models.allextpropnames
165
166
167
168
169
170
171
        self.intprop = {prop: SharedArray(nobs)
                        for prop in allintpropnames}
        self.extprop = {prop: SharedArray(nobs)
                        for prop in allextpropnames}
        self.chi2 = SharedArray(nobs)
        self.scaling = SharedArray(nobs)

172
        # We store the index as a float to work around python issue #10746
173
        self.index = SharedArray(nobs)
174
175

    @property
176
    def flux(self):
177
178
179
180
        """Returns a shared array containing the fluxes of the best fit for
        each observation.

        """
181
182
183
184
185
        return self._flux

    @flux.setter
    def flux(self, flux):
        self._flux = flux
186

187
    @property
188
    def intprop(self):
189
190
191
192
        """Returns a shared array containing the fluxes of the best fit for
        each observation.

        """
193
194
195
196
197
        return self._intprop

    @intprop.setter
    def intprop(self, intprop):
        self._intprop = intprop
198
199

    @property
200
    def extprop(self):
201
202
203
204
        """Returns a shared array containing the fluxes of the best fit for
        each observation.

        """
205
        return self._extprop
206

207
208
209
    @extprop.setter
    def extprop(self, extprop):
        self._extprop = extprop
210
211
212
213
214
215
216

    @property
    def chi2(self):
        """Returns a shared array containing the raw chi² of the best fit for
        each observation.

        """
217
        return self._chi2
218

219
220
221
222
    @chi2.setter
    def chi2(self, chi2):
        self._chi2 = chi2

223
224
225
226
227
228
    @property
    def index(self):
        """Returns a shared array containing the index of the best fit for each
        observation.

        """
229
        return self._index
230

231
232
233
234
    @index.setter
    def index(self, index):
        self._index = index

235
236
237
238
239
240
    @property
    def scaling(self):
        """Returns a shared array containing the scaling of the best fit for each
        observation.

        """
241
        return self._scaling
242
243
244
245
246

    @scaling.setter
    def scaling(self, scaling):
        self._scaling = scaling

247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
    @staticmethod
    def merge(results):
        """Merge a list of partial results computed on individual blocks of
        models.

        Parameters
        ----------
        results: list of BestResultsManager instances
            List of the partial results to be merged

        """
        if not (isinstance(results, list) and
                all(isinstance(result, BestResultsManager)
                    for result in results)):
            raise TypeError("A list of BestResultsManager is required.")

        if len(results) == 1:
            return results[0]

266
        best = np.argmin([result.chi2 for result in results], axis=0)
267
268

        merged = results[0]
269
270
        for iobs, bestidx in enumerate(best):
            for band in merged.flux:
271
272
                merged.flux[band][iobs] = \
                    results[bestidx].flux[band][iobs]
273
            for prop in merged.intprop:
274
275
                merged.intprop[prop][iobs] = \
                    results[bestidx].intprop[prop][iobs]
276
            for prop in merged.extprop:
277
278
                merged.extprop[prop][iobs] = \
                    results[bestidx].extprop[prop][iobs]
279
            merged.chi2[iobs] = results[bestidx].chi2[iobs]
280
            merged.scaling[iobs] = results[bestidx].scaling[iobs]
281
            merged.index[iobs] = results[bestidx].index[iobs]
282
283
284
285
286
287
288
289

        return merged

    def analyse_chi2(self):
        """Function to analyse the best chi^2 and find out what fraction of
         objects seems to be overconstrainted.

        """
290
        obs = [self.obs.table[obs].data for obs in self.obs.tofit]
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
        nobs = np.count_nonzero(np.isfinite(obs), axis=0)
        chi2_red = self.chi2 / (nobs - 1)
        # If low values of reduced chi^2, it means that the data are overfitted
        # Errors might be under-estimated or not enough valid data.
        print("\n{}% of the objects have chi^2_red~0 and {}% chi^2_red<0.5"
              .format(np.round((chi2_red < 1e-12).sum() / chi2_red.size, 1),
                      np.round((chi2_red < 0.5).sum() / chi2_red.size, 1)))


class ResultsManager(object):
    """This class contains the physical properties (best fit and bayesian) of
    the analysed objects. It is constructed from a ModelsManager instance,
    which provides the required information to initialise the instances of
    BestResultsManager and BayesResultsManager that store the results.

    """

    def __init__(self, models):
        self.conf = models.conf
        self.obs = models.obs
        self.params = models.params

        self.bayes = BayesResultsManager(models)
        self.best = BestResultsManager(models)

    @staticmethod
    def merge(results):
        """Merge a list of partial results computed on individual blocks of
        models.

        Parameters
        ----------
        results: list of ResultsManager instances
            List of the partial results to be merged

        """
        merged = results[0]
        merged.bayes = BayesResultsManager.merge([result.bayes
                                                  for result in results])
        merged.best = BestResultsManager.merge([result.best
                                                for result in results])

        return merged

    def save(self, filename):
        """Save the estimated values derived from the analysis of the PDF and
        the parameters associated with the best fit. A simple text file and a
        FITS file are generated.

        Parameters
        ----------
        filename:

        """
        table = Table()

        table.add_column(Column(self.obs.table['id'], name="id"))

349
        for prop in sorted(self.bayes.intmean):
350
            table.add_column(Column(self.bayes.intmean[prop],
351
                                    name="bayes."+prop))
352
            table.add_column(Column(self.bayes.interror[prop],
353
354
                                    name="bayes."+prop+"_err"))
        for prop in sorted(self.bayes.extmean):
355
            table.add_column(Column(self.bayes.extmean[prop],
356
                                    name="bayes."+prop))
357
            table.add_column(Column(self.bayes.exterror[prop],
358
                                    name="bayes."+prop+"_err"))
359
360

        table.add_column(Column(self.best.chi2, name="best.chi_square"))
361
        obs = [self.obs.table[obs].data for obs in self.obs.tofit]
362
363
364
365
        nobs = np.count_nonzero(np.isfinite(obs), axis=0)
        table.add_column(Column(self.best.chi2 / (nobs - 1),
                                name="best.reduced_chi_square"))

366
        for prop in sorted(self.best.intprop):
367
            table.add_column(Column(self.best.intprop[prop],
368
369
                                    name="best."+prop))
        for prop in sorted(self.best.extprop):
370
            table.add_column(Column(self.best.extprop[prop],
371
372
373
                                    name="best."+prop))

        for band in self.obs.bands:
374
375
376
377
            if band.startswith('line.'):
                unit = 'W/m^2'
            else:
                unit = 'mJy'
378
            table.add_column(Column(self.best.flux[band],
379
                                    name="best."+band, unit=unit))
380
381
382
383
384


        table.write("out/{}.txt".format(filename), format='ascii.fixed_width',
                    delimiter=None)
        table.write("out/{}.fits".format(filename), format='fits')