gl_space.py 13.5 KB
Newer Older
csongor's avatar
csongor committed
1
2
from __future__ import division

Jait Dixit's avatar
Jait Dixit committed
3
import itertools
csongor's avatar
csongor committed
4
5
6
7
8
9
10
import numpy as np
import pylab as pl
from matplotlib.colors import LogNorm as ln
from matplotlib.ticker import LogFormatter as lf

from d2o import STRATEGIES as DISTRIBUTION_STRATEGIES

11
from nifty.spaces.lm_space import LMSpace
12

13
from nifty.spaces.space import Space
14
15
from nifty.config import about, nifty_configuration as gc,\
                         dependency_injector as gdi
theos's avatar
theos committed
16
from gl_space_paradict import GLSpaceParadict
csongor's avatar
csongor committed
17
18
19
20
21
22
from nifty.nifty_random import random

gl = gdi.get('libsharp_wrapper_gl')

GL_DISTRIBUTION_STRATEGIES = DISTRIBUTION_STRATEGIES['global']

23
24

class GLSpace(Space):
csongor's avatar
csongor committed
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
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
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
    """
        ..                 __
        ..               /  /
        ..     ____ __  /  /
        ..   /   _   / /  /
        ..  /  /_/  / /  /_
        ..  \___   /  \___/  space class
        .. /______/

        NIFTY subclass for Gauss-Legendre pixelizations [#]_ of the two-sphere.

        Parameters
        ----------
        nlat : int
            Number of latitudinal bins, or rings.
        nlon : int, *optional*
            Number of longitudinal bins (default: ``2*nlat - 1``).
        dtype : numpy.dtype, *optional*
            Data type of the field values (default: numpy.float64).

        See Also
        --------
        hp_space : A class for the HEALPix discretization of the sphere [#]_.
        lm_space : A class for spherical harmonic components.

        Notes
        -----
        Only real-valued fields on the two-sphere are supported, i.e.
        `dtype` has to be either numpy.float64 or numpy.float32.

        References
        ----------
        .. [#] M. Reinecke and D. Sverre Seljebotn, 2013, "Libsharp - spherical
               harmonic transforms revisited";
               `arXiv:1303.4945 <http://www.arxiv.org/abs/1303.4945>`_
        .. [#] K.M. Gorski et al., 2005, "HEALPix: A Framework for
               High-Resolution Discretization and Fast Analysis of Data
               Distributed on the Sphere", *ApJ* 622..759G.

        Attributes
        ----------
        para : numpy.ndarray
            One-dimensional array containing the two numbers `nlat` and `nlon`.
        dtype : numpy.dtype
            Data type of the field values.
        discrete : bool
            Whether or not the underlying space is discrete, always ``False``
            for spherical spaces.
        vol : numpy.ndarray
            An array containing the pixel sizes.
    """

    def __init__(self, nlat, nlon=None, dtype=np.dtype('float64')):
        """
            Sets the attributes for a gl_space class instance.

            Parameters
            ----------
            nlat : int
                Number of latitudinal bins, or rings.
            nlon : int, *optional*
                Number of longitudinal bins (default: ``2*nlat - 1``).
            dtype : numpy.dtype, *optional*
                Data type of the field values (default: numpy.float64).

            Returns
            -------
            None

            Raises
            ------
            ImportError
                If the libsharp_wrapper_gl module is not available.
            ValueError
                If input `nlat` is invaild.

        """
        # check imports
        if not gc['use_libsharp']:
            raise ImportError(about._errors.cstring(
                "ERROR: libsharp_wrapper_gl not loaded."))

Jait Dixit's avatar
Jait Dixit committed
107
        # setup paradict
theos's avatar
theos committed
108
        self.paradict = GLSpaceParadict(nlat=nlat, nlon=nlon)
csongor's avatar
csongor committed
109

Jait Dixit's avatar
Jait Dixit committed
110
        # check and set data type
