nifty_power.py 30.4 KB
Newer Older
Marco Selig's avatar
Marco Selig committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
## NIFTY (Numerical Information Field Theory) has been developed at the
## Max-Planck-Institute for Astrophysics.
##
## Copyright (C) 2013 Max-Planck-Society
##
## Author: Marco Selig
## Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
##
## 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.
##
## 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.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.

"""
    ..                     __   ____   __
    ..                   /__/ /   _/ /  /_
    ..         __ ___    __  /  /_  /   _/  __   __
    ..       /   _   | /  / /   _/ /  /   /  / /  /
    ..      /  / /  / /  / /  /   /  /_  /  /_/  /
    ..     /__/ /__/ /__/ /__/    \___/  \___   /  power
    ..                                  /______/

    NIFTy offers a number of automatized routines for handling
    power spectra. It is possible to draw a field from a random distribution
    with a certain autocorrelation or, equivalently, with a certain
    power spectrum in its conjugate space (see :py:func:`field.random`). In
    NIFTy, it is usually assumed that such a field follows statistical
    homogeneity and isotropy. Fields which are only statistically homogeneous
    can also be created using the diagonal operator routine.

39
40
    At the moment, NIFTY offers several additional routines for power spectrum
    manipulation.
Marco Selig's avatar
Marco Selig committed
41
42
43

"""
from __future__ import division
44
from scipy.interpolate import interp1d as ip ## conflicts with sphinx's autodoc
45
#import numpy as np
46
from nifty.nifty_core import *
Marco Selig's avatar
Marco Selig committed
47
48
49
import smoothing as gs


50
51
##-----------------------------------------------------------------------------

52
def weight_power(domain,spec,power=1,pindex=None,pundex=None,**kwargs):
53
54
55
56
57
58
59
60
    """
        Weights a given power spectrum with the corresponding pixel volumes
        to a given power.

        Parameters
        ----------
        domain : space
            The space wherein valid arguments live.
Marco Selig's avatar
Marco Selig committed
61
        spec : {scalar, list, array, field, function}
62
63
            The power spectrum. A scalars is interpreted as a constant
            spectrum.
64
        pindex : ndarray, *optional*
65
66
            Indexing array giving the power spectrum index for each
            represented mode.
Marco Selig's avatar
Marco Selig committed
67
68
        pundex : ndarray, *optional*
            Unindexing array undoing power indexing.
69
70
71
72
73
74

        Returns
        -------
        spev : ndarray
            Weighted power spectrum.

75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
        Other Parameters
        ----------------
        log : bool, *optional*
            Flag specifying if the spectral binning is performed on logarithmic
            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*
            Number of used spectral bins; if given `log` is set to ``False``;
            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
            (default: None).            vmin : {scalar, list, ndarray, field}, *optional*
            Lower limit of the uniform distribution if ``random == "uni"``
            (default: 0).

93
94
95
96
97
98
99
100
101
102
103
        Raises
        ------
        TypeError
            If `domain` is no space.
        ValueError
            If `domain` is no harmonic space.

    """
    ## check domain
    if(not isinstance(domain,space)):
        raise TypeError(about._errors.cstring("ERROR: invalid input."))
104
    ## check implicit power indices
105
106
    if(pindex is None):
        try:
107
108
            domain.set_power_indices(**kwargs)
        except:
109
            raise ValueError(about._errors.cstring("ERROR: invalid input."))
110
111
112
113
        else:
            pindex = domain.power_indices.get("pindex")
            if(pundex is None):
                pundex = domain.power_indices.get("pundex")
Marco Selig's avatar
Marco Selig committed
114
115
116
117
            else:
                pundex = np.array(pundex,dtype=np.int)
                if(np.size(pundex)!=np.max(pindex,axis=None,out=None)+1):
                    raise ValueError(about._errors.cstring("ERROR: size mismatch ( "+str(np.size(pundex))+" <> "+str(np.max(pindex,axis=None,out=None)+1)+" )."))
