interpolate_power implemented; several minor adjustments.

nifty_core.py

@@ -2393,10 +2393,10 @@ class rg_space(space):
         if(not np.all(np.isfinite(spec))):
             about.warnings.cprint("WARNING: infinite value(s).")
         ## check positivity (excluding null)
-        if(np.any(spec==0)):
-            about.warnings.cprint("WARNING: nonpositive value(s).")
         if(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).")
 
         if(size is None):
             ## explicit kindex
@@ -3488,10 +3488,10 @@ class lm_space(space):
         if(not np.all(np.isfinite(spec))):
             about.warnings.cprint("WARNING: infinite value(s).")
         ## check positivity (excluding null)
-        if(np.any(spec==0)):
-            about.warnings.cprint("WARNING: nonpositive value(s).")
         if(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 = self.para[0]+1 ## lmax+1
         ## drop imaginary part
@@ -3689,7 +3689,7 @@ class lm_space(space):
             if(self.datatype==np.complex64):
                 x = gl.synalm_f(arg[1],lmax=lmax,mmax=lmax)
             else:
-                #x = gl.synalm(spec,lmax=lmax,mmax=lmax)
+                #x = gl.synalm(arg[1],lmax=lmax,mmax=lmax)
                 x = hp.synalm(arg[1],lmax=lmax,mmax=lmax)
             return x
 
@@ -4307,10 +4307,10 @@ class gl_space(space):
         if(not np.all(np.isfinite(spec))):
             about.warnings.cprint("WARNING: infinite value(s).")
         ## check positivity (excluding null)
-        if(np.any(spec==0)):
-            about.warnings.cprint("WARNING: nonpositive value(s).")
         if(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 = self.para[0] ## nlat
 
@@ -4932,10 +4932,10 @@ class hp_space(space):
         if(not np.all(np.isfinite(spec))):
             about.warnings.cprint("WARNING: infinite value(s).")
         ## check positivity (excluding null)
-        if(np.any(spec==0)):
-            about.warnings.cprint("WARNING: nonpositive value(s).")
         if(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 = 3*self.para[0] ## 3*nside
 

nifty_power.py

@@ -41,6 +41,7 @@
 
 """
 from __future__ import division
+from scipy.interpolate import interp1d as ip
 #import numpy as np
 from nifty_core import *
 import smoothing as gs
@@ -127,7 +128,7 @@ def weight_power(domain,spec,power=1,pindex=None,pundex=None,**kwargs):
         elif(len(pundex)!=np.size(domain.dim(split=True))):
             raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(len(pundex))+" <> "+str(np.size(domain.dim(split=True)))+" )."))
 
-    return np.real(domain.calc_weight(domain.enforce_power(spec,size=len(set(pindex.flatten(order='C'))))[pindex],power=power)[pundex])
+    return np.real(domain.calc_weight(domain.enforce_power(spec,size=np.max(pindex,axis=None,out=None)+1)[pindex],power=power)[pundex])
 
 ##-----------------------------------------------------------------------------
 
@@ -198,11 +199,47 @@ def smooth_power(spec,domain=None,kindex=None,mode="2s",exclude=1,sigma=-1,**kwa
                 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 explicit power indices
+        ## check power spectrum
+        spec = domain.enforce_power(spec,size=np.size(kindex))
+    ## check explicit kindex
     else:
-        kindex = np.array(kindex,dtype=domain.vol.dtype)
+        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:
+                spec = spec(kindex)
+            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]
+
     ## smoothing
     if(mode=="2s"):
         return gs.smooth_power_2s(spec,kindex,exclude=exclude,smooth_length=sigma)
@@ -381,7 +418,7 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None
 
         Raises
         ------
-        IndexError, TypeError, ValueError
+        Exception, IndexError, TypeError, ValueError
             If some input is invalid.
 
     """
@@ -391,7 +428,7 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None
     ## check domain
     if(domain is None):
         if(Sk is None):
-            raise TypeError(about._errors.cstring("ERROR: insufficient input."))
+            raise Exception(about._errors.cstring("ERROR: insufficient input."))
         else:
             domain = Sk.domain
     elif(not isinstance(domain,space)):
@@ -417,7 +454,7 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None
         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))+" )."))
-        kindex = np.array(kindex,dtype=domain.vol.dtype)
+        kindex = np.sort(np.real(np.array(kindex,dtype=domain.vol.dtype).flatten(order='C')),axis=0,kind="quicksort",order=None)
         rho = np.array(rho,dtype=np.int)
         if(pundex is None):
             ## quick pundex
@@ -509,3 +546,162 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None
 
 ##=============================================================================
 
+##-----------------------------------------------------------------------------
+
+def interpolate_power(spec,mode="linear",domain=None,kindex=None,newkindex=None,**kwargs):
+    """
+        Interpolates a given power spectrum at new k(-indices).
+
+        Parameters
+        ----------
+        spec : {scalar, array}
+            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:
+                spec = spec(kindex)
+            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
+
+##-----------------------------------------------------------------------------
+

smoothing.py

@@ -64,8 +64,8 @@ def smooth_power(power, k, exclude=1, smooth_length=None):
     pbinned = pbinned / counter
 
     nmirror = int(5 * smooth_length / dk) + 2
-    zpbinned = np.r_[(2 * pbinned[0] - pbinned[1:nmirror][::-1]), pbinned,
-                     (2 * pbinned[-1] - pbinned[-nmirror:-1][::-1])]
+    zpbinned = np.r_[np.exp(2 * np.log(pbinned[0]) - np.log(pbinned[1:nmirror][::-1])), pbinned,
+                     np.exp(2 * np.log(pbinned[-1])- np.log(pbinned[-nmirror:-1][::-1]))]
     zpbinned = np.maximum(0, zpbinned)
 
     tpbinned = np.fft.fftshift(np.fft.fft(zpbinned))
@@ -101,7 +101,7 @@ def smooth_power_bf(power, k, exclude=1, smooth_length=None):
         smooth_length = k[1]-k[0]
 
     nmirror = int(5*smooth_length/(k[1]-k[0]))+2
-    mpower = np.r_[(2*power[0]-power[1:nmirror][::-1]),power,(2*power[-1]-power[-nmirror:-1][::-1])]
+    mpower = np.r_[np.exp(2*np.log(power[0])-np.log(power[1:nmirror][::-1])),power,np.exp(2*np.log(power[-1])-np.log(power[-nmirror:-1][::-1]))]
     mk = np.r_[(2*k[0]-k[1:nmirror][::-1]),k,(2*k[-1]-k[-nmirror:-1][::-1])]
     mdk = np.r_[0.5*(mk[1]-mk[0]),0.5*(mk[2:]-mk[:-2]),0.5*(mk[-1]-mk[-2])]
 
@@ -143,7 +143,7 @@ def smooth_power_2s(power, k, exclude=1, smooth_length=None):
         smooth_length = k[1]-k[0]
 
     nmirror = int(5*smooth_length/(k[1]-k[0]))+2
-    mpower = np.r_[(2*power[0]-power[1:nmirror][::-1]),power,(2*power[-1]-power[-nmirror:-1][::-1])]
+    mpower = np.r_[np.exp(2*np.log(power[0])-np.log(power[1:nmirror][::-1])),power,np.exp(2*np.log(power[-1])-np.log(power[-nmirror:-1][::-1]))]
     mk = np.r_[(2*k[0]-k[1:nmirror][::-1]),k,(2*k[-1]-k[-nmirror:-1][::-1])]
     mdk = np.r_[0.5*(mk[1]-mk[0]),0.5*(mk[2:]-mk[:-2]),0.5*(mk[-1]-mk[-2])]