csongor's avatar
csongor committed
111
112
113
114
115
116
        dtype = np.dtype(dtype)
        if dtype not in [np.dtype('float32'), np.dtype('float64')]:
            about.warnings.cprint("WARNING: data type set to default.")
            dtype = np.dtype('float')
        self.dtype = dtype

Jait Dixit's avatar
Jait Dixit committed
117
118
        # GLSpace is not harmonic
        self._harmonic = False
csongor's avatar
csongor committed
119
120

    def copy(self):
121
122
123
        return GLSpace(nlat=self.paradict['nlat'],
                       nlon=self.paradict['nlon'],
                       dtype=self.dtype)
csongor's avatar
csongor committed
124
125
126
127
128

    @property
    def shape(self):
        return (np.int((self.paradict['nlat'] * self.paradict['nlon'])),)

129
130
131
132
    @property
    def vol(self):
        return np.sum(self.paradict['nlon'] * np.array(self.distances[0]))

Jait Dixit's avatar
Jait Dixit committed
133
    def weight(self, x, power=1, axes=None, inplace=False):
134
135
136
137
138
        # check if the axes provided are valid given the input shape
        if axes is not None and \
                not all(axis in range(len(x.shape)) for axis in axes):
            raise ValueError("ERROR: Provided axes does not match array shape")

Jait Dixit's avatar
Jait Dixit committed
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
        weight = np.array(list(
            itertools.chain.from_iterable(
                itertools.repeat(x ** power, self.paradict['nlon'])
                for x in gl.vol(self.paradict['nlat'])
            )
        ))

        if axes is not None:
            # reshape the weight array to match the input shape
            new_shape = np.ones(x.shape)
            for index in range(len(axes)):
                new_shape[index] = len(weight)
            weight = weight.reshape(new_shape)

        if inplace:
            x *= weight
            result_x = x
csongor's avatar
csongor committed
156
        else:
Jait Dixit's avatar
Jait Dixit committed
157
            result_x = x * weight
csongor's avatar
csongor committed
158

Jait Dixit's avatar
Jait Dixit committed
159
        return result_x
csongor's avatar
csongor committed
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
207

    def get_plot(self, x, title="", vmin=None, vmax=None, power=False,
                 unit="", norm=None, cmap=None, cbar=True, other=None,
                 legend=False, mono=True, **kwargs):
        """
            Creates a plot of field values according to the specifications
            given by the parameters.

            Parameters
            ----------
            x : numpy.ndarray
                Array containing the field values.

            Returns
            -------
            None

            Other parameters
            ----------------
            title : string, *optional*
                Title of the plot (default: "").
            vmin : float, *optional*
                Minimum value to be displayed (default: ``min(x)``).
            vmax : float, *optional*
                Maximum value to be displayed (default: ``max(x)``).
            power : bool, *optional*
                Whether to plot the power contained in the field or the field
                values themselves (default: False).
            unit : string, *optional*
                Unit of the field values (default: "").
            norm : string, *optional*
                Scaling of the field values before plotting (default: None).
            cmap : matplotlib.colors.LinearSegmentedColormap, *optional*
                Color map to be used for two-dimensional plots (default: None).
            cbar : bool, *optional*
                Whether to show the color bar or not (default: True).
            other : {single object, tuple of objects}, *optional*
                Object or tuple of objects to be added, where objects can be
                scalars, arrays, or fields (default: None).
            legend : bool, *optional*
                Whether to show the legend or not (default: False).
            mono : bool, *optional*
                Whether to plot the monopole or not (default: True).
            save : string, *optional*
                Valid file name where the figure is to be stored, by default
                the figure is not saved (default: False).

        """
208
209
        from nifty.field import Field

csongor's avatar
csongor committed
210
211
212
213
214
        try:
            x = x.get_full_data()
        except AttributeError:
            pass

Jait Dixit's avatar
Jait Dixit committed
215
        if (not pl.isinteractive()) and (not bool(kwargs.get("save", False))):
csongor's avatar
csongor committed
216
217
            about.warnings.cprint("WARNING: interactive mode off.")

Jait Dixit's avatar
Jait Dixit committed
218
        if (power):