118
119
120
121
122
123
124
    ## check explicit power indices
    else:
        pindex = np.array(pindex,dtype=np.int)
        if(not np.all(np.array(np.shape(pindex))==domain.dim(split=True))):
            raise ValueError(about._errors.cstring("ERROR: shape mismatch ( "+str(np.array(np.shape(pindex)))+" <> "+str(domain.dim(split=True))+" )."))
        if(pundex is None):
            ## quick pundex
Marco Selig's avatar
Marco Selig committed
125
126
127
128
129
            pundex = np.unique(pindex,return_index=True,return_inverse=False)[1]
        else:
            pundex = np.array(pundex,dtype=np.int)
            if(np.size(pundex)!=np.max(pindex,axis=None,out=None)+1):
                raise ValueError(about._errors.cstring("ERROR: size mismatch ( "+str(np.size(pundex))+" <> "+str(np.max(pindex,axis=None,out=None)+1)+" )."))
130

Marco Selig's avatar
Marco Selig committed
131
    return np.real(domain.calc_weight(domain.enforce_power(spec,size=np.max(pindex,axis=None,out=None)+1)[pindex],power=power).flatten(order='C')[pundex])
132

Marco Selig's avatar
Marco Selig committed
133
134
##-----------------------------------------------------------------------------

135
136
##-----------------------------------------------------------------------------

137
def smooth_power(spec,domain=None,kindex=None,mode="2s",exclude=1,sigma=-1,**kwargs):
Marco Selig's avatar
Marco Selig committed
138
    """
139
        Smoothes a power spectrum via convolution with a Gaussian kernel.
Marco Selig's avatar
Marco Selig committed
140

141
142
        Parameters
        ----------
Marco Selig's avatar
Marco Selig committed
143
        spec : {scalar, list, array, field, function}
144
145
146
147
148
149
150
151
            The power spectrum to be smoothed.
        domain : space, *optional*
            The space wherein the power spectrum is defined (default: None).
        kindex : ndarray, *optional*
            The array specifying the coordinate indices in conjugate space
            (default: None).
        mode : string, *optional*
            Specifies the smoothing mode (default: "2s") :
Marco Selig's avatar
Marco Selig committed
152

153
154
155
            - "ff" (smoothing in the harmonic basis using fast Fourier transformations)
            - "bf" (smoothing in the position basis by brute force)
            - "2s" (smoothing in the position basis restricted to a 2-`sigma` interval)
Marco Selig's avatar
Marco Selig committed
156

157
158
159
160
161
162
        exclude : scalar, *optional*
            Excludes the first power spectrum entries from smoothing, indicated by
            the given integer scalar (default = 1, the monopol is not smoothed).
        sigma : scalar, *optional*
            FWHM of Gaussian convolution kernel (default = -1, `sigma` is set
            automatically).
Marco Selig's avatar
Marco Selig committed
163

164
165
166
167
        Returns
        -------
        smoothspec : ndarray
            The smoothed power spectrum.
Marco Selig's avatar
Marco Selig committed
168

169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
        Other Parameters
        ----------------
        log : bool, *optional*
            Flag specifying if the spectral binning is performed on logarithmic
            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*
            Number of used spectral bins; if given `log` is set to ``False``;
            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
            (default: None).            vmin : {scalar, list, ndarray, field}, *optional*
            Lower limit of the uniform distribution if ``random == "uni"``
            (default: 0).
Marco Selig's avatar
Marco Selig committed
186

187
188
189
190
        Raises
        ------
        KeyError
            If `mode` is unsupported.
Marco Selig's avatar
Marco Selig committed
191

Marco Selig's avatar
Marco Selig committed
192
    """
193
194
195
196
197
198
199
200
201
    ## check implicit kindex
    if(kindex is None):
        if(isinstance(domain,space)):
            try:
                domain.set_power_indices(**kwargs)
            except:
                raise ValueError(about._errors.cstring("ERROR: invalid input."))
            else:
                kindex = domain.power_indices.get("kindex")
