diff --git a/.gitignore b/.gitignore
index ea98e22b728117941248f28ab64cf202827b1219..81fc51b2496005396047398c8364e2095e70309c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -15,4 +15,6 @@ build
demos/*
!demos/demo_faraday.py
!demos/demo_faraday_map.npy
-!demos/demo_excaliwir.py
\ No newline at end of file
+!demos/demo_excaliwir.py
+!demos/demo_wf1.py
+!demos/demo_wf2.py
\ No newline at end of file
diff --git a/demos/demo_wf1.py b/demos/demo_wf1.py
new file mode 100644
index 0000000000000000000000000000000000000000..e4035cb2fd365767a4a6c1db1b54cea69c8dece9
--- /dev/null
+++ b/demos/demo_wf1.py
@@ -0,0 +1,73 @@
+## 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/>.
+
+"""
+ .. __ ____ __
+ .. /__/ / _/ / /_
+ .. __ ___ __ / /_ / _/ __ __
+ .. / _ | / / / _/ / / / / / /
+ .. / / / / / / / / / /_ / /_/ /
+ .. /__/ /__/ /__/ /__/ \___/ \___ / demo
+ .. /______/
+
+ NIFTY demo applying a Wiener filter using conjugate gradient.
+
+"""
+from __future__ import division
+from nifty import * # version 0.6.0
+from nifty.nifty_tools import *
+
+
+
+# some signal space; e.g., a two-dimensional regular grid
+x_space = rg_space([256, 256]) # define signal 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
+
+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
+
+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
+
+m = D(j, tol=1E-4, note=True) # reconstruct map
+
+
+
+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
+
diff --git a/demos/demo_wf2.py b/demos/demo_wf2.py
new file mode 100644
index 0000000000000000000000000000000000000000..6c3def9a8f37e5e6c7fd485aec3bf0869d17f7d8
--- /dev/null
+++ b/demos/demo_wf2.py
@@ -0,0 +1,86 @@
+## 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/>.
+
+"""
+ .. __ ____ __
+ .. /__/ / _/ / /_
+ .. __ ___ __ / /_ / _/ __ __
+ .. / _ | / / / _/ / / / / / /
+ .. / / / / / / / / / /_ / /_/ /
+ .. /__/ /__/ /__/ /__/ \___/ \___ / demo
+ .. /______/
+
+ NIFTY demo applying a Wiener filter using steepest descent.
+
+"""
+from __future__ import division
+from nifty import * # version 0.6.0
+from nifty.nifty_tools import *
+
+
+
+# some signal space; e.g., a two-dimensional regular grid
+x_space = rg_space([256, 256]) # define signal 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
+
+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
+
+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
+
+
+def eggs(x):
+ """
+ Calculation of the information Hamiltonian and its gradient.
+
+ """
+ Dx = D.inverse_times(x)
+ H = 0.5 * Dx.dot(x) - j.dot(x) # compute information Hamiltonian
+ g = Dx - 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)
+
+
+
+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
+
diff --git a/nifty_core.py b/nifty_core.py
index 004d4e357a7779342a2c28ac168cb271476567ba..cbaf09919ffb2567299a25bf208cb947f2db1fbd 100644
--- a/nifty_core.py
+++ b/nifty_core.py
@@ -2909,6 +2909,11 @@ class rg_space(space):
dot = np.dot(np.conjugate(x.flatten(order='C')),y.flatten(order='C'),out=None)
if(np.isreal(dot)):
return np.real(dot)
+ elif(self.para[(np.size(self.para)-1)//2]!=2):
+ ## check imaginary part
+ if(dot.imag>self.epsilon**2*dot.real):
+ about.warnings.cprint("WARNING: discarding considerable imaginary part.")
+ return np.real(dot)
else:
return dot
@@ -8138,7 +8143,7 @@ class operator(object):
domain = self.domain
diag = diagonal_probing(self,function=self.times,domain=domain,target=target,random=random,ncpu=ncpu,nrun=nrun,nper=nper,var=var,save=save,path=path,prefix=prefix,**kwargs)(loop=loop)
if(diag is None):
- about.warnings.cprint("WARNING: 'NoneType' forwarded.")
+# about.warnings.cprint("WARNING: forwarding 'NoneType'.")
return None
## weight if ...
elif(not domain.discrete)and(bare):
@@ -8218,7 +8223,7 @@ class operator(object):
domain = self.target
diag = diagonal_probing(self,function=self.inverse_times,domain=domain,target=target,random=random,ncpu=ncpu,nrun=nrun,nper=nper,var=var,save=save,path=path,prefix=prefix,**kwargs)(loop=loop)
if(diag is None):
- about.warnings.cprint("WARNING: 'NoneType' forwarded.")
+# about.warnings.cprint("WARNING: forwarding 'NoneType'.")
return None
## weight if ...
elif(not domain.discrete)and(bare):
@@ -8319,6 +8324,7 @@ class operator(object):
domain = self.domain
diag = self.diag(bare=bare,domain=domain,target=target,var=False,**kwargs)
if(diag is None):
+ about.warnings.cprint("WARNING: forwarding 'NoneType'.")
return None
return field(domain,val=diag,target=target)
@@ -8386,6 +8392,7 @@ class operator(object):
domain = self.target
diag = self.inverse_diag(bare=bare,domain=domain,target=target,var=False,**kwargs)
if(diag is None):
+ about.warnings.cprint("WARNING: forwarding 'NoneType'.")
return None
return field(domain,val=diag,target=target)
@@ -8448,6 +8455,7 @@ class operator(object):
domain = self.domain
diag = self.diag(bare=False,domain=domain,var=False,**kwargs)
if(diag is None):
+ about.warnings.cprint("WARNING: forwarding 'NoneType'.")
return None
return diagonal_operator(domain=domain,diag=diag,bare=False)
@@ -8511,6 +8519,7 @@ class operator(object):
domain = self.target
diag = self.inverse_diag(bare=False,domain=domain,var=False,**kwargs)
if(diag is None):
+ about.warnings.cprint("WARNING: forwarding 'NoneType'.")
return None
return diagonal_operator(domain=domain,diag=diag,bare=False)
diff --git a/nifty_tools.py b/nifty_tools.py
index db7ce73f5d004f738466cd58b5183f013d55b77d..2803b5e9826f97cdb9dc4b19b6c4ca1e194a8872 100644
--- a/nifty_tools.py
+++ b/nifty_tools.py
@@ -753,7 +753,7 @@ class conjugate_gradient(object):
d = r+beta*d
delta = delta_*np.absolute(gamma)**0.5
- self.note.cflush("\niteration : %08u alpha = %3.1E beta = %3.1E delta = %3.1E"%(ii,alpha,beta,delta))
+ self.note.cflush("\niteration : %08u alpha = %3.1E beta = %3.1E delta = %3.1E"%(ii,np.real(alpha),np.real(beta),np.real(delta)))
if(ii==limii):
self.note.cprint("\n... quit.")
break
@@ -815,7 +815,7 @@ class conjugate_gradient(object):
d = s+beta*d ## conjugated gradient
delta = delta_*np.absolute(gamma)**0.5
- self.note.cflush("\niteration : %08u alpha = %3.1E beta = %3.1E delta = %3.1E"%(ii,alpha,beta,delta))
+ self.note.cflush("\niteration : %08u alpha = %3.1E beta = %3.1E delta = %3.1E"%(ii,np.real(alpha),np.real(beta),np.real(delta)))
if(ii==limii):
self.note.cprint("\n... quit.")
break