Passed
Pull Request — master (#4)
by
unknown
11:38
created

hawc_hal.HAL.HAL.display_fit()   B

Complexity

Conditions 3

Size

Total Lines 82
Code Lines 38

Duplication

Lines 0
Ratio 0 %

Importance

Changes 0
Metric Value
eloc 38
dl 0
loc 82
rs 8.968
c 0
b 0
f 0
cc 3
nop 2

How to fix   Long Method   

Long Method

Small methods make your code easier to understand, in particular if combined with a good name. Besides, if your method is small, finding a good name is usually much easier.

For example, if you find yourself adding comments to a method's body, this is usually a good sign to extract the commented part to a new method, and use the comment as a starting point when coming up with a good name for this new method.

Commonly applied refactorings include:

1
import collections
0 ignored issues
show
Coding Style introduced by
This module should have a docstring.

The coding style of this project requires that you add a docstring to this code element. Below, you find an example for methods:

class SomeClass:
    def some_method(self):
        """Do x and return foo."""

If you would like to know more about docstrings, we recommend to read PEP-257: Docstring Conventions.

Loading history...
Coding Style Naming introduced by
Module name "HAL" doesn't conform to snake_case naming style ('([a-z_][a-z0-9_]*)$' pattern)

This check looks for invalid names for a range of different identifiers.

You can set regular expressions to which the identifiers must conform if the defaults do not match your requirements.

If your project includes a Pylint configuration file, the settings contained in that file take precedence.

To find out more about Pylint, please refer to their site.

Loading history...
2
import numpy as np
0 ignored issues
show
introduced by
Unable to import 'numpy'
Loading history...
3
import healpy as hp
0 ignored issues
show
introduced by
Unable to import 'healpy'
Loading history...
4
import matplotlib.pyplot as plt
0 ignored issues
show
introduced by
Unable to import 'matplotlib.pyplot'
Loading history...
5
import copy
0 ignored issues
show
introduced by
standard import "import copy" should be placed before "import numpy as np"
Loading history...
6
from astropy.convolution import Gaussian2DKernel
0 ignored issues
show
introduced by
Unable to import 'astropy.convolution'
Loading history...
7
from astropy.convolution import convolve_fft as convolve
0 ignored issues
show
introduced by
Unable to import 'astropy.convolution'
Loading history...
8
from scipy.stats import poisson
0 ignored issues
show
introduced by
Unable to import 'scipy.stats'
Loading history...
9
10
from threeML.plugin_prototype import PluginPrototype
0 ignored issues
show
introduced by
Unable to import 'threeML.plugin_prototype'
Loading history...
11
from threeML.plugins.gammaln import logfactorial
0 ignored issues
show
introduced by
Unable to import 'threeML.plugins.gammaln'
Loading history...
12
from threeML.parallel import parallel_client
0 ignored issues
show
introduced by
Unable to import 'threeML.parallel'
Loading history...
13
from threeML.io.progress_bar import progress_bar
0 ignored issues
show
introduced by
Unable to import 'threeML.io.progress_bar'
Loading history...
14
15
from astromodels import Parameter
0 ignored issues
show
introduced by
Unable to import 'astromodels'
Loading history...
16
17
from maptree import map_tree_factory
0 ignored issues
show
Coding Style introduced by
Relative import 'maptree', should be 'hawc_hal.maptree'
Loading history...
18
from response import hawc_response_factory
0 ignored issues
show
Coding Style introduced by
Relative import 'response', should be 'hawc_hal.response'
Loading history...
19
from convolved_source import ConvolvedPointSource, \
0 ignored issues
show
Coding Style introduced by
Relative import 'convolved_source', should be 'hawc_hal.convolved_source'
Loading history...
20
    ConvolvedExtendedSource3D, ConvolvedExtendedSource2D, ConvolvedSourcesContainer
21
from healpix_handling import FlatSkyToHealpixTransform
0 ignored issues
show
Coding Style introduced by
Relative import 'healpix_handling', should be 'hawc_hal.healpix_handling'
Loading history...
22
from healpix_handling import SparseHealpix
0 ignored issues
show
Coding Style introduced by
Relative import 'healpix_handling', should be 'hawc_hal.healpix_handling'
Loading history...
23
from healpix_handling import get_gnomonic_projection
0 ignored issues
show
Coding Style introduced by
Relative import 'healpix_handling', should be 'hawc_hal.healpix_handling'
Loading history...
24
from psf_fast import PSFConvolutor
0 ignored issues
show
Coding Style introduced by
Relative import 'psf_fast', should be 'hawc_hal.psf_fast'
Loading history...
25
from log_likelihood import log_likelihood
0 ignored issues
show
Coding Style introduced by
Relative import 'log_likelihood', should be 'hawc_hal.log_likelihood'
Loading history...
26
from util import ra_to_longitude
0 ignored issues
show
Coding Style introduced by
Relative import 'util', should be 'hawc_hal.util'
Loading history...
27
28
29
class HAL(PluginPrototype):
0 ignored issues
show
Coding Style introduced by
This class should have a docstring.

The coding style of this project requires that you add a docstring to this code element. Below, you find an example for methods:

class SomeClass:
    def some_method(self):
        """Do x and return foo."""

If you would like to know more about docstrings, we recommend to read PEP-257: Docstring Conventions.

Loading history...
best-practice introduced by
Too many instance attributes (16/7)
Loading history...
30
31
    def __init__(self, name, maptree, response_file, roi, flat_sky_pixels_sizes=0.17):
0 ignored issues
show
best-practice introduced by
Too many arguments (6/5)
Loading history...
32
33
        # Store ROI
34
        self._roi = roi
35
36
        # Set up the flat-sky projection
37
38
        self._flat_sky_projection = roi.get_flat_sky_projection(flat_sky_pixels_sizes)
39
40
        # Read map tree (data)
41
42
        self._maptree = map_tree_factory(maptree, roi=roi)
43
44
        # Read detector response_file
45
46
        self._response = hawc_response_factory(response_file)
47
48
        # Use a renormalization of the background as nuisance parameter
49
        # NOTE: it is fixed to 1.0 unless the user explicitly sets it free (experimental)
50
        self._nuisance_parameters = collections.OrderedDict()
51
        self._nuisance_parameters['%s_bkg_renorm' % name] = Parameter('%s_bkg_renorm' % name, 1.0,
52
                                                                      min_value=0.5, max_value=1.5,
53
                                                                      delta=0.01,
54
                                                                      desc="Renormalization for background map",
55
                                                                      free=False,
56
                                                                      is_normalization=False)
57
58
        # Instance parent class
59
60
        super(HAL, self).__init__(name, self._nuisance_parameters)
61
62
        self._likelihood_model = None
63
64
        # These lists will contain the maps for the point sources
65
        self._convolved_point_sources = ConvolvedSourcesContainer()
66
        # and this one for extended sources
67
        self._convolved_ext_sources = ConvolvedSourcesContainer()
68
69
        # All energy/nHit bins are loaded in memory
70
        self._all_planes = list(self._maptree.analysis_bins_labels)
71
72
        # The active planes list always contains the list of *indexes* of the active planes
73
        self._active_planes = None
74
75
        # Set up the transformations from the flat-sky projection to Healpix, as well as the list of active pixels
76
        # (one for each energy/nHit bin). We make a separate transformation because different energy bins might have
77
        # different nsides
78
        self._active_pixels = collections.OrderedDict()
79
        self._flat_sky_to_healpix_transform = collections.OrderedDict()
80
81
        for bin_id in self._maptree:
82
83
            this_maptree = self._maptree[bin_id]
84
            this_nside = this_maptree.nside
85
            this_active_pixels = roi.active_pixels(this_nside)
86
87
            this_flat_sky_to_hpx_transform = FlatSkyToHealpixTransform(self._flat_sky_projection.wcs,
88
                                                                       'icrs',
89
                                                                       this_nside,
90
                                                                       this_active_pixels,
91
                                                                       (self._flat_sky_projection.npix_width,
92
                                                                        self._flat_sky_projection.npix_height),
93
                                                                       order='bilinear')
94
95
            self._active_pixels[bin_id] = this_active_pixels
96
            self._flat_sky_to_healpix_transform[bin_id] = this_flat_sky_to_hpx_transform
97
98
        # This will contain a list of PSF convolutors for extended sources, if there is any in the model
99
100
        self._psf_convolutors = None
101
102
        # Pre-compute the log-factorial factor in the likelihood, so we do not keep to computing it over and over
103
        # again.
104
        self._log_factorials = collections.OrderedDict()
105
106
        # We also apply a bias so that the numerical value of the log-likelihood stays small. This helps when
107
        # fitting with algorithms like MINUIT because the convergence criterium involves the difference between
108
        # two likelihood values, which would be affected by numerical precision errors if the two values are
109
        # too large
110
        self._saturated_model_like_per_maptree = collections.OrderedDict()
111
112
        # The actual computation is in a method so we can recall it on clone (see the get_simulated_dataset method)
113
        self._compute_likelihood_biases()
114
115
        # This will save a clone of self for simulations
116
        self._clone = None
117
118
    def _setup_psf_convolutors(self):
119
120
        central_response_bins = self._response.get_response_dec_bin(self._roi.ra_dec_center[1])
121
122
        self._psf_convolutors = collections.OrderedDict()
123
        for bin_id in central_response_bins:
124
            self._psf_convolutors[bin_id] = PSFConvolutor(central_response_bins[bin_id].psf, self._flat_sky_projection)
125
126
    def _compute_likelihood_biases(self):
127
128
        for bin_label in self._maptree:
129
130
            data_analysis_bin = self._maptree[bin_label]
131
132
            this_log_factorial = np.sum(logfactorial(data_analysis_bin.observation_map.as_partial()))
133
            self._log_factorials[bin_label] = this_log_factorial
134
135
            # As bias we use the likelihood value for the saturated model
136
            obs = data_analysis_bin.observation_map.as_partial()
137
            bkg = data_analysis_bin.background_map.as_partial()
138
139
            sat_model = np.clip(obs - bkg, 1e-50, None).astype(np.float64)
140
141
            self._saturated_model_like_per_maptree[bin_label] = log_likelihood(obs, bkg, sat_model) - this_log_factorial
142
143
    def get_saturated_model_likelihood(self):
144
        """
145
        Returns the likelihood for the saturated model (i.e. a model exactly equal to observation - background).
146
147
        :return:
148
        """
149
        return sum(self._saturated_model_like_per_maptree.values())
150
151
    def set_active_measurements(self, bin_id_min=None, bin_id_max=None, bin_list=None):
0 ignored issues
show
Coding Style introduced by
This method should have a docstring.

The coding style of this project requires that you add a docstring to this code element. Below, you find an example for methods:

class SomeClass:
    def some_method(self):
        """Do x and return foo."""

If you would like to know more about docstrings, we recommend to read PEP-257: Docstring Conventions.

Loading history...
152
153
        # Check for legal input
154
        if bin_id_min is not None:
155
156
            assert bin_id_max is not None, "If you provide a minimum bin, you also need to provide a maximum bin"
157
158
            # Make sure they are integers
159
            bin_id_min = int(bin_id_min)
160
            bin_id_max = int(bin_id_max)
161
162
            self._active_planes = []
163
            for this_bin in range(bin_id_min, bin_id_max + 1):
164
                this_bin = str(this_bin)
165
                if not this_bin in self._all_planes:
0 ignored issues
show
Unused Code introduced by
Consider changing "not this_bin in self._all_planes" to "this_bin not in self._all_planes"
Loading history...
166
167
                    raise ValueError("Bin %s it not contained in this response" % this_bin)
168
169
                self._active_planes.append(this_bin)
170
171
        else:
172
173
            assert bin_id_max is None, "If you provide a maximum bin, you also need to provide a minimum bin"
174
175
            assert bin_list is not None
176
177
            self._active_planes = []
178
179
            for this_bin in bin_list:
180
181
                if not this_bin in self._all_planes:
182
183
                    raise ValueError("Bin %s it not contained in this response" % this_bin)
184
185
                self._active_planes.append(this_bin)
186
187
    def display(self, verbose=False):
0 ignored issues
show
Coding Style introduced by
This method should have a docstring.

The coding style of this project requires that you add a docstring to this code element. Below, you find an example for methods:

class SomeClass:
    def some_method(self):
        """Do x and return foo."""

If you would like to know more about docstrings, we recommend to read PEP-257: Docstring Conventions.

Loading history...
188
189
        print("Region of Interest: ")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
190
        print("--------------------\n")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
191
        self._roi.display()
192
193
        print("")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
194
        print("Flat sky projection: ")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
195
        print("----------------------\n")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
196
197
        print("Width x height: %s x %s px" % (self._flat_sky_projection.npix_width,
198
                                              self._flat_sky_projection.npix_height))
199
        print("Pixel sizes: %s deg" % self._flat_sky_projection.pixel_size)
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
200
201
        print("")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
202
        print("Response: ")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
203
        print("----------\n")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
204
205
        self._response.display(verbose)
206
207
        print("")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
208
        print("Map Tree: ")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
209
        print("----------\n")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
210
211
        self._maptree.display()
212
213
        print("")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
214
        print("Active energy/nHit planes: ")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
215
        print("---------------------------\n")
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
216
        print(self._active_planes)
0 ignored issues
show
Unused Code Coding Style introduced by
There is an unnecessary parenthesis after print.
Loading history...
217
218
    def set_model(self, likelihood_model_instance):
219
        """
220
        Set the model to be used in the joint minimization. Must be a LikelihoodModel instance.
221
        """
222
223
        self._likelihood_model = likelihood_model_instance
224
225
        # Reset
226
        self._convolved_point_sources.reset()
227
        self._convolved_ext_sources.reset()
228
229
        # For each point source in the model, build the convolution class
230
231
        for source in self._likelihood_model.point_sources.values():
232
233
            this_convolved_point_source = ConvolvedPointSource(source, self._response, self._flat_sky_projection)
234
235
            self._convolved_point_sources.append(this_convolved_point_source)
236
237
        # Samewise for extended sources
238
239
        ext_sources = self._likelihood_model.extended_sources.values()
240
241
        if len(ext_sources) > 0:
0 ignored issues
show
Unused Code introduced by
Do not use len(SEQUENCE) to determine if a sequence is empty
Loading history...
242
243
            # We will need to convolve
244
245
            self._setup_psf_convolutors()
246
247
            for source in ext_sources:
248
249
                if source.spatial_shape.n_dim == 2:
250
251
                    this_convolved_ext_source = ConvolvedExtendedSource2D(source,
252
                                                                          self._response,
253
                                                                          self._flat_sky_projection)
254
255
                else:
256
257
                    this_convolved_ext_source = ConvolvedExtendedSource3D(source,
258
                                                                          self._response,
259
                                                                          self._flat_sky_projection)
260
261
                self._convolved_ext_sources.append(this_convolved_ext_source)
262
263
    def display_spectrum(self):
0 ignored issues
show
Coding Style introduced by
This method should have a docstring.

The coding style of this project requires that you add a docstring to this code element. Below, you find an example for methods:

class SomeClass:
    def some_method(self):
        """Do x and return foo."""

If you would like to know more about docstrings, we recommend to read PEP-257: Docstring Conventions.

Loading history...
Comprehensibility introduced by
This function exceeds the maximum number of variables (26/15).
Loading history...
264
265
        n_point_sources = self._likelihood_model.get_number_of_point_sources()
266
        n_ext_sources = self._likelihood_model.get_number_of_extended_sources()
267
268
        total_counts = np.zeros(len(self._active_planes), dtype=float)
269
        total_model = np.zeros_like(total_counts)
270
        model_only = np.zeros_like(total_counts)
271
        net_counts = np.zeros_like(total_counts)
272
        yerr_low = np.zeros_like(total_counts)
273
        yerr_high = np.zeros_like(total_counts)
274
275
        for i, energy_id in enumerate(self._active_planes):
276
277
            data_analysis_bin = self._maptree[energy_id]
278
279
            this_model_map_hpx = self._get_expectation(data_analysis_bin, energy_id, n_point_sources, n_ext_sources)
280
281
            this_model_tot = np.sum(this_model_map_hpx)
282
283
            this_data_tot = np.sum(data_analysis_bin.observation_map.as_partial())
284
            this_bkg_tot = np.sum(data_analysis_bin.background_map.as_partial())
285
286
            total_counts[i] = this_data_tot
287
            net_counts[i] = this_data_tot - this_bkg_tot
288
            model_only[i] = this_model_tot
289
290
            this_wh_model = this_model_tot + this_bkg_tot
291
            total_model[i] = this_wh_model
292
293
            if this_data_tot >= 50.0:
294
295
                # Gaussian limit
296
                # Under the null hypothesis the data are distributed as a Gaussian with mu = model
297
                # and sigma = sqrt(model)
298
                # NOTE: since we neglect the background uncertainty, the background is part of the
299
                # model
300
                yerr_low[i] = np.sqrt(this_data_tot)
301
                yerr_high[i] = np.sqrt(this_data_tot)
302
303
            else:
304
305
                # Low-counts
306
                # Under the null hypothesis the data are distributed as a Poisson distribution with
307
                # mean = model, plot the 68% confidence interval (quantile=[0.16,1-0.16]).
308
                # NOTE: since we neglect the background uncertainty, the background is part of the
309
                # model
310
                quantile = 0.16
311
                mean = this_wh_model
312
                y_low = poisson.isf(1-quantile, mu=mean)
313
                y_high = poisson.isf(quantile, mu=mean)
314
                yerr_low[i] = mean-y_low
315
                yerr_high[i] = y_high-mean
316
317
        residuals = (total_counts - total_model) / np.sqrt(total_model)
318
        residuals_err = [yerr_high / np.sqrt(total_model),
319
                         yerr_low / np.sqrt(total_model)]
320
321
        fig, subs = plt.subplots(2, 1, gridspec_kw={'height_ratios': [2, 1], 'hspace': 0})
322
323
        subs[0].errorbar(self._active_planes, net_counts, yerr=[yerr_high, yerr_low],
324
                         capsize=0,
325
                         color='black', label='Net counts', fmt='.')
326
327
        subs[0].plot(self._active_planes, model_only, label='Convolved model')
328
329
        subs[0].legend(bbox_to_anchor=(1.0, 1.0), loc="upper right",
330
                       numpoints=1)
331
332
        # Residuals
333
        subs[1].axhline(0, linestyle='--')
334
335
        subs[1].errorbar(
336
            self._active_planes, residuals,
337
            yerr=residuals_err,
338
            capsize=0, fmt='.'
339
        )
340
341
        y_limits = [min(net_counts[net_counts > 0]) / 2., max(net_counts) * 2.]
342
343
        subs[0].set_yscale("log", nonposy='clip')
344
        subs[0].set_ylabel("Counts per bin")
345
        subs[0].set_xticks([])
346
347
        subs[1].set_xlabel("Analysis bin")
348
        subs[1].set_ylabel(r"$\frac{{cts - mod - bkg}}{\sqrt{mod + bkg}}$")
349
        subs[1].set_xticks(self._active_planes)
350
        subs[1].set_xticklabels(self._active_planes)
351
352
        subs[0].set_ylim(y_limits)
353
354
        return fig
355
356
    def get_log_like(self):
357
        """
358
        Return the value of the log-likelihood with the current values for the
359
        parameters
360
        """
361
362
        n_point_sources = self._likelihood_model.get_number_of_point_sources()
363
        n_ext_sources = self._likelihood_model.get_number_of_extended_sources()
364
365
        # Make sure that no source has been added since we filled the cache
366
        assert n_point_sources == self._convolved_point_sources.n_sources_in_cache and \
367
               n_ext_sources == self._convolved_ext_sources.n_sources_in_cache, \
368
            "The number of sources has changed. Please re-assign the model to the plugin."
369
370
        # This will hold the total log-likelihood
371
        total_log_like = 0
372
373
        for bin_id in self._active_planes:
374
375
            data_analysis_bin = self._maptree[bin_id]
376
377
            this_model_map_hpx = self._get_expectation(data_analysis_bin, bin_id, n_point_sources, n_ext_sources)
378
379
            # Now compare with observation
380
            bkg_renorm = self._nuisance_parameters.values()[0].value
381
382
            obs = data_analysis_bin.observation_map.as_partial()  # type: np.array
383
            bkg = data_analysis_bin.background_map.as_partial() * bkg_renorm  # type: np.array
384
385
            this_pseudo_log_like = log_likelihood(obs,
386
                                                  bkg,
387
                                                  this_model_map_hpx)
388
389
            total_log_like += this_pseudo_log_like - self._log_factorials[bin_id] - self._saturated_model_like_per_maptree[bin_id]
0 ignored issues
show
Coding Style introduced by
This line is too long as per the coding-style (130/120).

This check looks for lines that are too long. You can specify the maximum line length.

Loading history...
390
391
        return total_log_like
392
393
    def write(self, response_file_name, map_tree_file_name):
394
        """
395
        Write this dataset to disk in HDF format.
396
397
        :param response_file_name: filename for the response
398
        :param map_tree_file_name: filename for the map tree
399
        :return: None
400
        """
401
402
        self._maptree.write(map_tree_file_name)
403
        self._response.write(response_file_name)
404
405
    def get_simulated_dataset(self, name):
0 ignored issues
show
Coding Style introduced by
This method should have a docstring.

The coding style of this project requires that you add a docstring to this code element. Below, you find an example for methods:

class SomeClass:
    def some_method(self):
        """Do x and return foo."""

If you would like to know more about docstrings, we recommend to read PEP-257: Docstring Conventions.

Loading history...
406
407
        # First get expectation under the current model and store them, if we didn't do it yet
408
409
        if self._clone is None:
410
411
            n_point_sources = self._likelihood_model.get_number_of_point_sources()
412
            n_ext_sources = self._likelihood_model.get_number_of_extended_sources()
413
414
            expectations = []
415
416
            for i, data_analysis_bin in enumerate(self._maptree):
417
418
                if i not in self._active_planes:
419
420
                    expectations.append(None)
421
422
                else:
423
424
                    expectations.append(self._get_expectation(data_analysis_bin, i,
425
                                                              n_point_sources, n_ext_sources) +
426
                                        data_analysis_bin.background_map.as_partial())
427
428
            if parallel_client.is_parallel_computation_active():
429
430
                # Do not clone, as the parallel environment already makes clones
431
432
                clone = self
433
434
            else:
435
436
                clone = copy.deepcopy(self)
437
438
            self._clone = (clone, expectations)
439
440
        # Substitute the observation and background for each data analysis bin
441
        for i, (data_analysis_bin, orig_data_analysis_bin) in enumerate(zip(self._clone[0]._maptree, self._maptree)):
0 ignored issues
show
Coding Style Best Practice introduced by
It seems like _maptree was declared protected and should not be accessed from this context.

Prefixing a member variable _ is usually regarded as the equivalent of declaring it with protected visibility that exists in other languages. Consequentially, such a member should only be accessed from the same class or a child class:

class MyParent:
    def __init__(self):
        self._x = 1;
        self.y = 2;

class MyChild(MyParent):
    def some_method(self):
        return self._x    # Ok, since accessed from a child class

class AnotherClass:
    def some_method(self, instance_of_my_child):
        return instance_of_my_child._x   # Would be flagged as AnotherClass is not
                                         # a child class of MyParent
Loading history...
Unused Code introduced by
The variable orig_data_analysis_bin seems to be unused.
Loading history...
442
443
            if i not in self._active_planes:
444
445
                continue
446
447
            else:
448
449
                # Active plane. Generate new data
450
                expectation = self._clone[1][i]
451
                new_data = np.random.poisson(expectation, size=(1, expectation.shape[0])).flatten()
452
453
                # Substitute data
454
                data_analysis_bin.observation_map.set_new_values(new_data)
455
456
        # Now change name and return
457
        self._clone[0]._name = name
0 ignored issues
show
Coding Style Best Practice introduced by
It seems like _name was declared protected and should not be accessed from this context.

Prefixing a member variable _ is usually regarded as the equivalent of declaring it with protected visibility that exists in other languages. Consequentially, such a member should only be accessed from the same class or a child class:

class MyParent:
    def __init__(self):
        self._x = 1;
        self.y = 2;

class MyChild(MyParent):
    def some_method(self):
        return self._x    # Ok, since accessed from a child class

class AnotherClass:
    def some_method(self, instance_of_my_child):
        return instance_of_my_child._x   # Would be flagged as AnotherClass is not
                                         # a child class of MyParent
Loading history...
458
        # Adjust the name of the nuisance parameter
459
        old_name = self._clone[0]._nuisance_parameters.keys()[0]
0 ignored issues
show
Coding Style Best Practice introduced by
It seems like _nuisance_parameters was declared protected and should not be accessed from this context.

Prefixing a member variable _ is usually regarded as the equivalent of declaring it with protected visibility that exists in other languages. Consequentially, such a member should only be accessed from the same class or a child class:

class MyParent:
    def __init__(self):
        self._x = 1;
        self.y = 2;

class MyChild(MyParent):
    def some_method(self):
        return self._x    # Ok, since accessed from a child class

class AnotherClass:
    def some_method(self, instance_of_my_child):
        return instance_of_my_child._x   # Would be flagged as AnotherClass is not
                                         # a child class of MyParent
Loading history...
460
        new_name = old_name.replace(self.name, name)
461
        self._clone[0]._nuisance_parameters[new_name] = self._clone[0]._nuisance_parameters.pop(old_name)
0 ignored issues
show
Coding Style Best Practice introduced by
It seems like _nuisance_parameters was declared protected and should not be accessed from this context.

Prefixing a member variable _ is usually regarded as the equivalent of declaring it with protected visibility that exists in other languages. Consequentially, such a member should only be accessed from the same class or a child class:

class MyParent:
    def __init__(self):
        self._x = 1;
        self.y = 2;

class MyChild(MyParent):
    def some_method(self):
        return self._x    # Ok, since accessed from a child class

class AnotherClass:
    def some_method(self, instance_of_my_child):
        return instance_of_my_child._x   # Would be flagged as AnotherClass is not
                                         # a child class of MyParent
Loading history...
462
463
        # Recompute biases
464
        self._clone[0]._compute_likelihood_biases()
0 ignored issues
show
Coding Style Best Practice introduced by
It seems like _compute_likelihood_biases was declared protected and should not be accessed from this context.

Prefixing a member variable _ is usually regarded as the equivalent of declaring it with protected visibility that exists in other languages. Consequentially, such a member should only be accessed from the same class or a child class:

class MyParent:
    def __init__(self):
        self._x = 1;
        self.y = 2;

class MyChild(MyParent):
    def some_method(self):
        return self._x    # Ok, since accessed from a child class

class AnotherClass:
    def some_method(self, instance_of_my_child):
        return instance_of_my_child._x   # Would be flagged as AnotherClass is not
                                         # a child class of MyParent
Loading history...
465
466
        return self._clone[0]
467
468
    def _get_expectation(self, data_analysis_bin, energy_bin_id, n_point_sources, n_ext_sources):
469
470
        # Compute the expectation from the model
471
472
        this_model_map = None
473
474
        for pts_id in range(n_point_sources):
475
476
            this_convolved_source = self._convolved_point_sources[pts_id]
477
478
            expectation_per_transit = this_convolved_source.get_source_map(energy_bin_id, tag=None)
479
480
            expectation_from_this_source = expectation_per_transit * data_analysis_bin.n_transits
481
482
            if this_model_map is None:
483
484
                # First addition
485
486
                this_model_map = expectation_from_this_source
487
488
            else:
489
490
                this_model_map += expectation_from_this_source
491
492
        # Now process extended sources
493
        if n_ext_sources > 0:
494
495
            this_ext_model_map = None
496
497
            for ext_id in range(n_ext_sources):
498
499
                this_convolved_source = self._convolved_ext_sources[ext_id]
500
501
                expectation_per_transit = this_convolved_source.get_source_map(energy_bin_id)
502
503
                if this_ext_model_map is None:
504
505
                    # First addition
506
507
                    this_ext_model_map = expectation_per_transit
508
509
                else:
510
511
                    this_ext_model_map += expectation_per_transit
512
513
            # Now convolve with the PSF
514
            if this_model_map is None:
515
                
0 ignored issues
show
Coding Style introduced by
Trailing whitespace
Loading history...
516
                # Only extended sources
517
            
0 ignored issues
show
Coding Style introduced by
Trailing whitespace
Loading history...
518
                this_model_map = (self._psf_convolutors[energy_bin_id].extended_source_image(this_ext_model_map) *
519
                                  data_analysis_bin.n_transits)
520
            
0 ignored issues
show
Coding Style introduced by
Trailing whitespace
Loading history...
521
            else:
522
523
                this_model_map += (self._psf_convolutors[energy_bin_id].extended_source_image(this_ext_model_map) *
524
                                   data_analysis_bin.n_transits)
525
526
527
        # Now transform from the flat sky projection to HEALPiX
528
529
        if this_model_map is not None:
530
531
            # First divide for the pixel area because we need to interpolate brightness
532
            this_model_map = this_model_map / self._flat_sky_projection.project_plane_pixel_area
533
534
            this_model_map_hpx = self._flat_sky_to_healpix_transform[energy_bin_id](this_model_map, fill_value=0.0)
535
536
            # Now multiply by the pixel area of the new map to go back to flux
537
            this_model_map_hpx *= hp.nside2pixarea(data_analysis_bin.nside, degrees=True)
538
539
        else:
540
541
            # No sources
542
543
            this_model_map_hpx = 0.0
544
545
        return this_model_map_hpx
546
547
    @staticmethod
548
    def _represent_healpix_map(fig, hpx_map, longitude, latitude, xsize, resolution, smoothing_kernel_sigma):
0 ignored issues
show
best-practice introduced by
Too many arguments (7/5)
Loading history...
549
550
        proj = get_gnomonic_projection(fig, hpx_map,
551
                                       rot=(longitude, latitude, 0.0),
552
                                       xsize=xsize,
553
                                       reso=resolution)
554
555
        if smoothing_kernel_sigma is not None:
556
557
            # Get the sigma in pixels
558
            sigma = smoothing_kernel_sigma * 60 / resolution
559
560
            proj = convolve(list(proj),
561
                            Gaussian2DKernel(sigma),
562
                            nan_treatment='fill',
563
                            preserve_nan=True)
564
565
        return proj
566
567
    def display_fit(self, smoothing_kernel_sigma=0.1):
0 ignored issues
show
Coding Style introduced by
This method should have a docstring.

The coding style of this project requires that you add a docstring to this code element. Below, you find an example for methods:

class SomeClass:
    def some_method(self):
        """Do x and return foo."""

If you would like to know more about docstrings, we recommend to read PEP-257: Docstring Conventions.

Loading history...
Comprehensibility introduced by
This function exceeds the maximum number of variables (24/15).
Loading history...
568
569
        n_point_sources = self._likelihood_model.get_number_of_point_sources()
570
        n_ext_sources = self._likelihood_model.get_number_of_extended_sources()
571
572
        # This is the resolution (i.e., the size of one pixel) of the image
573
        resolution = 3.0 # arcmin
574
575
        # The image is going to cover the diameter plus 20% padding
576
        xsize = self._get_optimal_xsize(resolution)
577
578
        n_active_planes = len(self._active_planes)
579
580
        fig, subs = plt.subplots(n_active_planes, 3, figsize=(8, n_active_planes * 2))
581
582
        with progress_bar(len(self._active_planes), title='Smoothing maps') as prog_bar:
583
584
            for i, plane_id in enumerate(self._active_planes):
585
586
                data_analysis_bin = self._maptree[plane_id]
587
588
                # Get the center of the projection for this plane
589
                this_ra, this_dec = self._roi.ra_dec_center
590
591
                this_model_map_hpx = self._get_expectation(data_analysis_bin, plane_id, n_point_sources, n_ext_sources)
592
593
                # Make a full healpix map for a second
594
                whole_map = SparseHealpix(this_model_map_hpx,
595
                                          self._active_pixels[plane_id],
596
                                          data_analysis_bin.observation_map.nside).as_dense()
597
598
                # Healpix uses longitude between -180 and 180, while R.A. is between 0 and 360. We need to fix that:
599
                longitude = ra_to_longitude(this_ra)
600
601
                # Declination is already between -90 and 90
602
                latitude = this_dec
603
604
                # Plot model
605
606
                proj_m = self._represent_healpix_map(fig, whole_map,
607
                                                     longitude, latitude,
608
                                                     xsize, resolution, smoothing_kernel_sigma)
609
610
                subs[i][0].imshow(proj_m, origin='lower')
611
612
                # Remove numbers from axis
613
                subs[i][0].axis('off')
614
615
                # Plot data map
616
                # Here we removed the background otherwise nothing is visible
617
                # Get background (which is in a way "part of the model" since the uncertainties are neglected)
618
                background_map = data_analysis_bin.background_map.as_dense()
619
                bkg_subtracted = data_analysis_bin.observation_map.as_dense() - background_map
620
621
                proj_d = self._represent_healpix_map(fig, bkg_subtracted,
622
                                                     longitude, latitude,
623
                                                     xsize, resolution, smoothing_kernel_sigma)
624
625
                subs[i][1].imshow(proj_d, origin='lower')
626
627
                # Remove numbers from axis
628
                subs[i][1].axis('off')
629
630
                # Now residuals
631
                res = proj_d - proj_m
632
                # proj_res = self._represent_healpix_map(fig, res,
633
                #                                        longitude, latitude,
634
                #                                        xsize, resolution, smoothing_kernel_sigma)
635
                subs[i][2].imshow(res, origin='lower')
636
637
                # Remove numbers from axis
638
                subs[i][2].axis('off')
639
640
                subs[i][0].set_title('model, bin {}'.format(data_analysis_bin.name))
641
                subs[i][1].set_title('excess, bin {}'.format(data_analysis_bin.name))
642
                subs[i][2].set_title('residuals, bin {}'.format(data_analysis_bin.name))
643
644
                prog_bar.increase()
645
646
        fig.set_tight_layout(True)
647
648
        return fig
649
650
    def _get_optimal_xsize(self, resolution):
651
652
        return 2.2 * self._roi.data_radius.to("deg").value / (resolution / 60.0)
653
654
    def display_stacked_image(self, smoothing_kernel_sigma=0.5):
0 ignored issues
show
Coding Style introduced by
This method should have a docstring.

The coding style of this project requires that you add a docstring to this code element. Below, you find an example for methods:

class SomeClass:
    def some_method(self):
        """Do x and return foo."""

If you would like to know more about docstrings, we recommend to read PEP-257: Docstring Conventions.

Loading history...
Comprehensibility introduced by
This function exceeds the maximum number of variables (19/15).
Loading history...
655
656
        # This is the resolution (i.e., the size of one pixel) of the image in arcmin
657
        resolution = 3.0
658
659
        # The image is going to cover the diameter plus 20% padding
660
        xsize = self._get_optimal_xsize(resolution)
661
662
        active_planes_bins = map(lambda x: self._maptree[x], self._active_planes)
0 ignored issues
show
introduced by
map/filter on lambda could be replaced by comprehension
Loading history...
663
664
        # Get the center of the projection for this plane
665
        this_ra, this_dec = self._roi.ra_dec_center
666
667
        # Healpix uses longitude between -180 and 180, while R.A. is between 0 and 360. We need to fix that:
668
        longitude = ra_to_longitude(this_ra)
669
670
        # Declination is already between -90 and 90
671
        latitude = this_dec
672
673
        total = None
674
675
        for i, data_analysis_bin in enumerate(active_planes_bins):
676
677
            # Plot data
678
            background_map = data_analysis_bin.background_map.as_dense()
679
            this_data = data_analysis_bin.observation_map.as_dense() - background_map
680
            idx = np.isnan(this_data)
681
            # this_data[idx] = hp.UNSEEN
682
683
            if i == 0:
684
685
                total = this_data
686
687
            else:
688
689
                # Sum only when there is no UNSEEN, so that the UNSEEN pixels will stay UNSEEN
690
                total[~idx] += this_data[~idx]
691
692
        delta_coord = (self._roi.data_radius.to("deg").value * 2.0) / 15.0
693
694
        fig, sub = plt.subplots(1,1)
0 ignored issues
show
Coding Style introduced by
Exactly one space required after comma
Loading history...
695
696
        proj = self._represent_healpix_map(fig, total, longitude, latitude, xsize, resolution, smoothing_kernel_sigma)
697
698
        sub.imshow(proj, origin='lower')
699
        sub.axis('off')
700
701
        hp.graticule(delta_coord, delta_coord)
702
703
        return fig
704
705
    def inner_fit(self):
706
        """
707
        This is used for the profile likelihood. Keeping fixed all parameters in the
708
        LikelihoodModel, this method minimize the logLike over the remaining nuisance
709
        parameters, i.e., the parameters belonging only to the model for this
710
        particular detector. If there are no nuisance parameters, simply return the
711
        logLike value.
712
        """
713
714
        return self.get_log_like()
715
716
    def get_number_of_data_points(self):
0 ignored issues
show
Coding Style introduced by
This method should have a docstring.

The coding style of this project requires that you add a docstring to this code element. Below, you find an example for methods:

class SomeClass:
    def some_method(self):
        """Do x and return foo."""

If you would like to know more about docstrings, we recommend to read PEP-257: Docstring Conventions.

Loading history...
717
718
        n_points = 0
719
720
        for bin_id in self._maptree:
721
            n_points += self._maptree[bin_id].observation_map.as_partial().shape[0]
722
723
        return n_points
724