From 7d9aaa5e09c979d55d26e7d4e9e2d1be5f7b9e61 Mon Sep 17 00:00:00 2001
From: Marco Selig <mselig@ncg-02.MPA-Garching.MPG.DE>
Date: Fri, 25 Apr 2014 16:52:27 +0200
Subject: [PATCH] demos updated.
---
demos/demo_excaliwir.py | 416 +++++++++++++++++-----------------------
demos/demo_faraday.py | 99 +++++-----
demos/demo_wf1.py | 37 ++--
demos/demo_wf2.py | 40 ++--
demos/demo_wf3.py | 1 -
nifty_core.py | 380 ++++++++++++++++++------------------
6 files changed, 459 insertions(+), 514 deletions(-)
diff --git a/demos/demo_excaliwir.py b/demos/demo_excaliwir.py
index 159c137cf..3aa48b892 100644
--- a/demos/demo_excaliwir.py
+++ b/demos/demo_excaliwir.py
@@ -28,270 +28,206 @@
.. /__/ /__/ /__/ /__/ \___/ \___ / demo
.. /______/
- NIFTY demo for (critical) Wiener filtering of Gaussian random signals.
+ NIFTY demo for (extended) critical Wiener filtering of Gaussian random signals.
"""
from __future__ import division
from nifty import *
-from scipy.sparse.linalg import LinearOperator as lo
-from scipy.sparse.linalg import cg
-note = notification()
+##=============================================================================
+class problem(object):
-##-----------------------------------------------------------------------------
+ def __init__(self, x_space, s2n=2, **kwargs):
+ """
+ Sets up a Wiener filter problem.
+
+ Parameters
+ ----------
+ x_space : space
+ Position space the signal lives in.
+ s2n : float, *optional*
+ Signal-to-noise ratio (default: 2).
-## spaces
-r1 = rg_space(512,1,zerocenter=False)
-r2 = rg_space(64,2)
-h = hp_space(16)
-g = gl_space(48)
-z = s_space = k = k_space = p = d_space = None
-
-## power spectrum (and more)
-power = kindex = rho = powerindex = powerundex = None
-
-## operators
-S = Sk = R = N = Nj = D = None
-
-## fields
-s = n = d = j = m = None
-
-## propagator class
-class propagator_operator(operator):
- """
- This is the information propagator from the Wiener filter formula.
- It is defined by its inverse. It is given the prior signal covariance S,
- the noise covariance N and the response R in para.
- """
- def _inverse_multiply(self,x):
- ## The inverse can be calculated directly
- S,N,R = self.para
- return S.inverse_times(x)+R.adjoint_times(N.inverse_times(R.times(x)))
-
- ## the inverse multiplication and multiplication with S modified to return 1d arrays
- _matvec = (lambda self,x: self.inverse_times(x).val.flatten())
- _precon = (lambda self,x: self.para[0].times(x).val.flatten())
-
- def _multiply(self,x):
"""
- the operator is defined by its inverse, so multiplication has to be
- done by inverting the inverse numerically using the conjugate gradient
- method from scipy
+ ## set signal space
+ self.z = x_space
+ ## set conjugate space
+ self.k = self.z.get_codomain()
+ self.k.set_power_indices(**kwargs)
+
+ ## set some power spectrum
+ self.power = (lambda k: 42 / (k + 1) ** 3)
+
+ ## define signal covariance
+ self.S = power_operator(self.k, spec=self.power, bare=True)
+ ## define projector to spectral bands
+ self.Sk = self.S.get_projection_operator()
+ ## generate signal
+ self.s = self.S.get_random_field(domain=self.z)
+
+ ## define response
+ self.R = response_operator(self.z, sigma=0.0, mask=1.0)
+ ## get data space
+ d_space = self.R.target
+
+ ## define noise covariance
+ self.N = diagonal_operator(d_space, diag=abs(s2n) * self.s.var(), bare=True)
+ ## define (plain) projector
+ self.Nj = projection_operator(d_space)
+ ## generate noise
+ n = self.N.get_random_field(domain=d_space)
+
+ ## compute data
+ self.d = self.R(self.s) + n
+
+ ## define information source
+ self.j = self.R.adjoint_times(self.N.inverse_times(self.d), target=self.k)
+ ## define information propagator
+ self.D = propagator_operator(S=self.S, N=self.N, R=self.R)
+
+ ## reserve map
+ self.m = None
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def solve(self, newspec=None):
"""
- A = lo(shape=tuple(self.dim()),matvec=self._matvec,dtype=self.domain.datatype) ## linear operator
- b = x.val.flatten()
- x_,info = cg(A,b,x0=None,tol=1.0E-5,maxiter=10*len(x),xtype=None,M=None,callback=None) ## conjugate gradient
- if(info==0):
- return x_
- else:
- note.cprint("NOTE: conjugate gradient failed.")
- return None
+ Solves the Wiener filter problem for a given power spectrum
+ reconstructing a signal estimate.
-##-----------------------------------------------------------------------------
+ Parameters
+ ----------
+ newspace : {scalar, list, array, field, function}, *optional*
+ Assumed power spectrum (default: k ** -2).
+ """
+ ## set (given) power spectrum
+ if(newspec is None):
+ newspec = np.r_[1, 1 / self.k.power_indices["kindex"][1:] ** 2] ## Laplacian
+ elif(newspec is False):
+ newspec = self.power ## assumed to be known
+ self.S.set_power(newspec, bare=True)
+ ## reconstruct map
+ self.m = self.D(self.j, W=self.S, tol=1E-3, note=False)
-##-----------------------------------------------------------------------------
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-def setup(space,s2n=3,nvar=None):
- """
- sets up the spaces, operators and fields
-
- Parameters
- ----------
- space : space
- the signal lives in `space`
- s2n : positive number, *optional*
- `s2n` is the signal to noise ratio (default: 3)
- nvar = positive number, *optional*
- the noise variance, `nvar` will be calculated according to
- `s2n` if not specified (default: None)
- """
- global z,s_space,k,k_space,p,d_space,power,kindex,rho,powerindex,powerundex,S,Sk,R,N,Nj,D,s,n,d,j,m
-
- ## signal space
- z = s_space = space
-
- ## conjugate space
- k = k_space = s_space.get_codomain()
- ## the power indices are calculated once and saved
- kindex,rho,powerindex,powerundex = k_space.get_power_indices()
-
- ## power spectrum
- power = (lambda kk: 42/(kk+1)**3)
-
- ## prior signal covariance operator (power operator)
- S = power_operator(k_space,spec=power)
- ## projection operator to the spectral bands
- Sk = S.get_projection_operator()
- ## the Gaussian random field generated from its power operator S
- s = S.get_random_field(domain=s_space,target=k_space)
-
- ## response
- R = response_operator(s_space,sigma=0,mask=1)
- ## data space
- p = d_space = R.target
-
- ## calculating the noise covariance
- if(nvar is None):
- svar = np.var(s.val) ## given unit response
- nvar = svar/s2n**2
- ## noise covariance operator
- N = diagonal_operator(d_space,diag=nvar,bare=True)
- ## Gaussian noise generated from its covariance N
- n = N.get_random_field(domain=d_space,target=d_space)
-
- ## data
- d = R(s)+n
+ def solve_critical(self, newspec=None, q=0, alpha=1, delta=1, epsilon=0):
+ """
+ Solves the (generalised) Wiener filter problem
+ reconstructing a signal estimate and a power spectrum.
+
+ Parameters
+ ----------
+ newspace : {scalar, list, array, field, function}, *optional*
+ Initial power spectrum (default: k ** -2).
+ q : {scalar, list, array}, *optional*
+ Spectral scale parameter of the assumed inverse-Gamme prior
+ (default: 0).
+ alpha : {scalar, list, array}, *optional*
+ Spectral shape parameter of the assumed inverse-Gamme prior
+ (default: 1).
+ delta : float, *optional*
+ First filter perception parameter (default: 1).
+ epsilon : float, *optional*
+ Second filter perception parameter (default: 0).
+
+ See Also
+ --------
+ infer_power
-##-----------------------------------------------------------------------------
+ """
+ ## set (initial) power spectrum
+ if(newspec is None):
+ newspec = np.r_[1, 1 / self.k.power_indices["kindex"][1:] ** 2] ## Laplacian
+ elif(newspec is False):
+ newspec = self.power ## assumed to be known
+ self.S.set_power(newspec, bare=True)
-##=============================================================================
+ ## pre-compute denominator
+ denominator = self.k.power_indices["rho"] + 2 * (alpha - 1 + abs(epsilon))
-def run(space=r1,s2n=3,nvar=None,**kwargs):
- """
- runs the demo of the generalised Wiener filter
-
- Parameters
- ----------
- space : space, *optional*
- `space` can be any space from nifty, that supports the plotting
- routine (default: r1 = rg_space(512,1,zerocenter=False))
- s2n : positive number, *optional*
- `s2n` is the signal to noise (default: 3)
- nvar : positive number, *optional*
- the noise variance, `nvar` will be calculated according to
- `s2n` if not specified (default: None)
- """
- global s_space,k_space,d_space,power,S,Sk,R,N,Nj,D,s,n,d,j,m
-
- ## setting up signal, noise, data and the operators S, N and R
- setup(space,s2n=s2n,nvar=nvar)
-
- ## information source
- j = R.adjoint_times(N.inverse_times(d))
-
- ## information propagator
- D = propagator_operator(s_space,sym=True,imp=True,para=[S,N,R])
-
- ## reconstructed map
- m = D(j)
- if(m is None):
- return None
-
- ## fields
- s.plot(title="signal",**kwargs)
-# n.cast_domain(s_space,newtarget=k_space)
-# n.plot(title="noise",**kwargs)
-# n.cast_domain(d_space,newtarget=d_space)
- d.cast_domain(s_space,newtarget=k_space)
- d.plot(title="data",vmin=np.min(s.val),vmax=np.max(s.val),**kwargs)
- d.cast_domain(d_space,newtarget=d_space)
- m.plot(title="reconstructed map",vmin=np.min(s.val),vmax=np.max(s.val),**kwargs)
-
- ## power spectrum
-# s.plot(title="power spectra",power=True,other=(m,power),mono=False)
-
- ## uncertainty
-# uncertainty = D.hat(bare=True,nrun=D.domain.dim()//4,target=k_space)
-# if(np.all(uncertainty.val>0)):
-# sqrt(uncertainty).plot(title="standard deviation",**kwargs)
+ ## iterate
+ iterating = True
+ while(iterating):
+
+ ## reconstruct map
+ self.m = self.D(self.j, W=self.S, tol=1E-3, note=False)
+ if(self.m is None):
+ break
+
+ ## reconstruct power spectrum
+ tr_B1 = self.Sk.pseudo_tr(self.m) ## == Sk(m).pseudo_dot(m)
+ tr_B2 = self.Sk.pseudo_tr(self.D, loop=True)
+
+ numerator = 2 * q + tr_B1 + abs(delta) * tr_B2 ## non-bare(!)
+ power = numerator / denominator
+
+ ## check convergence
+ dtau = np.log(power / self.S.get_power())
+ iterating = (np.max(np.abs(dtau)) > 2E-2)
+
+ ## update signal covariance
+ self.S.set_power(power, bare=False) ## auto-updates D
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def plot(self):
+ """
+ Produces plots.
+
+ """
+ ## plot signal
+ self.s.plot(title="signal")
+ ## plot data
+ try:
+ d_ = field(self.z, val=self.d.val, target=self.k)
+ d_.plot(title="data", vmin=self.s.min(), vmax=self.s.max())
+ except:
+ pass
+ ## plot map
+ if(self.m is None):
+ self.s.plot(power=True, mono=False, other=self.power)
+ else:
+ self.m.plot(title="reconstructed map", vmin=self.s.min(), vmax=self.s.max())
+ self.m.plot(power=True, mono=False, other=(self.power, self.S.get_power()))
##=============================================================================
##-----------------------------------------------------------------------------
-def run_critical(space=r2,s2n=3,nvar=None,q=1E-12,alpha=1,perception=[1,0],**kwargs):
- """
- runs the demo of the critical generalised Wiener filter
-
- Parameters
- ----------
- space : space, *optional*
- `space` can be any space from nifty, that supports the plotting
- routine (default: r2 = rg_space(64,2))
- s2n : positive number, *optional*
- `s2n` is the signal to noise (default: 3)
- nvar : positive number, *optional*
- the noise variance, `nvar` will be calculated according to
- `s2n` if not specified (default: None)
- q : positive number, *optional*
- `q` is the minimal power on all scales (default: 1E-12)
- alpha : a number >= 1, *optional*
- `alpha` = 1 means Jeffreys prior for the power spectrum (default: 1)
- perception : array of shape (2,1), *optional*
- perception[0] is delta, perception[1] is epsilon. They are tuning
- factors for the filter (default: [1,0])
-
- See Also
- --------
- infer_power
-
- """
- global s_space,k_space,d_space,power,rho,S,Sk,R,N,Nj,D,s,n,d,j,m
-
- ## setting up signal, noise, data and the operators S, N and R
- setup(space,s2n=s2n,nvar=nvar)
- if(perception[1] is None):
- perception[1] = rho/2*(perception[0]-1)
-
- ## information source
- j = R.adjoint_times(N.inverse_times(d))
-
- ## unknown power spectrum, the power operator is given an initial value
- S.set_power(42,bare=True) ## The answer is 42. I double checked.
- ## the power spectrum is drawn from the first guess power operator
- pk = S.get_power(bare=False) ## non-bare(!)
-
- ## information propagator
- D = propagator_operator(s_space,sym=True,imp=True,para=[S,N,R])
-
- ## iterative reconstruction of the power spectrum and the map
- iteration = 0
- while(True):
- iteration += 1
-
- ## the Wiener filter reconstruction using the current power spectrum
- m = D(j)
- if(m is None):
- return None
-
- ## measuring a new estimated power spectrum from the current reconstruction
- b1 = Sk.pseudo_tr(m) ## == Sk(m).pseudo_dot(m), but faster
- b2 = Sk.pseudo_tr(D,nrun=np.sqrt(Sk.domain.dim())//4) ## probing of the partial traces of D
- pk_new = (2*q+b1+perception[0]*b2)/(rho+2*(alpha-1+perception[1])) ## non-bare(!)
- pk_new = smooth_power(pk_new,domain=k_space,mode="2s",exclude=min(8,len(rho))) ## smoothing
- ## the power operator is given the new spectrum
- S.set_power(pk_new,bare=False) ## auto-updates D
-
- ## check convergence
- log_change = np.max(np.abs(log(pk_new/pk)))
- if(log_change<=0.01):
- note.cprint("NOTE: desired accuracy reached in iteration %u."%(iteration))
- break
- else:
- note.cprint("NOTE: log-change == %4.3f ( > 1%% ) in iteration %u."%(log_change,iteration))
- pk = np.copy(pk_new)
-
- ## fields
- s.plot(title="signal",**kwargs)
-# n.cast_domain(s_space,newtarget=k_space)
-# n.plot(title="noise",**kwargs)
-# n.cast_domain(d_space,newtarget=d_space)
- d.cast_domain(s_space,newtarget=k_space)
- d.plot(title="data",vmin=np.min(s.val),vmax=np.max(s.val),**kwargs)
- d.cast_domain(d_space,newtarget=d_space)
- m.plot(title="reconstructed map",vmin=np.min(s.val),vmax=np.max(s.val),**kwargs)
-
- ## power spectrum
- s.plot(title="power spectra",power=True,other=(S.get_power(),power),mono=False)
-
- ## uncertainty
-# uncertainty = D.hat(bare=True,nrun=D.domain.dim()//4,target=k_space)
-# if(np.all(uncertainty.val>0)):
-# sqrt(uncertainty).plot(title="standard deviation")
+if(__name__=="__main__"):
+# pl.close("all")
+
+ ## define signal space
+ x_space = rg_space(128)
+
+ ## setup problem
+ p = problem(x_space, log=True)
+ ## solve problem given some power spectrum
+ p.solve()
+ ## solve problem
+ p.solve_critical()
+
+ p.plot()
+
+ ## retrieve objects
+ k_space = p.k
+ power = p.power
+ S = p.S
+ Sk = p.Sk
+ s = p.s
+ R = p.R
+ d_space = p.R.target
+ N = p.N
+ Nj = p.Nj
+ d = p.d
+ j = p.j
+ D = p.D
+ m = p.m
##-----------------------------------------------------------------------------
diff --git a/demos/demo_faraday.py b/demos/demo_faraday.py
index b89b8a954..ee8037853 100644
--- a/demos/demo_faraday.py
+++ b/demos/demo_faraday.py
@@ -47,24 +47,11 @@ about.infos.off()
##-----------------------------------------------------------------------------
-## spaces
-h = hp_space(128)
-l = lm_space(383)
-g = gl_space(384) ## nlon == 767
-g_ = gl_space(384, nlon=768)
-r = rg_space([768, 384], dist=[1/360, 1/180])
-r_ = rg_space([256, 128], dist=[1/120, 1/60])
-
-## map
-m = field(h, val=np.load("demo_faraday_map.npy"))
-
-## projection operator
-Sk = None
+# (global) Faraday map
+m = field(hp_space(128), val=np.load("demo_faraday_map.npy"))
##-----------------------------------------------------------------------------
-
-
##=============================================================================
def run(projection=False, power=False):
@@ -74,47 +61,71 @@ def run(projection=False, power=False):
Parameters
----------
projection : bool, *optional*
- Whether to additionaly include projections in the demo or not. If
- ``projection == True`` the projection operator `Sk` will be
- defined. (default: False)
+ Whether to additionaly show projections or not (default: False).
power : bool, *optional*
- Whether to additionaly show power spectra in the demo or not.
- (default: False)
+ Whether to additionaly show power spectra or not (default: False).
"""
- global Sk
- ## start in hp_space
+ nicely = {"vmin":-4, "vmax":4, "cmap":ncmap.fm()}
+
+ # (a) representation on HEALPix grid
m0 = m
+ m0.plot(title=r"$m$ on a HEALPix grid", **nicely)
+ nicely.update({"cmap":ncmap.fm()}) # healpy bug workaround
- ## transform to lm_space
- m1 = m0.transform(l)
+ # (b) representation in spherical harmonics
+ k_space = m0.target # == lm_space(383, mmax=383)
+ m1 = m0.transform(k_space) # == m.transform()
+# m1.plot(title=r"$m$ in spherical harmonics")
+
+ if(power):
+ m1.plot(title=r"angular power spectrum of $m$", vmin=1E-2, vmax=1E+1, mono=False)
if(projection):
- ## define projection operator
- Sk = projection_operator(l)
- ## project quadrupole
+ # define projection operator
+ Sk = projection_operator(m1.domain)
+ # project quadrupole
m2 = Sk(m0, band=2)
+ m2.plot(title=r"angular quadrupole of $m$ on a HEALPix grid", **nicely)
- ## transform to gl_space
- m3 = m1.transform(g)
+ # (c) representation on Gauss-Legendre grid
+ y_space = m.target.get_codomain(coname="gl") # == gl_space(384, nlon=767)
+ m3 = m1.transform(y_space) # == m0.transform().transform(y_space)
+ m3.plot(title=r"$m$ on a spherical Gauss-Legendre grid", **nicely)
- ## transform to rg_space
- m4 = m1.transform(g_) ## auxiliary gl_space
- m4.cast_domain(r) ## rg_space cast
- m4.set_val(np.roll(m4.val[::-1, ::-1], g.nlon()//2, axis=1)) ## rearrange
- if(power):
- ## restrict to central window
- m5 = field(r_, val=m4[128:256, 256:512]).transform()
+ if(projection):
+ m4 = Sk(m3, band=2)
+ m4.plot(title=r"angular quadrupole of $m$ on a Gauss-Legendre grid", **nicely)
+
+ # (d) representation on regular grid
+ y_space = gl_space(384, nlon=768) # auxiliary gl_space
+ z_space = rg_space([768, 384], dist=[1/360, 1/180])
+ m5 = m1.transform(y_space)
+ m5.cast_domain(z_space)
+ m5.set_val(np.roll(m5.val[::-1, ::-1], y_space.nlon()//2, axis=1)) # rearrange value array
+ m5.plot(title=r"$m$ on a regular 2D grid", **nicely)
- ## plots
- m0.plot(title=r"$m$ on a HEALPix grid", vmin=-4, vmax=4, cmap=ncmap.fm())
if(power):
- m1.plot(title=r"angular power spectrum of $m$", vmin=1E-2, vmax=1E+1, mono=False)
+ m5.target.set_power_indices(log=False)
+ m5.plot(power=True, title=r"Fourier power spectrum of $m$", vmin=1E-3, vmax=1E+0, mono=False)
if(projection):
- m2.plot(title=r"quadrupole of $m$ on a HEALPix grid", vmin=-4, vmax=4, cmap=ncmap.fm())
- m3.plot(title=r"$m$ on a spherical Gauss-Legendre grid", vmin=-4, vmax=4, cmap=ncmap.fm())
- m4.plot(title=r"$m$ on a regular 2D grid", vmin=-4, vmax=4, cmap=ncmap.fm())
- if(power):
- m5.plot(title=r"(restricted, binned) Fourier power spectrum of $m$", vmin=1E-3, vmax=1E+0, mono=False, log=True)
+ m5.target.set_power_indices(log=False)
+ # define projection operator
+ Sk = projection_operator(m5.target)
+ # project quadrupole
+ m6 = Sk(m5, band=2)
+ m6.plot(title=r"Fourier quadrupole of $m$ on a regular 2D grid", **nicely)
##=============================================================================
+##-----------------------------------------------------------------------------
+
+if(__name__=="__main__"):
+# pl.close("all")
+
+ # run demo
+ run(projection=False, power=False)
+ # define projection operator
+ Sk = projection_operator(m.target)
+
+##-----------------------------------------------------------------------------
+
diff --git a/demos/demo_wf1.py b/demos/demo_wf1.py
index 580f25ac8..ac96c0b0b 100644
--- a/demos/demo_wf1.py
+++ b/demos/demo_wf1.py
@@ -32,36 +32,39 @@
"""
from __future__ import division
-from nifty import * # version 0.6.0
+from nifty import * # version 0.8.0
# some signal space; e.g., a two-dimensional regular grid
-x_space = rg_space([256, 256]) # define signal space
+x_space = rg_space([128, 128]) # define signal space
+#x_space = rg_space(512)
+#x_space = hp_space(32)
+#x_space = gl_space(96)
-k_space = x_space.get_codomain() # get conjugate space
+k_space = x_space.get_codomain() # get conjugate space
# some power spectrum
power = (lambda k: 42 / (k + 1) ** 3)
-S = power_operator(k_space, spec=power) # define signal covariance
-s = S.get_random_field(domain=x_space) # generate signal
+S = power_operator(k_space, spec=power) # define signal covariance
+s = S.get_random_field(domain=x_space) # generate signal
-R = response_operator(x_space, sigma=0.0, mask=1.0, assign=None) # define response
-d_space = R.target # get data space
+R = response_operator(x_space, sigma=0.0, mask=1.0, assign=None) # define response
+d_space = R.target # get data space
-# some noise variance; e.g., 100
-N = diagonal_operator(d_space, diag=100, bare=True) # define noise covariance
-n = N.get_random_field(domain=d_space) # generate noise
+# some noise variance; e.g., signal-to-noise ratio of 1
+N = diagonal_operator(d_space, diag=s.var(), bare=True) # define noise covariance
+n = N.get_random_field(domain=d_space) # generate noise
-d = R(s) + n # compute data
+d = R(s) + n # compute data
-j = R.adjoint_times(N.inverse_times(d)) # define information source
-D = propagator_operator(S=S, N=N, R=R) # define information propagator
+j = R.adjoint_times(N.inverse_times(d)) # define information source
+D = propagator_operator(S=S, N=N, R=R) # define information propagator
-m = D(j, tol=1E-4, note=True) # reconstruct map
+m = D(j, W=S, tol=1E-3, note=True) # reconstruct map
-s.plot(title="signal") # plot signal
+s.plot(title="signal") # plot signal
d_ = field(x_space, val=d.val, target=k_space)
-d_.plot(title="data", vmin=s.val.min(), vmax=s.val.max()) # plot data
-m.plot(title="reconstructed map", vmin=s.val.min(), vmax=s.val.max()) # plot map
+d_.plot(title="data", vmin=s.min(), vmax=s.max()) # plot data
+m.plot(title="reconstructed map", vmin=s.min(), vmax=s.max()) # plot map
diff --git a/demos/demo_wf2.py b/demos/demo_wf2.py
index 8d0b03500..2d9f81983 100644
--- a/demos/demo_wf2.py
+++ b/demos/demo_wf2.py
@@ -32,31 +32,31 @@
"""
from __future__ import division
-from nifty import * # version 0.6.0
+from nifty import * # version 0.8.0
# some signal space; e.g., a two-dimensional regular grid
-x_space = rg_space([256, 256]) # define signal space
+x_space = rg_space([128, 128]) # define signal space
-k_space = x_space.get_codomain() # get conjugate space
+k_space = x_space.get_codomain() # get conjugate space
# some power spectrum
power = (lambda k: 42 / (k + 1) ** 3)
-S = power_operator(k_space, spec=power) # define signal covariance
-s = S.get_random_field(domain=x_space) # generate signal
+S = power_operator(k_space, spec=power) # define signal covariance
+s = S.get_random_field(domain=x_space) # generate signal
-R = response_operator(x_space, sigma=0.0, mask=1.0, assign=None) # define response
-d_space = R.target # get data space
+R = response_operator(x_space, sigma=0.0, mask=1.0, assign=None) # define response
+d_space = R.target # get data space
-# some noise variance; e.g., 100
-N = diagonal_operator(d_space, diag=100, bare=True) # define noise covariance
-n = N.get_random_field(domain=d_space) # generate noise
+# some noise variance; e.g., signal-to-noise ratio of 1
+N = diagonal_operator(d_space, diag=s.var(), bare=True) # define noise covariance
+n = N.get_random_field(domain=d_space) # generate noise
-d = R(s) + n # compute data
+d = R(s) + n # compute data
-j = R.adjoint_times(N.inverse_times(d)) # define information source
-D = propagator_operator(S=S, N=N, R=R) # define information propagator
+j = R.adjoint_times(N.inverse_times(d)) # define information source
+D = propagator_operator(S=S, N=N, R=R) # define information propagator
def eggs(x):
@@ -65,16 +65,16 @@ def eggs(x):
"""
DIx = D.inverse_times(x)
- H = 0.5 * DIx.dot(x) - j.dot(x) # compute information Hamiltonian
- g = DIx - j # compute its gradient
+ H = 0.5 * DIx.dot(x) - j.dot(x) # compute information Hamiltonian
+ g = DIx - j # compute its gradient
return H, g
-m = field(x_space, target=k_space) # reconstruct map
-m,convergence = steepest_descent(eggs=eggs, note=True)(m, tol=1E-4, clevel=3)
+m = field(x_space, target=k_space) # reconstruct map
+m, convergence = steepest_descent(eggs=eggs, note=True)(m, tol=1E-3, clevel=3)
-s.plot(title="signal") # plot signal
+s.plot(title="signal") # plot signal
d_ = field(x_space, val=d.val, target=k_space)
-d_.plot(title="data", vmin=s.val.min(), vmax=s.val.max()) # plot data
-m.plot(title="reconstructed map", vmin=s.val.min(), vmax=s.val.max()) # plot map
+d_.plot(title="data", vmin=s.min(), vmax=s.max()) # plot data
+m.plot(title="reconstructed map", vmin=s.min(), vmax=s.max()) # plot map
diff --git a/demos/demo_wf3.py b/demos/demo_wf3.py
index 5b6038056..b9215591e 100644
--- a/demos/demo_wf3.py
+++ b/demos/demo_wf3.py
@@ -33,7 +33,6 @@
"""
from __future__ import division
from nifty import * # version 0.7.0
-from nifty.nifty_explicit import *
# some signal space; e.g., a one-dimensional regular grid
diff --git a/nifty_core.py b/nifty_core.py
index 1312131be..e39ce06c3 100644
--- a/nifty_core.py
+++ b/nifty_core.py
@@ -1100,73 +1100,6 @@ class space(object):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-# def get_power_index(self,irreducible=False): ## TODO: remove in future version
-# """
-# **DEPRECATED** Provides the indexing array of the power spectrum.
-#
-# Provides either an array giving for each component of a field the
-# corresponding index of a power spectrum (if ``irreducible==False``)
-# or two arrays containing the scales of the modes and the numbers of
-# modes with this scale (if ``irreducible==True``).
-#
-# Parameters
-# ----------
-# irreducible : bool, *optional*
-# Whether to return two arrays containing the scales and
-# corresponding number of represented modes (if True) or the
-# indexing array (if False) (default: False).
-#
-# Returns
-# -------
-# kindex : numpy.ndarray
-# Scale of each band, returned only if ``irreducible==True``.
-# rho : numpy.ndarray
-# Number of modes per scale represented in the discretization,
-# returned only if ``irreducible==True``.
-# pindex : numpy.ndarray
-# Indexing array giving the power spectrum index for each
-# represented mode, returned only if ``irreducible==False``.
-#
-# Notes
-# -----
-# The indexing array is of the same shape as a field living in this
-# space and contains the indices of the associated bands.
-# kindex and rho are each one-dimensional arrays.
-# """
-# about.warnings.cprint("WARNING: 'get_power_index' is deprecated.")
-# raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'get_power_index'."))
-
- def get_power_undex(self,pindex=None): ## TODO: remove in future version
- """
- **DEPRECATED** Provides the Unindexing array for an indexed power spectrum.
-
- Parameters
- ----------
- pindex : numpy.ndarray, *optional*
- Indexing array giving the power spectrum index for each
- represented mode.
-
- Returns
- -------
- pundex : numpy.ndarray
- Unindexing array undoing power indexing.
-
- Notes
- -----
- Indexing with the unindexing array undoes the indexing with the
- indexing array; i.e., ``power == power[pindex].flatten()[pundex]``.
-
- See Also
- --------
- get_power_index
-
- """
- about.warnings.cprint("WARNING: 'get_power_undex' is deprecated.")
- if(pindex is None):
- pindex = self.get_power_index(irreducible=False)
-# return list(np.unravel_index(np.unique(pindex,return_index=True,return_inverse=False)[1],pindex.shape,order='C')) ## < version 0.4
- return np.unique(pindex,return_index=True,return_inverse=False)[1]
-
def set_power_indices(self,**kwargs):
"""
Sets the (un)indexing objects for spectral indexing internally.
@@ -1985,14 +1918,6 @@ class point_space(space):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-# def get_power_index(self,irreducible=False): ## TODO: remove in future version
-# """
-# **DEPRECATED** Raises an error since the power spectrum is
-# ill-defined for point spaces.
-# """
-# about.warnings.cprint("WARNING: 'get_power_index' is deprecated.")
-# raise AttributeError(about._errors.cstring("ERROR: power spectra ill-defined for (unstructured) point spaces."))
-
def set_power_indices(self,**kwargs):
"""
Raises
@@ -2502,45 +2427,6 @@ class rg_space(space):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-# def get_power_index(self,irreducible=False): ## TODO: remove in future version
-# """
-# **DEPRECATED** Provides the indexing array of the power spectrum.
-#
-# Provides either an array giving for each component of a field the
-# corresponding index of a power spectrum (if ``irreducible==False``)
-# or two arrays containing the scales of the modes and the numbers of
-# modes with this scale (if ``irreducible==True``).
-#
-# Parameters
-# ----------
-# irreducible : bool, *optional*
-# Whether to return two arrays containing the scales and
-# corresponding number of represented modes (if True) or the
-# indexing array (if False) (default: False).
-#
-# Returns
-# -------
-# kindex : numpy.ndarray
-# Scale of each band, returned only if ``irreducible==True``.
-# rho : numpy.ndarray
-# Number of modes per scale represented in the discretization,
-# returned only if ``irreducible==True``.
-# pindex : numpy.ndarray
-# Indexing array giving the power spectrum index for each
-# represented mode, returned only if ``irreducible==False``.
-#
-# Notes
-# -----
-# The indexing array is of the same shape as a field living in this
-# space and contains the indices of the associated bands.
-# kindex and rho are each one-dimensional arrays.
-# """
-# about.warnings.cprint("WARNING: 'get_power_index' is deprecated.")
-# if(self.fourier):
-# return gp.get_power_index(self.para[:(np.size(self.para)-1)//2],self.vol,self.para[-((np.size(self.para)-1)//2):].astype(np.bool),irred=irreducible,fourier=self.fourier) ## nontrivial
-# else:
-# raise AttributeError(about._errors.cstring("ERROR: power spectra indexing ill-defined."))
-
def set_power_indices(self,**kwargs):
"""
Sets the (un)indexing objects for spectral indexing internally.
@@ -3613,46 +3499,6 @@ class lm_space(space):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
- def get_power_index(self,irreducible=False): ## TODO: remove in future version
- """
- **DEPRECATED** Provides the indexing array of the power spectrum.
-
- Provides either an array giving for each component of a field the
- corresponding index of a power spectrum (if ``irreducible==False``)
- or two arrays containing the scales of the modes and the numbers of
- modes with this scale (if ``irreducible==True``).
-
- Parameters
- ----------
- irreducible : bool, *optional*
- Whether to return two arrays containing the scales and
- corresponding number of represented modes (if True) or the
- indexing array (if False) (default: False).
-
- Returns
- -------
- kindex : numpy.ndarray
- Scale of each band, returned only if ``irreducible==True``.
- rho : numpy.ndarray
- Number of modes per scale represented in the discretization,
- returned only if ``irreducible==True``.
- pindex : numpy.ndarray
- Indexing array giving the power spectrum index for each
- represented mode, returned only if ``irreducible==False``.
-
- Notes
- -----
- The indexing array is of the same shape as a field living in this
- space and contains the indices of the associated bands.
- kindex and rho are each one-dimensional arrays.
- """
- about.warnings.cprint("WARNING: 'get_power_index' is deprecated.")
- if(irreducible):
- ind = np.arange(self.para[0]+1)
- return ind,2*ind+1
- else:
- return hp.Alm.getlm(self.para[0],i=None)[0] ## l of (l,m)
-
def set_power_indices(self,**kwargs):
"""
Sets the (un)indexing objects for spectral indexing internally.
@@ -4451,15 +4297,6 @@ class gl_space(space):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-# def get_power_index(self,irreducible=False): ## TODO: remove in future version
-# """
-# **DEPRECATED** Raises an error since the power spectrum for a field on the sphere
-# is defined via the spherical harmonics components and not its
-# position-space representation.
-# """
-# about.warnings.cprint("WARNING: 'get_power_index' is deprecated.")
-# raise AttributeError(about._errors.cstring("ERROR: power spectra indexing ill-defined."))
-
def set_power_indices(self,**kwargs):
"""
Raises
@@ -5094,15 +4931,6 @@ class hp_space(space):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-# def get_power_index(self,irreducible=False): ## TODO: remove in future version
-# """
-# **DEPRECATED** Raises an error since the power spectrum for a field on the sphere
-# is defined via the spherical harmonics components and not its
-# position-space representation.
-# """
-# about.warnings.cprint("WARNING: 'get_power_index' is deprecated.")
-# raise AttributeError(about._errors.cstring("ERROR: power spectra indexing ill-defined."))
-
def set_power_indices(self,**kwargs):
"""
Raises
@@ -5636,14 +5464,6 @@ class nested_space(space):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-# def get_power_index(self,irreducible=False): ## TODO: remove in future version
-# """
-# **DEPRECATED** Raises an error since there is no canonical
-# definition for the power spectrum on a generic product space.
-# """
-# about.warnings.cprint("WARNING: 'get_power_index' is deprecated.")
-# raise AttributeError(about._errors.cstring("ERROR: power spectra indexing ill-defined."))
-
def set_power_indices(self,**kwargs):
"""
Raises
@@ -6332,22 +6152,21 @@ class field(object):
def dim(self,split=False):
"""
- Computes the (array) dimension of the underlying space.
+ Computes the (array) dimension of the underlying space.
- Parameters
- ----------
- split : bool
- Sets the output to be either split up per axis or
- in form of total number of field entries in all
- dimensions (default=False)
+ Parameters
+ ----------
+ split : bool
+ Sets the output to be either split up per axis or
+ in form of total number of field entries in all
+ dimensions (default=False)
- Returns
- -------
- dim : {scalar, ndarray}
- Dimension of space.
+ Returns
+ -------
+ dim : {scalar, ndarray}
+ Dimension of space.
"""
-
return self.domain.dim(split=split)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -6893,6 +6712,183 @@ class field(object):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ def min(self,ignore=False,**kwargs):
+ """
+ Returns the minimum of the field values.
+
+ Parameters
+ ----------
+ ignore : bool
+ Whether to ignore NANs or not (default: False).
+
+ Returns
+ -------
+ amin : {scalar, ndarray}
+ Minimum field value.
+
+ See Also
+ --------
+ np.amin, np.nanmin
+
+ """
+ if(ignore):
+ return np.nanmin(self.val,**kwargs)
+ else:
+ return np.amin(self.val,**kwargs)
+
+ def max(self,ignore=False,**kwargs):
+ """
+ Returns the maximum of the field values.
+
+ Parameters
+ ----------
+ ignore : bool
+ Whether to ignore NANs or not (default: False).
+
+ Returns
+ -------
+ amax : {scalar, ndarray}
+ Maximum field value.
+
+ See Also
+ --------
+ np.amax, np.nanmax
+
+ """
+ if(ignore):
+ return np.nanmax(self.val,**kwargs)
+ else:
+ return np.amax(self.val,**kwargs)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def med(self,**kwargs):
+ """
+ Returns the median of the field values.
+
+ Returns
+ -------
+ med : scalar
+ Median field value.
+
+ See Also
+ --------
+ np.median
+
+ """
+ return np.median(self.val,**kwargs)
+
+ def mean(self,**kwargs):
+ """
+ Returns the mean of the field values.
+
+ Returns
+ -------
+ mean : scalar
+ Mean field value.
+
+ See Also
+ --------
+ np.mean
+
+ """
+ return np.mean(self.val,**kwargs)
+
+ def std(self,**kwargs):
+ """
+ Returns the standard deviation of the field values.
+
+ Returns
+ -------
+ std : scalar
+ Standard deviation of the field values.
+
+ See Also
+ --------
+ np.std
+
+ """
+ return np.std(self.val,**kwargs)
+
+ def var(self,**kwargs):
+ """
+ Returns the variance of the field values.
+
+ Returns
+ -------
+ var : scalar
+ Variance of the field values.
+
+ See Also
+ --------
+ np.var
+
+ """
+ return np.var(self.val,**kwargs)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+ def argmin(self,split=True,**kwargs):
+ """
+ Returns the index of the minimum field value.
+
+ Parameters
+ ----------
+ split : bool
+ Whether to split (unravel) the flat index or not; does not
+ apply to multiple indices along some axis (default: True).
+
+ Returns
+ -------
+ ind : {integer, tuple, array}
+ Index of the minimum field value being an integer for
+ one-dimensional fields, a tuple for multi-dimensional fields,
+ and an array in case minima along some axis are requested.
+
+ See Also
+ --------
+ np.argmax, np.argmin
+
+ """
+ ind = np.argmin(self.val,**kwargs)
+ if(split)and(np.isscalar(ind)):
+ ind = np.unravel_index(ind,self.val.shape,order='C')
+ if(len(ind)==1):
+ return ind[0]
+ return ind
+
+ def argmax(self,split=True,**kwargs):
+ """
+ Returns the index of the maximum field value.
+
+ Parameters
+ ----------
+ split : bool
+ Whether to split (unravel) the flat index or not; does not
+ apply to multiple indices along some axis (default: True).
+
+ Returns
+ -------
+ ind : {integer, tuple, array}
+ Index of the maximum field value being an integer for
+ one-dimensional fields, a tuple for multi-dimensional fields,
+ and an array in case maxima along some axis are requested.
+
+ See Also
+ --------
+ np.argmax, np.argmin
+
+ """
+ ind = np.argmax(self.val,**kwargs)
+ if(split)and(np.isscalar(ind)):
+ ind = np.unravel_index(ind,self.val.shape,order='C')
+ if(len(ind)==1):
+ return ind[0]
+ return ind
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
def __len__(self):
return int(self.dim(split=True)[0])
--
GitLab