202

203
204
        else:
            raise TypeError(about._errors.cstring("ERROR: insufficient input."))
205
206
207
        ## check power spectrum
        spec = domain.enforce_power(spec,size=np.size(kindex))
    ## check explicit kindex
208
    else:
209
210
211
212
213
214
215
        kindex = np.sort(np.real(np.array(kindex,dtype=None).flatten(order='C')),axis=0,kind="quicksort",order=None)

        ## check power spectrum
        if(isinstance(spec,field)):
            spec = spec.val.astype(kindex.dtype)
        elif(callable(spec)):
            try:
216
                spec = np.array(spec(kindex),dtype=kindex.dtype)
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
            except:
                TypeError(about._errors.cstring("ERROR: invalid power spectra function.")) ## exception in ``spec(kindex)``
        elif(np.isscalar(spec)):
            spec = np.array([spec],dtype=kindex.dtype)
        else:
            spec = np.array(spec,dtype=kindex.dtype)
        ## drop imaginary part
        spec = np.real(spec)
        ## check finiteness and positivity (excluding null)
        if(not np.all(np.isfinite(spec))):
            raise ValueError(about._errors.cstring("ERROR: infinite value(s)."))
        elif(np.any(spec<0)):
            raise ValueError(about._errors.cstring("ERROR: nonpositive value(s)."))
        elif(np.any(spec==0)):
            about.warnings.cprint("WARNING: nonpositive value(s).")
        size = np.size(kindex)
        ## extend
        if(np.size(spec)==1):
            spec = spec*np.ones(size,dtype=spec.dtype,order='C')
        ## size check
        elif(np.size(spec)<size):
            raise ValueError(about._errors.cstring("ERROR: size mismatch ( "+str(np.size(spec))+" < "+str(size)+" )."))
        elif(np.size(spec)>size):
            about.warnings.cprint("WARNING: power spectrum cut to size ( == "+str(size)+" ).")
            spec = spec[:size]

243
    ## smoothing
Marco Selig's avatar
Marco Selig committed
244
    if(mode=="2s"):
245
        return gs.smooth_power_2s(spec,kindex,exclude=exclude,smooth_length=sigma)
Marco Selig's avatar
Marco Selig committed
246
    elif(mode=="ff"):
247
        return gs.smooth_power(spec,kindex,exclude=exclude,smooth_length=sigma)
Marco Selig's avatar
Marco Selig committed
248
    elif(mode=="bf"):
249
        return gs.smooth_power_bf(spec,kindex,exclude=exclude,smooth_length=sigma)
Marco Selig's avatar
Marco Selig committed
250
251
    else:
        raise KeyError(about._errors.cstring("ERROR: unsupported mode '"+str(mode)+"'."))
Marco Selig's avatar
Marco Selig committed
252
253
254

##-----------------------------------------------------------------------------

255
256
257
258
259
260
261
262
263
264
265
266
267
##=============================================================================

def _calc_laplace(kindex): ## > computes Laplace operator and integrand
    ## finite differences
    l = np.r_[0,0,np.log(kindex[2:]/kindex[1])]
    dl1 = l[1:]-l[:-1]
    dl2 = l[2:]-l[:-2]
    if(np.any(dl1[1:]==0))or(np.any(dl2==0)):
        raise ValueError(about._errors.cstring("ERROR: too finely divided harmonic grid."))
    ## operator(s)
    klim = len(kindex)
    L = np.zeros((klim,klim))
    I = np.zeros(klim)
Marco Selig's avatar
Marco Selig committed
268
    for jj in xrange(2,klim-1): ## leave out {0,1,kmax}
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
        L[jj,jj-1] = 2/(dl2[jj-1]*dl1[jj-1])
        L[jj,jj] = -2/dl2[jj-1]*(1/dl1[jj]+1/dl1[jj-1])
        L[jj,jj+1] = 2/(dl2[jj-1]*dl1[jj])
        I[jj] = dl2[jj-1]/2
    return L,I

