nifty_random.py 12.9 KB
Newer Older
1
2
# NIFTY (Numerical Information Field Theory) has been developed at the
# Max-Planck-Institute for Astrophysics.
Ultimanet's avatar
Ultimanet committed
3
##
4
# Copyright (C) 2013 Max-Planck-Society
Ultimanet's avatar
Ultimanet committed
5
##
6
7
# Author: Marco Selig
# Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
Ultimanet's avatar
Ultimanet committed
8
##
9
10
11
12
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
Ultimanet's avatar
Ultimanet committed
13
##
14
15
16
17
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# See the GNU General Public License for more details.
Ultimanet's avatar
Ultimanet committed
18
##
19
20
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
Ultimanet's avatar
Ultimanet committed
21
22

import numpy as np
23
from keepers import about
Ultimanet's avatar
Ultimanet committed
24
25
from nifty_mpi_data import distributed_data_object

26
27
# -----------------------------------------------------------------------------

Ultimanet's avatar
Ultimanet committed
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

class random(object):
    """
        ..                                          __
        ..                                        /  /
        ..       _____   ____ __   __ ___    ____/  /  ______    __ ____ ___
        ..     /   __/ /   _   / /   _   | /   _   / /   _   | /   _    _   |
        ..    /  /    /  /_/  / /  / /  / /  /_/  / /  /_/  / /  / /  / /  /
        ..   /__/     \______| /__/ /__/  \______|  \______/ /__/ /__/ /__/  class

        NIFTY (static) class for pseudo random number generators.

    """
    __init__ = None

43
    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Ultimanet's avatar
Ultimanet committed
44
45

    @staticmethod
46
    def parse_arguments(domain, **kwargs):
Ultimanet's avatar
Ultimanet committed
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
        """
            Analyses the keyword arguments for supported or necessary ones.

            Parameters
            ----------
            domain : space
                Space wherein the random field values live.
            random : string, *optional*
                Specifies a certain distribution to be drawn from using a
                pseudo random number generator. Supported distributions are:

                - "pm1" (uniform distribution over {+1,-1} or {+1,+i,-1,-i}
                - "gau" (normal distribution with zero-mean and a given
                    standard deviation or variance)
                - "syn" (synthesizes from a given power spectrum)
                - "uni" (uniform distribution over [vmin,vmax[)

            dev : {scalar, list, ndarray, field}, *optional*
                Standard deviation of the normal distribution if
                ``random == "gau"`` (default: None).
            var : {scalar, list, ndarray, field}, *optional*
                Variance of the normal distribution (outranks the standard
                deviation) if ``random == "gau"`` (default: None).
            spec : {scalar, list, array, field, function}, *optional*
                Power spectrum for ``random == "syn"`` (default: 1).
            size : integer, *optional*
                Number of irreducible bands for ``random == "syn"``
                (default: None).
            pindex : numpy.ndarray, *optional*
                Indexing array giving the power spectrum index of each band
                (default: None).
            kindex : numpy.ndarray, *optional*
                Scale of each irreducible band (default: None).
            vmax : {scalar, list, ndarray, field}, *optional*
                Upper limit of the uniform distribution if ``random == "uni"``
                (default: 1).

            Returns
            -------
            arg : list
                Ordered list of arguments (to be processed in
                ``get_random_values`` of the domain).

            Other Parameters
            ----------------
            codomain : nifty.space, *optional*
                A compatible codomain for power indexing (default: None).
            log : bool, *optional*
95
96
                Flag specifying if the spectral binning is performed on
                logarithmic
Ultimanet's avatar
Ultimanet committed
97
98
99
100
                scale or not; if set, the number of used bins is set
                automatically (if not given otherwise); by default no binning
                is done (default: None).
            nbin : integer, *optional*
101
102
                Number of used spectral bins; if given `log` is set to
                ``False``;
Ultimanet's avatar
Ultimanet committed
103
104
105
106
107
                integers below the minimum of 3 induce an automatic setting;
                by default no binning is done (default: None).
            binbounds : {list, array}, *optional*
                User specific inner boundaries of the bins, which are preferred
                over the above parameters; by default no binning is done
108
109
                (default: None).
            vmin : {scalar, list, ndarray, field}, *optional*
Ultimanet's avatar
Ultimanet committed
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
                Lower limit of the uniform distribution if ``random == "uni"``
                (default: 0).

            Raises
            ------
            KeyError
                If the `random` key is not supporrted.

        """
        if "random" in kwargs:
            key = kwargs.get("random")
        else:
            return None

        if key == "pm1":
Ultima's avatar
Ultima committed
125
            return {'random': key}
Ultimanet's avatar
Ultimanet committed
126
127
128

        elif key == "gau":
            mean = kwargs.get('mean', None)
