results.py 11.6 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
31
32
33
34
35
36
37
38
39
40
41


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):
        self.nobs = len(models.obs)
        self.propertiesnames = models.propertiesnames
        self.massproportional = models.massproportional.\
                                   intersection(models.propertiesnames)
        self.nproperties = len(models.propertiesnames)

        # 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.
42
43
44
        self._means = SharedArray((self.nobs, self.nproperties))
        self._errors = SharedArray((self.nobs, self.nproperties))
        self._weights = SharedArray((self.nobs))
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

    @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.")

        means = np.array([result.means for result in results])
        errors = np.array([result.errors for result in results])
64
        weights = np.array([result.weights for result in results])[..., None]
65
66

        merged = results[0]
67
68
        merged._means = SharedArray((merged.nobs, merged.nproperties))
        merged._errors = SharedArray((merged.nobs, merged.nproperties))
69
70
71
        merged._weights = None

        sumweights = np.sum(weights, axis=0)
72
73
74

        merged.means[:] = np.sum(means * weights, axis=0) / sumweights

75
        # We compute the merged standard deviation by combining the standard
76
        # deviations for each block. See http://stats.stackexchange.com/a/10445
77
        # where the number of datapoints has been substituted with the weights.
78
79
80
81
        # In short we exploit the fact that Var(X) = E(Var(X)) + Var(E(X)).
        merged.errors[:] = np.sqrt(np.sum(weights * (errors**2 +
                                          (means-merged.means)**2), axis=0) /
                                   sumweights)
82

83
84
        for i, variable in enumerate(merged.propertiesnames):
            if variable.endswith('_log'):
85
                merged.errors[:, i] = np.maximum(0.02, merged.errors[:, i])
86
            else:
87
88
                merged.errors[:, i] = np.maximum(0.05 * merged.means[:, i],
                                                 merged.errors[:, i])
89

90
91
92
93
94
95
96
97
        return merged

    @property
    def means(self):
        """Returns a shared array containing the weighted means for each
        physical property and each observation.

        """
98
        return self._means.data
99
100
101
102
103
104
105

    @property
    def errors(self):
        """Returns a shared array containing the weighted standard deviation
         for each physical property and each observation.

        """
106
        return self._errors.data
107
108
109
110
111
112
113

    @property
    def weights(self):
        """Returns a shared array containing the some of the weights for each
        each observation.

        """
114
        return self._weights.data
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141


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
        self.nbands = len(models.obs.bands)
        self.nobs = len(models.obs)
        self.propertiesnames = models.allpropertiesnames
        self.massproportional = models.massproportional
        self.nproperties = len(models.allpropertiesnames)

        self._fluxes_shape = (self.nobs, self.nbands)
        self._properties_shape = (self.nobs, self.nproperties)

        # 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.
142
143
144
        self._fluxes = SharedArray(self._fluxes_shape)
        self._properties = SharedArray(self._properties_shape)
        self._chi2 = SharedArray(self.nobs)
145
        # We store the index as a float to work around python issue #10746
146
        self._index = SharedArray(self.nobs)
147
148
149
150
151
152
153

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

        """
154
        return self._fluxes.data
155
156
157
158
159
160
161

    @property
    def properties(self):
        """Returns a shared array containing the physical properties of the
        best fit for each observation.

        """
162
        return self._properties.data
163
164
165
166
167
168
169

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

        """
170
        return self._chi2.data
171
172
173
174
175
176
177

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

        """
178
        return self._index.data
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

    @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]

        fluxes = np.array([result.fluxes for result in results])
        properties = np.array([result.properties for result in results])
        chi2 = np.array([result.chi2 for result in results])
        index = np.array([result.index for result in results])

        merged = results[0]
205
206
207
        merged._fluxes = SharedArray((merged.nobs, merged.nbands))
        merged._properties = SharedArray((merged.nobs, merged.nproperties))
        merged._chi2 = SharedArray(merged.nobs)
208
        # We store the index as a float to work around python issue #10746
209
        merged._index = SharedArray(merged.nobs)
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
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
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305

        for iobs, bestidx in enumerate(np.argmin(chi2, axis=0)):
            merged.fluxes[iobs, :] = fluxes[bestidx, iobs, :]
            merged.properties[iobs, :] = properties[bestidx, iobs, :]
            merged.chi2[iobs] = chi2[bestidx, iobs]
            merged.index[iobs] = index[bestidx, iobs]

        return merged

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

        """
        obs = [self.obs.table[obs].data for obs in self.obs.bands]
        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"))

        for idx, name in enumerate(self.bayes.propertiesnames):
            table.add_column(Column(self.bayes.means[:, idx],
                                    name="bayes."+name))
            table.add_column(Column(self.bayes.errors[:, idx],
                                    name="bayes."+name+"_err"))

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

        for idx, name in enumerate(self.best.propertiesnames):
            table.add_column(Column(self.best.properties[:, idx],
                                    name="best."+name))

        for idx, name in enumerate(self.obs.bands):
            table.add_column(Column(self.best.fluxes[:, idx],
                                    name="best."+name, unit='mJy'))

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