def _calc_inverse(tk,var,kindex,rho,b1,Amem): ## > computes the inverse Hessian `A` and `b2`
    ## operator `T` from Eq.(B8) times 2
    if(Amem is None):
        L,I = _calc_laplace(kindex)
        #T2 = 2*np.dot(L.T,np.dot(np.diag(I/var,k=0),L,out=None),out=None) # Eq.(B8) * 2
        if(np.isscalar(var)):
            Amem = np.dot(L.T,np.dot(np.diag(I,k=0),L,out=None),out=None)
            T2 = 2/var*Amem
        else:
            Amem = np.dot(np.diag(np.sqrt(I),k=0),L,out=None)
            T2 = 2*np.dot(Amem.T,np.dot(np.diag(1/var,k=0),Amem,out=None),out=None)
    elif(np.isscalar(var)):
        T2 = 2/var*Amem
    else:
        T2 = 2*np.dot(Amem.T,np.dot(np.diag(1/var,k=0),Amem,out=None),out=None)
    b2 = b1+np.dot(T2,tk,out=None)
    ## inversion
    return np.linalg.inv(T2+np.diag(b2,k=0)),b2,Amem
Marco Selig's avatar
Marco Selig committed
293

294
295
def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None,rho=None,q=1E-42,alpha=1,perception=(1,0),smoothness=False,var=100,bare=True,**kwargs):
    """
Marco Selig's avatar
Marco Selig committed
296
        Infers the power spectrum.
Marco Selig's avatar
Marco Selig committed
297

Marco Selig's avatar
Marco Selig committed
298
        Given a map the inferred power spectrum is equal to ``m.power()``; given
299
300
301
        an uncertainty a power spectrum is inferred according to the "critical"
        filter formula, which can be extended by a smoothness prior. For
        details, see references below.
Marco Selig's avatar
Marco Selig committed
302

303
304
305
        Parameters
        ----------
        m : field
Marco Selig's avatar
Marco Selig committed
306
            Map for which the power spectrum is inferred.
307
308
309
310
311
312
313
314
315
316
317
318
319
        domain : space
            The space wherein the power spectrum is defined, can be retrieved
            from `Sk.domain` (default: None).
        Sk : projection_operator
            Projection operator specifying the pseudo trace for all projection
            bands, can be initialized from `domain` and `pindex`
            (default: None).
        D : operator, *optional*
            Operator expressing the uncertainty of the map `m`, its diagonal
            `D.hathat` in the `domain` suffices (default: 0).
        pindex : numpy.ndarray, *optional*
            Indexing array giving the power spectrum index for each
            represented mode (default: None).
Marco Selig's avatar
Marco Selig committed
320
321
        pundex : ndarray, *optional*
            Unindexing array undoing power indexing.
322
323
324
325
326
327
328
329
330
331
332
333
        kindex : numpy.ndarray, *optional*
            Scale corresponding to each band in the power spectrum
            (default: None).
        rho : numpy.ndarray, *optional*
            Number of modes per scale (default: None).
        q : {scalar, list, array}, *optional*
            Spectral scale parameter of the assumed inverse-Gamme prior
            (default: 1E-42).
        alpha : {scalar, list, array}, *optional*
            Spectral shape parameter of the assumed inverse-Gamme prior
            (default: 1).
        perception : {tuple, list, array}, *optional*
Marco Selig's avatar
Marco Selig committed
334
335
            Tuple specifying the filter perception (delta,epsilon)
            (default: (1,0)).
336
337
338
339
340
341
342
343
344
        smoothness : bool, *optional*
            Indicates whether the smoothness prior is used or not
            (default: False).
        var : {scalar, list, array}, *optional*
            Variance of the assumed spectral smoothness prior (default: 100).
        bare : bool, *optional*
            Indicates whether the power spectrum entries returned are "bare"
            or not (mandatory for the correct incorporation of volume weights)
            (default: True).
Marco Selig's avatar
Marco Selig committed
345

346
347
348
        Returns
        -------
        pk : numpy.ndarray
Marco Selig's avatar
Marco Selig committed
349
            The inferred power spectrum, weighted according to the `bare` flag.
350
351
352
353

        Other Parameters
        ----------------
        random : string, *optional*
Marco Selig's avatar
Marco Selig committed
354
355
            The distribution from which the probes for the diagonal probing are
            drawn, supported distributions are (default: "pm1"):
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380

            - "pm1" (uniform distribution over {+1,-1} or {+1,+i,-1,-i})
            - "gau" (normal distribution with zero-mean and unit-variance)

        ncpu : int, *optional*
            The number of CPUs to be used for parallel probing (default: 2).
        nrun : int, *optional*
            The number of probes to be evaluated; if ``nrun < ncpu ** 2``, it
            will be set to ``ncpu ** 2`` (default: 8).
        nper : int, *optional*
            This number specifies how many probes will be evaluated by one
            worker. Afterwards a new worker will be created to evaluate a chunk
            of `nper` probes; it is recommended to stay with the default value
            (default: None).
        save : bool, *optional*
            If `save` is True, then the probing results will be written to the
            hard disk instead of being saved in the RAM; this is recommended
            for high dimensional fields whose probes would otherwise fill up
            the memory (default: False).
        path : string, *optional*
            The path, where the probing results are saved, if `save` is True
            (default: "tmp").
        prefix : string, *optional*
            A prefix for the saved probing results; the saved results will be
            named using that prefix and an 8-digit number (default: "").
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
        log : bool, *optional*
            Flag specifying if the spectral binning is performed on logarithmic
            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*
            Number of used spectral bins; if given `log` is set to ``False``;
            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
            (default: None).            vmin : {scalar, list, ndarray, field}, *optional*
            Lower limit of the uniform distribution if ``random == "uni"``
            (default: 0).
396
397
398
399

        Notes
        -----
        The general approach to inference of unknown power spectra is detailed
Marco Selig's avatar
Marco Selig committed
400
401
402
403
404
405
        in [#]_, where the "critical" filter formula, Eq.(37b), used here is
        derived, and the implications of a certain choise of the perception
        tuple (delta,epsilon) are discussed.
        The further incorporation of a smoothness prior as detailed in [#]_,
        where the underlying formula(s), Eq.(27), of this implementation are
        derived and discussed in terms of their applicability.
406
407
408
409
410
411
412
413
414
415
416
417

        References
        ----------
        .. [#] T.A. Ensslin and M. Frommert, "Reconstruction of signals with
            unknown spectra in information field theory with parameter
            uncertainty", Physical Review E, 2011,
            10.1103/PhysRevD.83.105014;
            `arXiv:1002.2928 <http://www.arxiv.org/abs/1002.2928>`_
        .. [#] N. Opermann et. al., "Reconstruction of Gaussian and log-normal
            fields with spectral smoothness", Physical Review E, 2013,
            10.1103/PhysRevE.87.032136;
            `arXiv:1210.6866 <http://www.arxiv.org/abs/1210.6866>`_
Marco Selig's avatar
Marco Selig committed
418

419
420
        Raises
        ------
421
        Exception, IndexError, TypeError, ValueError
422
            If some input is invalid.
Marco Selig's avatar
Marco Selig committed
423

424
425
426
427
428
429
430
    """
    ## check map
    if(not isinstance(m,field)):
        raise TypeError(about._errors.cstring("ERROR: invalid input."))
    ## check domain
    if(domain is None):
        if(Sk is None):
