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

Merge branch 'master' of gitlab.mpcdf.mpg.de:ift/NIFTy

parents 66c0072b e0f1f873
......@@ -2,6 +2,11 @@
setup.cfg
.idea
.DS_Store
*.pyc
*.html
.document
.svn/
*.csv
# from https://github.com/github/gitignore/blob/master/Python.gitignore
......@@ -95,4 +100,4 @@ ENV/
.spyderproject
# Rope project settings
.ropeproject
\ No newline at end of file
.ropeproject
......@@ -47,7 +47,9 @@ test_mpi_fftw_hdf5:
- ci/install_pyfftw.sh
- ci/install_h5py.sh
- python setup.py build_ext --inplace
- nosetests -vv --with-coverage --cover-package=nifty --cover-branches
- mpiexec --allow-run-as-root -n 2 nosetests -x
- mpiexec --allow-run-as-root -n 4 nosetests -x
- nosetests -x --with-coverage --cover-package=nifty --cover-branches
- >
coverage report | grep TOTAL | awk '{ print "TOTAL: "$6; }'
......
import numpy as np
from nifty import (VL_BFGS, DiagonalOperator, FFTOperator, Field,
LinearOperator, PowerSpace, RelaxedNewton, RGSpace,
SteepestDescent, create_power_operator, exp, log, sqrt)
from nifty.library.critical_filter import CriticalPowerEnergy
from nifty.library.wiener_filter import WienerFilterEnergy
import plotly.graph_objs as go
import plotly.offline as pl
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
np.random.seed(42)
def plot_parameters(m, t, p, p_d):
x = log(t.domain[0].kindex)
m = fft.adjoint_times(m)
m = m.val.get_full_data().real
t = t.val.get_full_data().real
p = p.val.get_full_data().real
p_d = p_d.val.get_full_data().real
pl.plot([go.Heatmap(z=m)], filename='map.html', auto_open=False)
pl.plot([go.Scatter(x=x, y=t), go.Scatter(x=x, y=p),
go.Scatter(x=x, y=p_d)], filename="t.html", auto_open=False)
class AdjointFFTResponse(LinearOperator):
def __init__(self, FFT, R, default_spaces=None):
super(AdjointFFTResponse, self).__init__(default_spaces)
self._domain = FFT.target
self._target = R.target
self.R = R
self.FFT = FFT
def _times(self, x, spaces=None):
return self.R(self.FFT.adjoint_times(x))
def _adjoint_times(self, x, spaces=None):
return self.FFT(self.R.adjoint_times(x))
@property
def domain(self):
return self._domain
@property
def target(self):
return self._target
@property
def unitary(self):
return False
if __name__ == "__main__":
distribution_strategy = 'not'
# Set up position space
s_space = RGSpace([128, 128])
# s_space = HPSpace(32)
# Define harmonic transformation and associated harmonic space
fft = FFTOperator(s_space)
h_space = fft.target[0]
# Set up power space
p_space = PowerSpace(h_space, logarithmic=True,
distribution_strategy=distribution_strategy)
# Choose the prior correlation structure and defining correlation operator
p_spec = (lambda k: (.5 / (k + 1) ** 3))
S = create_power_operator(h_space, power_spectrum=p_spec,
distribution_strategy=distribution_strategy)
# Draw a sample sh from the prior distribution in harmonic space
sp = Field(p_space, val=p_spec,
distribution_strategy=distribution_strategy)
sh = sp.power_synthesize(real_signal=True)
# Choose the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01)
Instrument = DiagonalOperator(s_space, diagonal=1.)
# Instrument._diagonal.val[200:400, 200:400] = 0
# Instrument._diagonal.val[64:512-64, 64:512-64] = 0
# Add a harmonic transformation to the instrument
R = AdjointFFTResponse(fft, Instrument)
noise = 1.
N = DiagonalOperator(s_space, diagonal=noise, bare=True)
n = Field.from_random(domain=s_space,
random_type='normal',
std=sqrt(noise),
mean=0)
# Create mock data
d = R(sh) + n
# The information source
j = R.adjoint_times(N.inverse_times(d))
realized_power = log(sh.power_analyze(binbounds=p_space.binbounds))
data_power = log(fft(d).power_analyze(binbounds=p_space.binbounds))
d_data = d.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html', auto_open=False)
# Minimization strategy
def convergence_measure(a_energy, iteration): # returns current energy
x = a_energy.value
print(x, iteration)
minimizer1 = RelaxedNewton(convergence_tolerance=1e-4,
convergence_level=1,
iteration_limit=5,
callback=convergence_measure)
minimizer2 = VL_BFGS(convergence_tolerance=1e-10,
convergence_level=1,
iteration_limit=30,
callback=convergence_measure,
max_history_length=20)
minimizer3 = SteepestDescent(convergence_tolerance=1e-4,
iteration_limit=100,
callback=convergence_measure)
# Set starting position
flat_power = Field(p_space, val=1e-8)
m0 = flat_power.power_synthesize(real_signal=True)
t0 = Field(p_space, val=log(1./(1+p_space.kindex)**2))
for i in range(500):
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
# Initialize non-linear Wiener Filter energy
map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
# Solve the Wiener Filter analytically
D0 = map_energy.curvature
m0 = D0.inverse_times(j)
# Initialize power energy with updated parameters
power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0,
smoothness_prior=10., samples=3)
(power_energy, convergence) = minimizer2(power_energy)
# Set new power spectrum
t0.val = power_energy.position.val.real
# Plot current estimate
print(i)
if i % 5 == 0:
plot_parameters(m0, t0, log(sp), data_power)
# -*- coding: utf-8 -*-
from nifty import *
if __name__ == "__main__":
nifty_configuration['default_distribution_strategy'] = 'fftw'
# Setting up parameters |\label{code:wf_parameters}|
correlation_length = 1. # Typical distance over which the field is correlated
field_variance = 2. # Variance of field in position space
response_sigma = 0.02 # Smoothing length of response (in same unit as L)
signal_to_noise = 100 # The signal to noise ratio
np.random.seed(43) # Fixing the random seed
def power_spectrum(k): # Defining the power spectrum
a = 4 * correlation_length * field_variance**2
return a / (1 + k * correlation_length) ** 4
# Setting up the geometry |\label{code:wf_geometry}|
L = 2. # Total side-length of the domain
N_pixels = 128 # Grid resolution (pixels per axis)
#signal_space = RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
signal_space = HPSpace(16)
harmonic_space = FFTOperator.get_default_codomain(signal_space)
fft = FFTOperator(harmonic_space, target=signal_space, target_dtype=np.float)
power_space = PowerSpace(harmonic_space)
# Creating the mock signal |\label{code:wf_mock_signal}|
S = create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = Field(power_space, val=power_spectrum)
mock_signal = fft(mock_power.power_synthesize(real_signal=True))
# Setting up an exemplary response
mask = Field(signal_space, val=1.)
N10 = int(N_pixels/10)
#mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
R = ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}|
data_domain = R.target[0]
R_harmonic = ComposedOperator([fft, R], default_spaces=[0, 0])
# Setting up the noise covariance and drawing a random noise realization
N = DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise, bare=True)
noise = Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(exp(mock_signal)) + noise #|\label{code:wf_mock_data}|
# Wiener filter
m0 = Field(harmonic_space, val=0.j)
energy = library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, N, S)
minimizer1 = VL_BFGS(convergence_tolerance=1e-5,
iteration_limit=3000,
#callback=convergence_measure,
max_history_length=20)
minimizer2 = RelaxedNewton(convergence_tolerance=1e-5,
iteration_limit=10,
#callback=convergence_measure
)
minimizer3 = SteepestDescent(convergence_tolerance=1e-5, iteration_limit=1000)
# me1 = minimizer1(energy)
# me2 = minimizer2(energy)
# me3 = minimizer3(energy)
# m1 = fft(me1[0].position)
# m2 = fft(me2[0].position)
# m3 = fft(me3[0].position)
#
# # Probing the variance
# class Proby(DiagonalProberMixin, Prober): pass
# proby = Proby(signal_space, probe_count=100)
# proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z))))
#
# sm = SmoothingOperator(signal_space, sigma=0.02)
# variance = sm(proby.diagonal.weight(-1))
#Plotting #|\label{code:wf_plotting}|
#plotter = plotting.RG2DPlotter(color_map=plotting.colormaps.PlankCmap())
plotter = plotting.HealpixPlotter(color_map=plotting.colormaps.PlankCmap())
plotter.figure.xaxis = plotting.Axis(label='Pixel Index')
plotter.figure.yaxis = plotting.Axis(label='Pixel Index')
plotter.plot.zmax = 5; plotter.plot.zmin = -5
## plotter(variance, path = 'variance.html')
# #plotter.plot.zmin = exp(mock_signal.min());
# plotter(mock_signal.real, path='mock_signal.html')
# plotter(Field(signal_space, val=np.log(data.val.get_full_data().real).reshape(signal_space.shape)),
# path = 'log_of_data.html')
#
# plotter(m1.real, path='m_LBFGS.html')
# plotter(m2.real, path='m_Newton.html')
# plotter(m3.real, path='m_SteepestDescent.html')
#
# -*- coding: utf-8 -*-
import numpy as np
import nifty as ift
from nifty import plotting
from keepers import Repository
if __name__ == "__main__":
signal_to_noise = 1.5 # The signal to noise ratioa
# Setting up parameters |\label{code:wf_parameters}|
correlation_length_1 = 1. # Typical distance over which the field is correlated
field_variance_1 = 2. # Variance of field in position space
response_sigma_1 = 0.05 # Smoothing length of response (in same unit as L)
def power_spectrum_1(k): # note: field_variance**2 = a*k_0/4.
a = 4 * correlation_length_1 * field_variance_1**2
return a / (1 + k * correlation_length_1) ** 4.
# Setting up the geometry |\label{code:wf_geometry}|
L_1 = 2. # Total side-length of the domain
N_pixels_1 = 512 # Grid resolution (pixels per axis)
signal_space_1 = ift.RGSpace([N_pixels_1], distances=L_1/N_pixels_1)
harmonic_space_1 = ift.FFTOperator.get_default_codomain(signal_space_1)
fft_1 = ift.FFTOperator(harmonic_space_1, target=signal_space_1,
domain_dtype=np.complex, target_dtype=np.complex)
power_space_1 = ift.PowerSpace(harmonic_space_1, distribution_strategy='fftw')
mock_power_1 = ift.Field(power_space_1, val=power_spectrum_1,
distribution_strategy='not')
# Setting up parameters |\label{code:wf_parameters}|
correlation_length_2 = 1. # Typical distance over which the field is correlated
field_variance_2 = 2. # Variance of field in position space
response_sigma_2 = 0.01 # Smoothing length of response (in same unit as L)
def power_spectrum_2(k): # note: field_variance**2 = a*k_0/4.
a = 4 * correlation_length_2 * field_variance_2**2
return a / (1 + k * correlation_length_2) ** 2.5
# Setting up the geometry |\label{code:wf_geometry}|
L_2 = 2. # Total side-length of the domain
N_pixels_2 = 512 # Grid resolution (pixels per axis)
signal_space_2 = ift.RGSpace([N_pixels_2], distances=L_2/N_pixels_2)
harmonic_space_2 = ift.FFTOperator.get_default_codomain(signal_space_2)
fft_2 = ift.FFTOperator(harmonic_space_2, target=signal_space_2,
domain_dtype=np.complex, target_dtype=np.complex)
power_space_2 = ift.PowerSpace(harmonic_space_2, distribution_strategy='not')
mock_power_2 = ift.Field(power_space_2, val=power_spectrum_2,
distribution_strategy='not')
fft = ift.ComposedOperator((fft_1, fft_2))
mock_power = ift.Field(domain=(power_space_1, power_space_2),
val=np.outer(mock_power_1.val.get_full_data(),
mock_power_2.val.get_full_data()),
distribution_strategy='not')
diagonal = mock_power.power_synthesize(spaces=(0, 1), mean=1, std=0,
real_signal=False,
distribution_strategy='fftw')**2
S = ift.DiagonalOperator(domain=(harmonic_space_1, harmonic_space_2),
diagonal=diagonal)
np.random.seed(10)
mock_signal = fft(mock_power.power_synthesize(real_signal=True,
distribution_strategy='fftw'))
# Setting up a exemplary response
N1_10 = int(N_pixels_1/10)
mask_1 = ift.Field(signal_space_1, val=1., distribution_strategy='fftw')
mask_1.val[N1_10*7:N1_10*9] = 0.
N2_10 = int(N_pixels_2/10)
mask_2 = ift.Field(signal_space_2, val=1., distribution_strategy='not')
mask_2.val[N2_10*7:N2_10*9] = 0.
R = ift.ResponseOperator((signal_space_1, signal_space_2),
sigma=(response_sigma_1, response_sigma_2),
exposure=(mask_1, mask_2)) #|\label{code:wf_response}|
data_domain = R.target
R_harmonic = ift.ComposedOperator([fft, R], default_spaces=(0, 1, 0, 1))
# Setting up the noise covariance and drawing a random noise realization
N = ift.DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise,
bare=True,
distribution_strategy='fftw')
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise),
mean=0,
distribution_strategy='fftw')
data = R(mock_signal) + noise #|\label{code:wf_mock_data}|
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data))
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic)
wiener_curvature._InvertibleOperatorMixin__inverter.convergence_tolerance = 1e-3
m_k = wiener_curvature.inverse_times(j) #|\label{code:wf_wiener_filter}|
m = fft(m_k)
# Probing the variance
class Proby(ift.DiagonalProberMixin, ift.Prober): pass
proby = Proby((signal_space_1, signal_space_2), probe_count=100)
proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z))))
# sm = SmoothingOperator(signal_space, sigma=0.02)
# variance = sm(proby.diagonal.weight(-1))
variance = proby.diagonal.weight(-1)
repo = Repository('repo_100.h5')
repo.add(mock_signal, 'mock_signal')
repo.add(data, 'data')
repo.add(m, 'm')
repo.add(variance, 'variance')
repo.commit()
plot_space = ift.RGSpace((N_pixels_1, N_pixels_2))
plotter = plotting.RG2DPlotter(color_map=plotting.colormaps.PlankCmap())
plotter.figure.xaxis = ift.plotting.Axis(label='Pixel Index')
plotter.figure.yaxis = ift.plotting.Axis(label='Pixel Index')
plotter.plot.zmin = 0.
plotter.plot.zmax = 3.
sm = ift.SmoothingOperator(plot_space, sigma=0.03)
plotter(ift.log(ift.sqrt(sm(ift.Field(plot_space, val=variance.val.real)))), path='uncertainty.html')
plotter.plot.zmin = np.real(mock_signal.min());
plotter.plot.zmax = np.real(mock_signal.max());
plotter(ift.Field(plot_space, val=mock_signal.val.real), path='mock_signal.html')
plotter(ift.Field(plot_space, val=data.val.get_full_data().real), path = 'data.html')
plotter(ift.Field(plot_space, val=m.val.real), path = 'map.html')
# -*- coding: utf-8 -*-
import nifty as ift
from nifty import plotting
import numpy as np
from keepers import Repository
if __name__ == "__main__":
ift.nifty_configuration['default_distribution_strategy'] = 'fftw'
# Setting up parameters |\label{code:wf_parameters}|
correlation_length_scale = 1. # Typical distance over which the field is correlated
fluctuation_scale = 2. # Variance of field in position space
response_sigma = 0.05 # Smoothing length of response (in same unit as L)
signal_to_noise = 1.5 # The signal to noise ratio
np.random.seed(43) # Fixing the random seed
def power_spectrum(k): # Defining the power spectrum
a = 4 * correlation_length_scale * fluctuation_scale**2
return a / (1 + (k * correlation_length_scale)**2) ** 2
# Setting up the geometry |\label{code:wf_geometry}|
L = 2. # Total side-length of the domain
N_pixels = 512 # Grid resolution (pixels per axis)
signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
harmonic_space = ift.FFTOperator.get_default_codomain(signal_space)
fft = ift.FFTOperator(harmonic_space, target=signal_space, target_dtype=np.float)
power_space = ift.PowerSpace(harmonic_space)
# Creating the mock signal |\label{code:wf_mock_signal}|
S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = ift.Field(power_space, val=power_spectrum)
mock_signal = fft(mock_power.power_synthesize(real_signal=True))
# Setting up an exemplary response
mask = ift.Field(signal_space, val=1.)
N10 = int(N_pixels/10)
mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}|
data_domain = R.target[0]
R_harmonic = ift.ComposedOperator([fft, R], default_spaces=[0, 0])
# Setting up the noise covariance and drawing a random noise realization
N = ift.DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise, bare=True)
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(mock_signal) + noise #|\label{code:wf_mock_data}|
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data))
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic)
m_k = wiener_curvature.inverse_times(j) #|\label{code:wf_wiener_filter}|
m = fft(m_k)
# Probing the uncertainty |\label{code:wf_uncertainty_probing}|
class Proby(ift.DiagonalProberMixin, ift.Prober): pass
proby = Proby(signal_space, probe_count=800)
proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z)))) #|\label{code:wf_variance_fft_wrap}|
sm = ift.SmoothingOperator(signal_space, sigma=0.03)
variance = ift.sqrt(sm(proby.diagonal.weight(-1))) #|\label{code:wf_variance_weighting}|
repo = Repository('repo_800.h5')
repo.add(mock_signal, 'mock_signal')
repo.add(data, 'data')
repo.add(m, 'm')
repo.add(variance, 'variance')
repo.commit()
# Plotting #|\label{code:wf_plotting}|
plotter = plotting.RG2DPlotter(color_map=plotting.colormaps.PlankCmap())
plotter.figure.xaxis = ift.plotting.Axis(label='Pixel Index')
plotter.figure.yaxis = ift.plotting.Axis(label='Pixel Index')
plotter.plot.zmax = variance.max(); plotter.plot.zmin = 0
plotter(variance, path = 'uncertainty.html')
plotter.plot.zmax = mock_signal.max(); plotter.plot.zmin = mock_signal.min()
plotter(mock_signal, path='mock_signal.html')
plotter(ift.Field(signal_space, val=data.val), path='data.html')
plotter(m, path='map.html')
from nifty import *
#import plotly.offline as pl
#import plotly.graph_objs as go
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
if __name__ == "__main__":
distribution_strategy = 'not'
#Setting up physical constants
#total length of Interval or Volume the field lives on, e.g. in meters
L = 2.
#typical distance over which the field is correlated (in same unit as L)
correlation_length = 0.1
#variance of field in position space sqrt(<|s_x|^2>) (in unit of s)
field_variance = 2.
#smoothing length that response (in same unit as L)
response_sigma = 0.1
#defining resolution (pixels per dimension)
N_pixels = 512
#Setting up derived constants
k_0 = 1./correlation_length
#note that field_variance**2 = a*k_0/4. for this analytic form of power
#spectrum
a = field_variance**2/k_0*4.
pow_spec = (lambda k: a / (1 + k/k_0) ** 4)
pixel_width = L/N_pixels
# Setting up the geometry
s_space = RGSpace([N_pixels, N_pixels], distances = pixel_width)
fft = FFTOperator(s_space)
h_space = fft.target[0]
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
# Creating the mock data
S = create_power_operator(h_space, power_spectrum=pow_spec,
distribution_strategy=distribution_strategy)
sp = Field(p_space, val=pow_spec,
distribution_strategy=distribution_strategy)
sh = sp.power_synthesize(real_signal=True)
ss = fft.inverse_times(sh)
R = SmoothingOperator(s_space, sigma=response_sigma)
signal_to_noise = 1
N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
n = Field.from_random(domain=s_space,
random_type='normal',
std=ss.std()/np.sqrt(signal_to_noise),
mean=0)
d = R(ss) + n
# Wiener filter
j = R.adjoint_times(N.inverse_times(d))
D = PropagatorOperator(S=S, N=N, R=R)
m = D(j)
d_data = d.val.get_full_data().real
m_data = m.val.get_full_data().real
ss_data = ss.val.get_full_data().real
</