Ultima's avatar
Ultima committed
129
130
131
132
            std = kwargs.get('std', None)
            return {'random': key,
                    'mean': mean,
                    'std': std}
Ultimanet's avatar
Ultimanet committed
133
134
135

        elif key == "syn":
            pindex = kwargs.get('pindex', None)
136
            kindex = kwargs.get('kindex', None)
Ultimanet's avatar
Ultimanet committed
137
138
139
140
141
142
            size = kwargs.get('size', None)
            log = kwargs.get('log', False)
            nbin = kwargs.get('nbin', None)
            binbounds = kwargs.get('binbounds', None)
            spec = kwargs.get('spec', 1)
            codomain = kwargs.get('codomain', None)
143

144
145
146
147
148
            # check which domain should be taken for powerindexing
            if domain.check_codomain(codomain) and codomain.harmonic:
                harmonic_domain = codomain
            elif domain.harmonic:
                harmonic_domain = domain
Ultimanet's avatar
Ultimanet committed
149
            else:
150
                harmonic_domain = domain.get_codomain()
151

152
153
            # building kpack
            if pindex is not None and kindex is not None:
154
                pindex = distributed_data_object(pindex,
Ultimanet's avatar
Ultimanet committed
155
156
157
158
                                                 distribution_strategy='fftw')
                kpack = [pindex, kindex]
            else:
                kpack = None
159

160
161
162
163
164
            # simply put size and kindex into enforce_power
            # if one or both are None, enforce power will fix that
            spec = harmonic_domain.enforce_power(spec,
                                                 size=size,
                                                 kindex=kindex)
Ultimanet's avatar
Ultimanet committed
165

Ultima's avatar
Ultima committed
166
167
168
169
170
171
172
            return {'random': key,
                    'spec': spec,
                    'kpack': kpack,
                    'harmonic_domain': harmonic_domain,
                    'log': log,
                    'nbin': nbin,
                    'binbounds': binbounds}
Ultimanet's avatar
Ultimanet committed
173
174

        elif key == "uni":
Ultima's avatar
Ultima committed
175
176
177
178
179
            vmin = domain.dtype.type(kwargs.get('vmin', 0))
            vmax = domain.dtype.type(kwargs.get('vmax', 1))
            return {'random': key,
                    'vmin': vmin,
                    'vmax': vmax}
Ultimanet's avatar
Ultimanet committed
180
181

        else:
182
183
            raise KeyError(about._errors.cstring(
                "ERROR: unsupported random key '" + str(key) + "'."))
Ultimanet's avatar
Ultimanet committed
184

185
    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Ultimanet's avatar
Ultimanet committed
186
187

    @staticmethod
188
    def pm1(dtype=np.dtype('int'), shape=1):
Ultimanet's avatar
Ultimanet committed
189
190
191
192
193
194
        """
            Generates random field values according to an uniform distribution
            over {+1,-1} or {+1,+i,-1,-i}, respectively.

            Parameters
            ----------
195
196
            dtype : type, *optional*
                Data type of the field values (default: int).
Ultimanet's avatar
Ultimanet committed
197
198
199
200
201
202
203
204
205
            shape : {integer, tuple, list, ndarray}, *optional*
                Split up dimension of the space (default: 1).

            Returns
            -------
            x : ndarray
                Random field values (with correct dtype and shape).

        """
206
        size = np.prod(shape, axis=0, dtype=np.dtype('int'), out=None)
Ultimanet's avatar
Ultimanet committed
207

Ultima's avatar
Ultima committed
208
        if issubclass(dtype.type, np.complexfloating):
209
210
211
212
            x = np.array([1 + 0j, 0 + 1j, -1 + 0j, 0 - 1j],
                         dtype=dtype)[np.random.randint(4,
                                                        high=None,
                                                        size=size)]
Ultimanet's avatar
Ultimanet committed
213
        else:
214
            x = 2 * np.random.randint(2, high=None, size=size) - 1
Ultimanet's avatar
Ultimanet committed
215

Ultima's avatar
Ultima committed
216
        return x.astype(dtype).reshape(shape)
Ultimanet's avatar
Ultimanet committed
217

218
    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Ultimanet's avatar
Ultimanet committed
219
220

    @staticmethod
Ultima's avatar
Ultima committed
221
    def gau(dtype=np.dtype('float64'), shape=(1,), mean=None, std=None):
Ultimanet's avatar
Ultimanet committed
222
223
224
225
226
        """
            Generates random field values according to a normal distribution.

            Parameters
            ----------
227
228
            dtype : type, *optional*
                Data type of the field values (default: float64).
Ultimanet's avatar
Ultimanet committed
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
            shape : {integer, tuple, list, ndarray}, *optional*
                Split up dimension of the space (default: 1).
            mean : {scalar, ndarray}, *optional*
                Mean of the normal distribution (default: 0).
            dev : {scalar, ndarray}, *optional*
                Standard deviation of the normal distribution (default: 1).
            var : {scalar, ndarray}, *optional*
                Variance of the normal distribution (outranks the standard
                deviation) (default: None).

            Returns
            -------
            x : ndarray
                Random field values (with correct dtype and shape).

            Raises
            ------
            ValueError
                If the array dimension of `mean`, `dev` or `var` mismatch with
                `shape`.

        """