431
            raise Exception(about._errors.cstring("ERROR: insufficient input."))
432
433
434
435
        else:
            domain = Sk.domain
    elif(not isinstance(domain,space)):
        raise TypeError(about._errors.cstring("ERROR: invalid input."))
436
437
    ## check implicit power indices
    if(pindex is None)or(kindex is None)or(rho is None):
438
        try:
439
            domain.set_power_indices(**kwargs)
440
        except:
441
            raise ValueError(about._errors.cstring("ERROR: invalid input."))
442
        else:
443
            pindex = domain.power_indices.get("pindex")
Marco Selig's avatar
Marco Selig committed
444
            pundex = domain.power_indices.get("pundex")
445
446
            kindex = domain.power_indices.get("kindex")
            rho = domain.power_indices.get("rho")
447
    ## check explicit power indices
448
449
    else:
        pindex = np.array(pindex,dtype=np.int)
450
451
        if(not np.all(np.array(np.shape(pindex))==domain.dim(split=True))):
            raise ValueError(about._errors.cstring("ERROR: shape mismatch ( "+str(np.array(np.shape(pindex)))+" <> "+str(domain.dim(split=True))+" )."))
452
        kindex = np.sort(np.real(np.array(kindex,dtype=domain.vol.dtype).flatten(order='C')),axis=0,kind="quicksort",order=None)
