Commit baf3a7b8 authored by Pumpe, Daniel (dpumpe)'s avatar Pumpe, Daniel (dpumpe)

Merge branch 'units' into 'master'

Wiener filter example with units and a little plot routine.

See merge request !59
parents 8ce837d0 9612be8f
Pipeline #11299 passed with stages
in 30 minutes and 49 seconds
from nifty import *
#import plotly.offline as pl
#import plotly.graph_objs as go
import plotly.offline as pl
import plotly.graph_objs as go
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
np.random.seed(42)
if __name__ == "__main__":
distribution_strategy = 'fftw'
distribution_strategy = 'not'
# Setting up the geometry
s_space = RGSpace([512, 512], dtype=np.float64)
......@@ -50,8 +51,8 @@ if __name__ == "__main__":
d_data = d.val.get_full_data().real
m_data = m.val.get_full_data().real
ss_data = ss.val.get_full_data().real
# if rank == 0:
# pl.plot([go.Heatmap(z=d_data)], filename='data.html')
# pl.plot([go.Heatmap(z=m_data)], filename='map.html')
# pl.plot([go.Heatmap(z=ss_data)], filename='map_orig.html')
if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html')
pl.plot([go.Heatmap(z=m_data)], filename='map.html')
pl.plot([go.Heatmap(z=ss_data)], filename='map_orig.html')
#
......@@ -53,7 +53,7 @@ class WienerFilterEnergy(Energy):
if __name__ == "__main__":
distribution_strategy = 'fftw'
distribution_strategy = 'not'
# Set up spaces and fft transformation
s_space = RGSpace([512, 512], dtype=np.float)
......
from nifty import *
from mpi4py import MPI
import plotly.offline as py
import plotly.graph_objs as go
comm = MPI.COMM_WORLD
rank = comm.rank
def plot_maps(x, name):
trace = [None]*len(x)
keys = x.keys()
field = x[keys[0]]
domain = field.domain[0]
shape = len(domain.shape)
max_n = domain.shape[0]*domain.distances[0]
step = domain.distances[0]
x_axis = np.arange(0, max_n, step)
if shape == 1:
for ii in xrange(len(x)):
trace[ii] = go.Scatter(x= x_axis, y=x[keys[ii]].val.get_full_data(), name=keys[ii])
fig = go.Figure(data=trace)
py.plot(fig, filename=name)
elif shape == 2:
for ii in xrange(len(x)):
py.plot([go.Heatmap(z=x[keys[ii]].val.get_full_data())], filename=keys[ii])
else:
raise TypeError("Only 1D and 2D field plots are supported")
def plot_power(x, name):
layout = go.Layout(
xaxis=dict(
type='log',
autorange=True
),
yaxis=dict(
type='log',
autorange=True
)
)
trace = [None]*len(x)
keys = x.keys()
field = x[keys[0]]
domain = field.domain[0]
x_axis = domain.kindex
for ii in xrange(len(x)):
trace[ii] = go.Scatter(x= x_axis, y=x[keys[ii]].val.get_full_data(), name=keys[ii])
fig = go.Figure(data=trace, layout=layout)
py.plot(fig, filename=name)
np.random.seed(42)
if __name__ == "__main__":
distribution_strategy = 'not'
# setting spaces
npix = np.array([500]) # number of pixels
total_volume = 1. # total length
# setting signal parameters
lambda_s = .05 # signal correlation length
sigma_s = 10. # signal variance
#setting response operator parameters
length_convolution = .025
exposure = 2000.
# calculating parameters
k_0 = 4. / (2 * np.pi * lambda_s)
a_s = sigma_s ** 2. * lambda_s * total_volume
# creation of spaces
x1 = RGSpace(npix, dtype=np.float64, distances=total_volume / npix,
zerocenter=False)
k1 = RGRGTransformation.get_codomain(x1)
p1 = PowerSpace(harmonic_domain=k1, log=False, dtype=np.float64)
# creating Power Operator with given spectrum
spec = (lambda k: a_s / (1 + (k / k_0) ** 2) ** 2)
p_field = Field(p1, val=spec)
S_op = create_power_operator(k1, spec)
# creating FFT-Operator and Response-Operator with Gaussian convolution
Fft_op = FFTOperator(domain=x1, target=k1,
domain_dtype=np.float64,
target_dtype=np.complex128)
R_op = ResponseOperator(x1, sigma=length_convolution,
exposure=exposure)
# drawing a random field
sk = p_field.power_synthesize(real_signal=True, mean=0.)
s = Fft_op.inverse_times(sk)
signal_to_noise = 1
N_op = DiagonalOperator(R_op.target, diagonal=s.var()/signal_to_noise, bare=True)
n = Field.from_random(domain=R_op.target,
random_type='normal',
std=s.std()/np.sqrt(signal_to_noise),
mean=0)
d = R_op(s) + n
# Wiener filter
j = R_op.adjoint_times(N_op.inverse_times(d))
D = PropagatorOperator(S=S_op, N=N_op, R=R_op)
m = D(j)
z={}
z["signal"] = s
z["reconstructed_map"] = m
z["data"] = d
z["lambda"] = R_op(s)
plot_maps(z, "Wiener_filter.html")
......@@ -38,3 +38,6 @@ from projection_operator import ProjectionOperator
from propagator_operator import PropagatorOperator
from composed_operator import ComposedOperator
from response_operator import ResponseOperator
from response_operator import ResponseOperator
\ No newline at end of file
from nifty import Field
from nifty.field_types.field_array import FieldArray
from nifty.operators.linear_operator import LinearOperator
from nifty.operators.smoothing_operator import SmoothingOperator
import numpy as np
class ResponseOperator(LinearOperator):
def __init__(self, domain,
sigma=1., exposure=1., implemented=True,
unitary=False):
self._domain = self._parse_domain(domain)
self._target = self._parse_domain(FieldArray(self._domain[0].shape,
dtype=np.float64))
self._sigma = sigma
self._implemented = implemented
self._unitary = unitary
self._kernel = SmoothingOperator(self._domain,
sigma=self._sigma)
self._exposure = exposure
@property
def domain(self):
return self._domain
@property
def target(self):
return self._target
@property
def implemented(self):
return self._implemented
@property
def unitary(self):
return self._unitary
def _times(self, x, spaces):
res = self._kernel.times(x)
res = self._exposure * res
res = res.weight(power=1)
# removing geometric information
return Field(self._target, val=res.val)
def _adjoint_times(self, x, spaces):
# setting correct spaces
res = x*self._exposure
res = Field(self.domain, val=res.val)
res = self._kernel.adjoint_times(res)
return res
......@@ -24,6 +24,6 @@ def create_power_operator(domain, power_spectrum, dtype=None,
f = fp.power_synthesize(mean=1, std=0, real_signal=False)
power_operator = DiagonalOperator(domain, diagonal=f)
power_operator = DiagonalOperator(domain, diagonal=f, bare=True)
return power_operator
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