Ultima's avatar
Ultima committed
251
        size = np.prod(shape)
Ultimanet's avatar
Ultimanet committed
252

Ultima's avatar
Ultima committed
253
254
        if issubclass(dtype.type, np.complexfloating):
            x = np.empty(size, dtype=dtype)
255
256
            x.real = np.random.normal(loc=0, scale=np.sqrt(0.5), size=size)
            x.imag = np.random.normal(loc=0, scale=np.sqrt(0.5), size=size)
Ultimanet's avatar
Ultimanet committed
257
        else:
258
            x = np.random.normal(loc=0, scale=1, size=size)
Ultimanet's avatar
Ultimanet committed
259

Ultima's avatar
Ultima committed
260
261
262
263
264
        if std is not None:
            if np.size(std) == 1:
                x *= np.abs(std)
            elif np.size(std) == size:
                x *= np.absolute(std).flatten()
Ultimanet's avatar
Ultimanet committed
265
            else:
266
                raise ValueError(about._errors.cstring(
Ultima's avatar
Ultima committed
267
                    "ERROR: dimension mismatch ( " + str(np.size(std)) +
268
                    " <> " + str(size) + " )."))
Ultima's avatar
Ultima committed
269
270
271

        if mean is not None:
            if np.size(mean) == 1:
Ultimanet's avatar
Ultimanet committed
272
                x += mean
Ultima's avatar
Ultima committed
273
            elif np.size(mean) == size:
Ultimanet's avatar
Ultimanet committed
274
275
                x += np.array(mean).flatten(order='C')
            else:
276
277
278
                raise ValueError(about._errors.cstring(
                    "ERROR: dimension mismatch ( " + str(np.size(mean)) +
                    " <> " + str(size) + " )."))
Ultimanet's avatar
Ultimanet committed
279

Ultima's avatar
Ultima committed
280
        return x.astype(dtype).reshape(shape)
Ultimanet's avatar
Ultimanet committed
281

282
    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Ultimanet's avatar
Ultimanet committed
283
284

    @staticmethod
285
    def uni(dtype=np.dtype('float64'), shape=1, vmin=0, vmax=1):
Ultimanet's avatar
Ultimanet committed
286
287
288
289
290
291
        """
            Generates random field values according to an uniform distribution
            over [vmin,vmax[.

            Parameters
            ----------
292
293
            dtype : type, *optional*
                Data type of the field values (default: float64).
Ultimanet's avatar
Ultimanet committed
294
295
296
297
298
299
300
301
302
303
304
305
306
307
            shape : {integer, tuple, list, ndarray}, *optional*
                Split up dimension of the space (default: 1).

            vmin : {scalar, list, ndarray, field}, *optional*
                Lower limit of the uniform distribution (default: 0).
            vmax : {scalar, list, ndarray, field}, *optional*
                Upper limit of the uniform distribution (default: 1).

            Returns
            -------
            x : ndarray
                Random field values (with correct dtype and shape).

        """
308
309
        size = np.prod(shape, axis=0, dtype=np.dtype('int'), out=None)
        if(np.size(vmin) > 1):
Ultimanet's avatar
Ultimanet committed
310
            vmin = np.array(vmin).flatten(order='C')
311
        if(np.size(vmax) > 1):
Ultimanet's avatar
Ultimanet committed
312
313
            vmax = np.array(vmax).flatten(order='C')

314
315
316
317
318
319
320
321
        if(dtype in [np.dtype('complex64'), np.dtype('complex128')]):
            x = np.empty(size, dtype=dtype, order='C')
            x.real = (vmax - vmin) * np.random.random(size=size) + vmin
            x.imag = (vmax - vmin) * np.random.random(size=size) + vmin
        elif(dtype in [np.dtype('int8'), np.dtype('int16'), np.dtype('int32'),
                       np.dtype('int64')]):
            x = np.random.randint(
                min(vmin, vmax), high=max(vmin, vmax), size=size)
Ultimanet's avatar
Ultimanet committed
322
        else:
323
            x = (vmax - vmin) * np.random.random(size=size) + vmin
Ultimanet's avatar
Ultimanet committed
324

325
        return x.astype(dtype).reshape(shape, order='C')
Ultimanet's avatar
Ultimanet committed
326

327
    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Ultimanet's avatar
Ultimanet committed
328
329
330
331

    def __repr__(self):
        return "<nifty_core.random>"

332
# -----------------------------------------------------------------------------