453
454
455
        rho = np.array(rho,dtype=np.int)
        if(pundex is None):
            ## quick pundex
Marco Selig's avatar
Marco Selig committed
456
457
458
            pundex = np.unique(pindex,return_index=True,return_inverse=False)[1]
        else:
            pundex = np.array(pundex,dtype=np.int)
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
    ## check projection operator
    if(Sk is None):
        Sk = projection_operator(domain,assign=pindex)
    elif(not isinstance(Sk,projection_operator))or(not hasattr(Sk,"pseudo_tr")):
        raise TypeError(about._errors.cstring("ERROR: invalid input."))
    elif(Sk.domain<>domain):
        raise ValueError(about._errors.cstring("ERROR: invalid input."))
    ## check critical parameters
    if(not np.isscalar(q)):
        q = np.array(q,dtype=domain.vol.dtype).flatten()
        if(np.size(q)<>np.size(kindex)):
            raise ValueError(about._errors.cstring("ERROR: invalid input."))
    if(not np.isscalar(alpha)):
        alpha = np.array(alpha,dtype=domain.vol.dtype).flatten()
        if(np.size(alpha)<>np.size(kindex)):
            raise ValueError(about._errors.cstring("ERROR: invalid input."))
    ## check perception (delta,epsilon)
    if(perception is None):
        perception = (1,0) ## critical perception
    elif(not isinstance(perception,(tuple,list,np.ndarray))):
        raise TypeError(about._errors.cstring("ERROR: invalid input."))
    elif(len(perception)<2):
        raise IndexError(about._errors.cstring("ERROR: invalid input."))
    if(perception[1] is None):
        perception[1] = rho/2*(perception[0]-1) ## critical epsilon
    ## check smothness variance
    if(not np.isscalar(var)):
        var = np.array(var,dtype=domain.vol.dtype).flatten()
        if(np.size(var)<>np.size(kindex)):
            raise ValueError(about._errors.cstring("ERROR: invalid input."))
Marco Selig's avatar
Marco Selig committed
489

490
491
492
493
494
495
    ## trace(s) of B
    trB1 = Sk.pseudo_tr(m) ## == Sk(m).pseudo_dot(m), but faster
    if(perception[0]==0)or(D is None)or(D==0):
        trB2 = 0
    else:
        trB2 = Sk.pseudo_tr(D,**kwargs) ## probing of the partial traces of D
