Commit 3711b36b authored by Marco Selig's avatar Marco Selig

several fixes; demo_wf added.

parent 74237b2f
......@@ -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
## 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
## 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
......@@ -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)
......
......@@ -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
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment