Commit 0389c909 authored by Theo Steininger's avatar Theo Steininger Committed by Pumpe, Daniel (dpumpe)

Revert "Merge branch 'units' into 'master'

This reverts commit baf3a7b8, reversing
changes made to 8ce837d0.
parent 8dea6eb7
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 = 'not'
distribution_strategy = 'fftw'
# Setting up the geometry
s_space = RGSpace([512, 512], dtype=np.float64)
......@@ -51,8 +50,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 = 'not'
distribution_strategy = 'fftw'
# 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,6 +38,3 @@ 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
# 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 = res.weight(power=-1)
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, bare=True)
power_operator = DiagonalOperator(domain, diagonal=f)
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