496
497
498
        if(np.any(trB2<0)):
            about.warnings.cprint("WARNING: nonpositive value(s) in tr[DSk] reset to 0.")
            trB2 = np.minimum(0,trB2)
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
    ## power spectrum
    numerator = 2*q+trB1+perception[0]*trB2 ## non-bare(!)
    denominator1 = rho+2*(alpha-1+perception[1])

    if(smoothness):
        if(not domain.discrete):
            numerator = weight_power(domain,numerator,power=-1,pindex=pindex,pundex=pundex)
        pk = numerator/denominator1 ## bare(!)

        ## smoothness prior
        tk = np.log(pk)
        Amemory = None
        var_ = var*1.1 # temporally increasing the variance
        breakinfo = False
        while(var_>=var): # slowly lowering the variance
            absdelta = 1
            while(absdelta>1E-3): # solving with fixed variance
                ## solution of A delta = b1 - b2
                Ainverse,denominator2,Amemory = _calc_inverse(tk,var_,kindex,rho,denominator1,Amemory)
                delta = np.dot(Ainverse,numerator/pk-denominator2,out=None)
                if(np.abs(delta).max()>absdelta): # increasing variance when speeding up
                    var_ *= 1.1
                absdelta = np.abs(delta).max()
                tk += min(1,0.1/absdelta)*delta # adaptive step width
                pk *= np.exp(min(1,0.1/absdelta)*delta) # adaptive step width
            var_ /= 1.1 # lowering the variance when converged
            if(var_<var):
                if(breakinfo): # making sure there's one iteration with the correct variance
                    break
                var_ = var
                breakinfo = True

        ## weight if ...
        if(not domain.discrete)and(not bare):
            pk = weight_power(domain,pk,power=1,pindex=pindex,pundex=pundex) ## non-bare(!)
Marco Selig's avatar
Marco Selig committed
534

535
536
537
538
539
    else:
        pk = numerator/denominator1 ## non-bare(!)
        ## weight if ...
        if(not domain.discrete)and(not bare):
            pk = weight_power(domain,pk,power=1,pindex=pindex,pundex=pundex) ## non-bare(!)
Marco Selig's avatar
Marco Selig committed
540

541
    return pk
Marco Selig's avatar
Marco Selig committed
542

543
##=============================================================================
Marco Selig's avatar
Marco Selig committed
544

545
546
547
548
549
550
551
552
##-----------------------------------------------------------------------------

def interpolate_power(spec,mode="linear",domain=None,kindex=None,newkindex=None,**kwargs):
    """
        Interpolates a given power spectrum at new k(-indices).

        Parameters
        ----------
Marco Selig's avatar
Marco Selig committed
553
        spec : {scalar, list, array, field, function}
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
            The power spectrum. A scalars is interpreted as a constant
            spectrum.
        mode : string
            String specifying the interpolation scheme, supported
            schemes are (default: "linear"):

            - "linear"
            - "nearest"
            - "zero"
            - "slinear"
            - "quadratic"
            - "cubic"

        domain : space, *optional*
            The space wherein the power spectrum is defined (default: None).
        kindex : numpy.ndarray, *optional*
            Scales corresponding to each band in the old power spectrum;
            can be retrieved from `domain` (default: None).
        newkindex : numpy.ndarray, *optional*
            Scales corresponding to each band in the new power spectrum;
            can be retrieved from `domain` if `kindex` is given
            (default: None).

        Returns
        -------
        newspec : numpy.ndarray
            The interpolated power spectrum.

        Other Parameters
        ----------------
        log : bool, *optional*
            Flag specifying if the spectral binning is performed on logarithmic
            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*
            Number of used spectral bins; if given `log` is set to ``False``;
            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
            (default: None).            vmin : {scalar, list, ndarray, field}, *optional*
            Lower limit of the uniform distribution if ``random == "uni"``
            (default: 0).

        See Also
        --------
        scipy.interpolate.interp1d

        Raises
        ------
        Exception, IndexError, TypeError, ValueError
            If some input is invalid.
        ValueError
            If an interpolation is flawed.

    """
    ## check implicit kindex
    if(kindex is None):
        if(isinstance(domain,space)):
            try:
                domain.set_power_indices(**kwargs)
            except:
                raise ValueError(about._errors.cstring("ERROR: invalid input."))
            else:
                kindex = domain.power_indices.get("kindex")
        else:
            raise TypeError(about._errors.cstring("ERROR: insufficient input."))
        ## check power spectrum
        spec = domain.enforce_power(spec,size=np.size(kindex))
        ## check explicit newkindex
        if(newkindex is None):
            raise Exception(about._errors.cstring("ERROR: insufficient input."))
        else:
            newkindex = np.sort(np.real(np.array(newkindex,dtype=domain.vol.dtype).flatten(order='C')),axis=0,kind="quicksort",order=None)
    ## check explicit kindex
    else:
        kindex = np.sort(np.real(np.array(kindex,dtype=None).flatten(order='C')),axis=0,kind="quicksort",order=None)

        ## check power spectrum
        if(isinstance(spec,field)):
            spec = spec.val.astype(kindex.dtype)
        elif(callable(spec)):
            try:
