Commit 4cc0f412 authored by M Selig's avatar M Selig
Browse files

Merge pull request #9 from mselig/develop

version update
parents d4d18b1c 40a174f8
......@@ -99,7 +99,7 @@ Requirements
Download
........
The latest release is tagged **v0.7.0** and is available as a source package
The latest release is tagged **v0.8.0** and is available as a source package
at `<https://github.com/mselig/nifty/tags>`_. The current version can be
obtained by cloning the repository::
......
......@@ -28,270 +28,207 @@
.. /__/ /__/ /__/ /__/ \___/ \___ / 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()
##-----------------------------------------------------------------------------
## 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
class problem(object):
## 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 __init__(self, x_space, s2n=2, **kwargs):
"""
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)))
Sets up a Wiener filter problem.
## 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())
Parameters
----------
x_space : space
Position space the signal lives in.
s2n : float, *optional*
Signal-to-noise ratio (default: 2).
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
##-----------------------------------------------------------------------------
##-----------------------------------------------------------------------------
def setup(space,s2n=3,nvar=None):
"""
sets up the spaces, operators and fields
Solves the Wiener filter problem for a given power spectrum
reconstructing a signal estimate.
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)
newspace : {scalar, list, array, field, function}, *optional*
Assumed power spectrum (default: k ** -2).
"""
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
## 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 run(space=r1,s2n=3,nvar=None,**kwargs):
def solve_critical(self, newspec=None, q=0, alpha=1, delta=1, epsilon=0):
"""
runs the demo of the generalised Wiener filter
Solves the (generalised) Wiener filter problem
reconstructing a signal estimate and a power spectrum.
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
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).
## 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])
See Also
--------
infer_power
## reconstructed map
m = D(j)
if(m is None):
return None
"""
## 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))
## 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
## 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)
## 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)
## power spectrum
# s.plot(title="power spectra",power=True,other=(m,power),mono=False)
numerator = 2 * q + tr_B1 + abs(delta) * tr_B2 ## non-bare(!)
power = numerator / denominator
## 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)
## check convergence
dtau = log(power / self.S.get_power(), base=self.S.get_power())
iterating = (np.max(np.abs(dtau)) > 2E-2)
print max(np.abs(dtau))
##=============================================================================
## update signal covariance
self.S.set_power(power, bare=False) ## auto-updates D
##-----------------------------------------------------------------------------
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def run_critical(space=r2,s2n=3,nvar=None,q=1E-12,alpha=1,perception=[1,0],**kwargs):
def plot(self):
"""
runs the demo of the critical generalised Wiener filter
Produces plots.
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])
"""
## 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()))
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
##-----------------------------------------------------------------------------
......@@ -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)
##-----------------------------------------------------------------------------
......@@ -32,11 +32,14 @@
"""
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
......@@ -49,8 +52,8 @@ s = S.get_random_field(domain=x_space) # generate
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
# 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
......@@ -58,10 +61,10 @@ d = R(s) + n # compute
j = R.adjoint_times(N.inverse_times(d)) # define information source
D = propagator_operator(S=S, N=N, R=R) # define information propagator