csongor's avatar
csongor committed
219
220
            x = self.calc_power(x)

Jait Dixit's avatar
Jait Dixit committed
221
222
223
224
            fig = pl.figure(num=None, figsize=(6.4, 4.8), dpi=None,
                            facecolor="none",
                            edgecolor="none", frameon=False,
                            FigureClass=pl.Figure)
csongor's avatar
csongor committed
225
226
227
            ax0 = fig.add_axes([0.12, 0.12, 0.82, 0.76])

            xaxes = np.arange(self.para[0], dtype=np.int)
Jait Dixit's avatar
Jait Dixit committed
228
            if (vmin is None):
csongor's avatar
csongor committed
229
                vmin = np.min(x[:mono].tolist(
Jait Dixit's avatar
Jait Dixit committed
230
231
232
                ) + (xaxes * (2 * xaxes + 1) * x)[1:].tolist(), axis=None,
                              out=None)
            if (vmax is None):
csongor's avatar
csongor committed
233
                vmax = np.max(x[:mono].tolist(
Jait Dixit's avatar
Jait Dixit committed
234
235
                ) + (xaxes * (2 * xaxes + 1) * x)[1:].tolist(), axis=None,
                              out=None)
csongor's avatar
csongor committed
236
            ax0.loglog(xaxes[1:], (xaxes * (2 * xaxes + 1) * x)[1:], color=[0.0,
Jait Dixit's avatar
Jait Dixit committed
237
238
239
240
241
242
243
244
245
246
247
                                                                            0.5,
                                                                            0.0],
                       label="graph 0", linestyle='-', linewidth=2.0, zorder=1)
            if (mono):
                ax0.scatter(0.5 * (xaxes[1] + xaxes[2]), x[0], s=20,
                            color=[0.0, 0.5, 0.0], marker='o',
                            cmap=None, norm=None, vmin=None, vmax=None,
                            alpha=None, linewidths=None, verts=None, zorder=1)

            if (other is not None):
                if (isinstance(other, tuple)):
csongor's avatar
csongor committed
248
249
                    other = list(other)
                    for ii in xrange(len(other)):
Jait Dixit's avatar
Jait Dixit committed
250
                        if (isinstance(other[ii], Field)):
csongor's avatar
csongor committed
251
252
253
                            other[ii] = other[ii].power(**kwargs)
                        else:
                            other[ii] = self.enforce_power(other[ii])
Jait Dixit's avatar
Jait Dixit committed
254
                elif (isinstance(other, Field)):
csongor's avatar
csongor committed
255
256
257
258
259
                    other = [other.power(**kwargs)]
                else:
                    other = [self.enforce_power(other)]
                imax = max(1, len(other) - 1)
                for ii in xrange(len(other)):
Jait Dixit's avatar
Jait Dixit committed
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
                    ax0.loglog(xaxes[1:],
                               (xaxes * (2 * xaxes + 1) * other[ii])[1:],
                               color=[max(0.0, 1.0 - (2 * ii / imax) ** 2),
                                      0.5 * ((2 * ii - imax) / imax)
                                      ** 2, max(0.0, 1.0 - (
                                   2 * (ii - imax) / imax) ** 2)],
                               label="graph " + str(ii + 1), linestyle='-',
                               linewidth=1.0, zorder=-ii)
                    if (mono):
                        ax0.scatter(0.5 * (xaxes[1] + xaxes[2]), other[ii][0],
                                    s=20,
                                    color=[max(0.0, 1.0 - (2 * ii / imax) ** 2),
                                           0.5 * ((2 * ii - imax) / imax) ** 2,
                                           max(
                                               0.0, 1.0 - (
                                               2 * (ii - imax) / imax) ** 2)],
                                    marker='o', cmap=None, norm=None, vmin=None,
                                    vmax=None, alpha=None, linewidths=None,
                                    verts=None, zorder=-ii)
                if (legend):