639
                spec = np.array(spec(kindex),dtype=kindex.dtype)
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
            except:
                TypeError(about._errors.cstring("ERROR: invalid power spectra function.")) ## exception in ``spec(kindex)``
        elif(np.isscalar(spec)):
            spec = np.array([spec],dtype=kindex.dtype)
        else:
            spec = np.array(spec,dtype=kindex.dtype)
        ## drop imaginary part
        spec = np.real(spec)
        ## check finiteness and positivity (excluding null)
        if(not np.all(np.isfinite(spec))):
            raise ValueError(about._errors.cstring("ERROR: infinite value(s)."))
        elif(np.any(spec<0)):
            raise ValueError(about._errors.cstring("ERROR: nonpositive value(s)."))
        elif(np.any(spec==0)):
            about.warnings.cprint("WARNING: nonpositive value(s).")
        size = np.size(kindex)
        ## extend
        if(np.size(spec)==1):
            spec = spec*np.ones(size,dtype=spec.dtype,order='C')
        ## size check
        elif(np.size(spec)<size):
            raise ValueError(about._errors.cstring("ERROR: size mismatch ( "+str(np.size(spec))+" < "+str(size)+" )."))
        elif(np.size(spec)>size):
            about.warnings.cprint("WARNING: power spectrum cut to size ( == "+str(size)+" ).")
            spec = spec[:size]

        ## check implicit newkindex
        if(newkindex is None):
            if(isinstance(domain,space)):
                try:
                    domain.set_power_indices(**kwargs)
                except:
                    raise ValueError(about._errors.cstring("ERROR: invalid input."))
                else:
                    newkindex = domain.power_indices.get("kindex")
            else:
                raise TypeError(about._errors.cstring("ERROR: insufficient input."))
        ## check explicit newkindex
        else:
            newkindex = np.sort(np.real(np.array(newkindex,dtype=None).flatten(order='C')),axis=0,kind="quicksort",order=None)

    ## check bounds
    if(kindex[0]<0)or(newkindex[0]<0):
        raise ValueError(about._errors.cstring("ERROR: invalid input."))
    if(np.any(newkindex>kindex[-1])):
        about.warnings.cprint("WARNING: interpolation beyond upper bound.")
        ## continuation extension by point mirror
        nmirror = np.size(kindex)-np.searchsorted(kindex,2*kindex[-1]-newkindex[-1],side='left')+1
        spec = np.r_[spec,np.exp(2*np.log(spec[-1])-np.log(spec[-nmirror:-1][::-1]))]
        kindex = np.r_[kindex,(2*kindex[-1]-kindex[-nmirror:-1][::-1])]
    ## interpolation
    newspec = ip(kindex,spec,kind=mode,axis=0,copy=True,bounds_error=True,fill_value=np.NAN)(newkindex)
    ## check new power spectrum
    if(not np.all(np.isfinite(newspec))):
        raise ValueError(about._errors.cstring("ERROR: infinite value(s)."))
    elif(np.any(newspec<0)):
        raise ValueError(about._errors.cstring("ERROR: nonpositive value(s)."))
    elif(np.any(newspec==0)):
        about.warnings.cprint("WARNING: nonpositive value(s).")

    return newspec

##-----------------------------------------------------------------------------