nifty_power.py 31.7 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
## 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
    ..                                  /______/

Marco Selig's avatar
Marco Selig committed
31
    NIFTY offers a number of automatized routines for handling
Marco Selig's avatar
Marco Selig committed
32
33
34
    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
Marco Selig's avatar
Marco Selig committed
35
    NIFTY, it is usually assumed that such a field follows statistical
Marco Selig's avatar
Marco Selig committed
36
37
38
    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
Marco Selig's avatar
Marco Selig committed
44
from scipy.interpolate import interp1d as ip ## FIXME: conflicts with sphinx's autodoc
45
#import numpy as np
46
from 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

Marco Selig's avatar
Marco Selig committed
294
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=True,var=10,force=False,bare=True,**kwargs):
295
    """
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
        smoothness : bool, *optional*
            Indicates whether the smoothness prior is used or not
Marco Selig's avatar
Marco Selig committed
338
            (default: True).
339
        var : {scalar, list, array}, *optional*
Marco Selig's avatar
Marco Selig committed
340
            Variance of the assumed spectral smoothness prior (default: 10).
Marco Selig's avatar
Marco Selig committed
341
        force : bool, *optional*, *experimental*
342
            Indicates whether smoothness is to be enforced or not
Marco Selig's avatar
Marco Selig committed
343
            (default: False).
344
345
346
347
        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
348

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

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

            - "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: "").
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
        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).
399
400
401
402

        Notes
        -----
        The general approach to inference of unknown power spectra is detailed
Marco Selig's avatar
Marco Selig committed
403
404
405
406
        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 [#]_,
407
        where the underlying formula(s), Eq.(26), of this implementation are
Marco Selig's avatar
Marco Selig committed
408
        derived and discussed in terms of their applicability.
409
410
411
412
413
414
415
416
417
418
419
420

        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
421

422
423
        Raises
        ------
424
        Exception, IndexError, TypeError, ValueError
425
            If some input is invalid.
Marco Selig's avatar
Marco Selig committed
426

427
428
429
430
431
432
433
    """
    ## 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):
434
            raise Exception(about._errors.cstring("ERROR: insufficient input."))
435
436
437
438
        else:
            domain = Sk.domain
    elif(not isinstance(domain,space)):
        raise TypeError(about._errors.cstring("ERROR: invalid input."))
439
440
    ## check implicit power indices
    if(pindex is None)or(kindex is None)or(rho is None):
441
        try:
442
            domain.set_power_indices(**kwargs)
443
        except:
444
            raise ValueError(about._errors.cstring("ERROR: invalid input."))
445
        else:
446
            pindex = domain.power_indices.get("pindex")
Marco Selig's avatar
Marco Selig committed
447
            pundex = domain.power_indices.get("pundex")
448
449
            kindex = domain.power_indices.get("kindex")
            rho = domain.power_indices.get("rho")
450
    ## check explicit power indices
451
452
    else:
        pindex = np.array(pindex,dtype=np.int)
453
454
        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))+" )."))
455
        kindex = np.sort(np.real(np.array(kindex,dtype=domain.vol.dtype).flatten(order='C')),axis=0,kind="quicksort",order=None)
456
457
458
        rho = np.array(rho,dtype=np.int)
        if(pundex is None):
            ## quick pundex
Marco Selig's avatar
Marco Selig committed
459
460
461
            pundex = np.unique(pindex,return_index=True,return_inverse=False)[1]
        else:
            pundex = np.array(pundex,dtype=np.int)
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
489
490
491
    ## 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
492

493
494
495
496
497
498
    ## 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
499
500
        if(np.any(trB2<0)):
            about.warnings.cprint("WARNING: nonpositive value(s) in tr[DSk] reset to 0.")
Marco Selig's avatar
Marco Selig committed
501
            trB2 = np.maximum(0,trB2)
502
503
504
505
506
507
    ## power spectrum
    numerator = 2*q+trB1+perception[0]*trB2 ## non-bare(!)
    denominator1 = rho+2*(alpha-1+perception[1])

    if(smoothness):
        if(not domain.discrete):
Marco Selig's avatar
Marco Selig committed
508
            numerator = weight_power(domain,numerator,power=-1,pindex=pindex,pundex=pundex) ## bare(!)
509
510

        ## smoothness prior
Marco Selig's avatar
Marco Selig committed
511
        permill = 0
Marco Selig's avatar
Marco Selig committed
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
        divergence = 1
        while(divergence):
            pk = numerator/denominator1 ## bare(!)
            tk = np.log(pk)
            Amemory = None
            var_ = var*1.1 # temporally increasing the variance
            var_OLD = -1
            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
Marco Selig's avatar
Marco Selig committed
531
                var_ /= 1.1+permill # lowering the variance when converged
Marco Selig's avatar
Marco Selig committed
532
533
534
535
536
537
538
539
540
541
542
543
544
                if(var_<var):
                    if(breakinfo): # making sure there's one iteration with the correct variance
                        break
                    var_ = var
                    breakinfo = True
                elif(var_==var_OLD):
                    if(divergence==3):
                        pot = int(np.log10(var_))
                        var = int(1+var_*10**-pot)*10**pot
                        about.warnings.cprint("WARNING: smoothness variance increased ( var = "+str(var)+" ).")
                        break
                    else:
                        divergence += 1
Marco Selig's avatar
Marco Selig committed
545
546
547
548
549
550
551
552
                        if(force):
                            permill = 0.001
                elif(force)and(var_/var_OLD>1.001):
                        permill = 0
                        pot = int(np.log10(var_))
                        var = int(1+var_*10**-pot)*10**pot
                        about.warnings.cprint("WARNING: smoothness variance increased ( var = "+str(var)+" ).")
                        break
Marco Selig's avatar
Marco Selig committed
553
554
555
556
                else:
                    var_OLD = var_
            if(breakinfo):
                break
557
558
559
560

        ## 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
561

562
563
564
    else:
        pk = numerator/denominator1 ## non-bare(!)
        ## weight if ...
Marco Selig's avatar
Marco Selig committed
565
566
        if(not domain.discrete)and(bare):
            pk = weight_power(domain,pk,power=-1,pindex=pindex,pundex=pundex) ## bare(!)
Marco Selig's avatar
Marco Selig committed
567

568
    return pk
Marco Selig's avatar
Marco Selig committed
569

570
##=============================================================================
Marco Selig's avatar
Marco Selig committed
571

572
573
574
575
576
577
578
579
##-----------------------------------------------------------------------------

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
580
        spec : {scalar, list, array, field, function}
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
639
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
            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:
666
                spec = np.array(spec(kindex),dtype=kindex.dtype)
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
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
            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

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