csongor's avatar
csongor committed
280
281
282
283
284
285
286
287
288
289
                    ax0.legend()

            ax0.set_xlim(xaxes[1], xaxes[-1])
            ax0.set_xlabel(r"$l$")
            ax0.set_ylim(vmin, vmax)
            ax0.set_ylabel(r"$l(2l+1) C_l$")
            ax0.set_title(title)

        else:
            x = self.cast(x)
Jait Dixit's avatar
Jait Dixit committed
290
            if (vmin is None):
csongor's avatar
csongor committed
291
                vmin = np.min(x, axis=None, out=None)
Jait Dixit's avatar
Jait Dixit committed
292
            if (vmax is None):
csongor's avatar
csongor committed
293
                vmax = np.max(x, axis=None, out=None)
Jait Dixit's avatar
Jait Dixit committed
294
            if (norm == "log") and (vmin <= 0):
csongor's avatar
csongor committed
295
296
297
                raise ValueError(about._errors.cstring(
                    "ERROR: nonpositive value(s)."))

Jait Dixit's avatar
Jait Dixit committed
298
299
300
301
            fig = pl.figure(num=None, figsize=(8.5, 5.4), dpi=None,
                            facecolor="none",
                            edgecolor="none", frameon=False,
                            FigureClass=pl.Figure)
csongor's avatar
csongor committed
302
303
304
305
306
            ax0 = fig.add_axes([0.02, 0.05, 0.96, 0.9])

            lon, lat = gl.bounds(self.para[0], nlon=self.para[1])
            lon = (lon - np.pi) * 180 / np.pi
            lat = (lat - np.pi / 2) * 180 / np.pi
Jait Dixit's avatar
Jait Dixit committed
307
            if (norm == "log"):
csongor's avatar
csongor committed
308
309
310
                n_ = ln(vmin=vmin, vmax=vmax)
            else:
                n_ = None
Jait Dixit's avatar
Jait Dixit committed
311
312
313
314
315
            sub = ax0.pcolormesh(lon, lat, np.roll(
                np.array(x).reshape((self.para[0], self.para[1]), order='C'),
                self.para[
                    1] // 2, axis=1)[::-1, ::-1], cmap=cmap, norm=n_, vmin=vmin,
                                 vmax=vmax)
csongor's avatar
csongor committed
316
317
318
319
            ax0.set_xlim(-180, 180)
            ax0.set_ylim(-90, 90)
            ax0.set_aspect("equal")
            ax0.axis("off")
Jait Dixit's avatar
Jait Dixit committed
320
321
            if (cbar):
                if (norm == "log"):
csongor's avatar
csongor committed
322
323
324
325
326
327
328
                    f_ = lf(10, labelOnlyBase=False)
                    b_ = sub.norm.inverse(np.linspace(0, 1, sub.cmap.N + 1))
                    v_ = np.linspace(sub.norm.vmin, sub.norm.vmax, sub.cmap.N)
                else:
                    f_ = None
                    b_ = None
                    v_ = None
Jait Dixit's avatar
Jait Dixit committed
329
330
331
332
333
334
335
336
337
                cb0 = fig.colorbar(sub, ax=ax0, orientation="horizontal",
                                   fraction=0.1, pad=0.05, shrink=0.5,
                                   aspect=25, ticks=[
                        vmin, vmax], format=f_, drawedges=False, boundaries=b_,
                                   values=v_)
                cb0.ax.text(0.5, -1.0, unit, fontdict=None, withdash=False,
                            transform=cb0.ax.transAxes,
                            horizontalalignment="center",
                            verticalalignment="center")
csongor's avatar
csongor committed
338
339
            ax0.set_title(title)

Jait Dixit's avatar
Jait Dixit committed
340
341
342
343
344
        if (bool(kwargs.get("save", False))):
            fig.savefig(str(kwargs.get("save")), dpi=None, facecolor="none",
                        edgecolor="none", orientation="portrait",
                        papertype=None, format=None, transparent=False,
                        bbox_inches=None, pad_inches=0.1)
csongor's avatar
csongor committed
345
346
347
            pl.close(fig)
        else:
            fig.canvas.draw()