Commit e7d60d48 authored by Theo Steininger's avatar Theo Steininger

Merge branch 'master' into mpitests

parents 226bf0d1 5d5a2701
Pipeline #14854 failed with stage
in 7 minutes and 58 seconds
......@@ -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
from nifty import *
from nifty.library.wiener_filter import WienerFilterEnergy
from nifty.library.critical_filter import CriticalPowerEnergy
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)
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')
pl.plot([go.Scatter(x=x,y=t), go.Scatter(x=x ,y=p), go.Scatter(x=x, y=p_d)], filename="t.html")
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]
# Setting up power space
p_space = PowerSpace(h_space, logarithmic=True,
distribution_strategy=distribution_strategy)
# Choosing 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)
# Drawing 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)
# Choosing 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
#Adding 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)
# Creating the mock data
d = R(sh) + n
# The information source
j = R.adjoint_times(N.inverse_times(d))
realized_power = log(sh.power_analyze(logarithmic=p_space.config["logarithmic"],
nbin=p_space.config["nbin"]))
data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"],
nbin=p_space.config["nbin"]))
d_data = d.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html')
# minimization strategy
def convergence_measure(a_energy, iteration): # returns current energy
x = a_energy.value
print (x, iteration)
minimizer1 = RelaxedNewton(convergence_tolerance=1e-2,
convergence_level=2,
iteration_limit=3,
callback=convergence_measure)
minimizer2 = VL_BFGS(convergence_tolerance=0,
iteration_limit=7,
callback=convergence_measure,
max_history_length=3)
# Setting 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)
# Initializing the nonlinear Wiener Filter energy
map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
# Solving the Wiener Filter analytically
D0 = map_energy.curvature
m0 = D0.inverse_times(j)
# Initializing the power energy with updated parameters
power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, smoothness_prior=10., samples=3)
(power_energy, convergence) = minimizer1(power_energy)
# Setting new power spectrum
t0.val = power_energy.position.val.real
# Plotting current estimate
print i
if i%50 == 0:
plot_parameters(m0,t0,log(sp), data_power)
from nifty import *
import plotly.offline as pl
import plotly.graph_objs as go
from nifty.library.wiener_filter import *
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
np.random.seed(42)
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]
# Setting up power space
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
# Choosing the prior correlation structure and defining correlation operator
p_spec = (lambda k: (42 / (k + 1) ** 3))
S = create_power_operator(h_space, power_spectrum=p_spec,
distribution_strategy=distribution_strategy)
# Drawing 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)
ss = fft.adjoint_times(sh)
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.05)
Instrument = DiagonalOperator(s_space, diagonal=1.)
# Instrument._diagonal.val[200:400, 200:400] = 0
#Adding a harmonic transformation to the instrument
R = AdjointFFTResponse(fft, Instrument)
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)
# Creating the mock data
d = R(sh) + n
j = R.adjoint_times(N.inverse_times(d))
# Choosing the minimization strategy
def convergence_measure(energy, iteration): # returns current energy
x = energy.value
print (x, iteration)
# minimizer = SteepestDescent(convergence_tolerance=0,
# iteration_limit=50,
# callback=convergence_measure)
minimizer = RelaxedNewton(convergence_tolerance=0,
iteration_limit=1,
callback=convergence_measure)
#
# minimizer = VL_BFGS(convergence_tolerance=0,
# iteration_limit=50,
# callback=convergence_measure,
# max_history_length=3)
#
inverter = ConjugateGradient(convergence_level=3,
convergence_tolerance=1e-5,
preconditioner=None)
# Setting starting position
m0 = Field(h_space, val=.0)
# Initializing the Wiener Filter energy
energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S, inverter=inverter)
D0 = energy.curvature
# Solving the problem analytically
m0 = D0.inverse_times(j)
sample_variance = Field(sh.domain,val=0. + 0j)
sample_mean = Field(sh.domain,val=0. + 0j)
# sampling the uncertainty map
n_samples = 1
for i in range(n_samples):
sample = sugar.generate_posterior_sample(m0,D0)
sample_variance += sample**2
sample_mean += sample
variance = sample_variance/n_samples - (sample_mean/n_samples)
import numpy as np
from nifty import RGSpace, PowerSpace, Field, FFTOperator, ComposedOperator,\
SmoothingOperator, DiagonalOperator, create_power_operator
from nifty.library import WienerFilterCurvature
from nifty import *
#import plotly.offline as pl
#import plotly.graph_objs as go
......@@ -10,36 +14,37 @@ 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
distribution_strategy = 'fftw'
# 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)
# 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)
# 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)
# smoothing length of response (in same unit as L)
response_sigma = 0.1
#defining resolution (pixels per dimension)
# defining resolution (pixels per dimension)
N_pixels = 512
#Setting up derived constants
# 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
# 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
pixel_length = L/N_pixels
# Setting up the geometry
s_space = RGSpace([N_pixels, N_pixels], distances = pixel_width)
fft = FFTOperator(s_space)
s_space = RGSpace([N_pixels, N_pixels], distances=pixel_length)
fft = FFTOperator(s_space, domain_dtype=np.float, target_dtype=np.complex)
h_space = fft.target[0]
inverse_fft = FFTOperator(h_space, target=s_space,
domain_dtype=np.complex, target_dtype=np.float)
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
# Creating the mock data
S = create_power_operator(h_space, power_spectrum=pow_spec,
......@@ -51,6 +56,7 @@ if __name__ == "__main__":
ss = fft.inverse_times(sh)
R = SmoothingOperator(s_space, sigma=response_sigma)
R_harmonic = ComposedOperator([inverse_fft, R], default_spaces=[0, 0])
signal_to_noise = 1
N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
......@@ -62,15 +68,10 @@ if __name__ == "__main__":
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
# 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')
j = R_harmonic.adjoint_times(N.inverse_times(d))
wiener_curvature = WienerFilterCurvature(S=S, N=N, R=R_harmonic)
m = wiener_curvature.inverse_times(j)
m_s = inverse_fft(m)
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
np.random.seed(42)
class WienerFilterEnergy(Energy):
def __init__(self, position, D, j):
# in principle not necessary, but useful in order to make the signature
# explicit
super(WienerFilterEnergy, self).__init__(position)
self.D = D
self.j = j
def at(self, position):
return self.__class__(position, D=self.D, j=self.j)
@property
def value(self):
D_inv_x = self.D_inverse_x()
H = 0.5 * D_inv_x.vdot(self.position) - self.j.vdot(self.position)
return H.real
@property
def gradient(self):
D_inv_x = self.D_inverse_x()
g = D_inv_x - self.j
return_g = g.copy_empty(dtype=np.float)
return_g.val = g.val.real
return return_g
@property
def curvature(self):
class Dummy(object):
def __init__(self, x):
self.x = x
def inverse_times(self, *args, **kwargs):
return self.x.times(*args, **kwargs)
my_dummy = Dummy(self.D)
return my_dummy
@memo
def D_inverse_x(self):
return D.inverse_times(self.position)
if __name__ == "__main__":
distribution_strategy = 'not'
# Set up spaces and fft transformation
s_space = RGSpace([512, 512])
fft = FFTOperator(s_space)
h_space = fft.target[0]
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
# create the field instances and power operator
pow_spec = (lambda k: (42 / (k + 1) ** 3))
S = create_power_operator(h_space, power_spectrum=pow_spec,
distribution_strategy=distribution_strategy)
sp = Field(p_space, val=lambda z: pow_spec(z)**(1./2),
distribution_strategy=distribution_strategy)
sh = sp.power_synthesize(real_signal=True)
ss = fft.inverse_times(sh)
# model the measurement process
R = SmoothingOperator(s_space, sigma=0.01)
# R = DiagonalOperator(s_space, diagonal=1.)
# R._diagonal.val[200:400, 200:400] = 0
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)
# create mock data
d = R(ss) + n
# set up reconstruction objects
j = R.adjoint_times(N.inverse_times(d))
D = PropagatorOperator(S=S, N=N, R=R)
def distance_measure(energy, iteration):
x = energy.position
print (iteration, ((x-ss).norm()/ss.norm()).real)
# minimizer = SteepestDescent(convergence_tolerance=0,
# iteration_limit=50,
# callback=distance_measure)
minimizer = RelaxedNewton(convergence_tolerance=0,
iteration_limit=2,
callback=distance_measure)
# minimizer = VL_BFGS(convergence_tolerance=0,
# iteration_limit=50,
# callback=distance_measure,
# max_history_length=3)
m0 = Field(s_space, val=1.)
energy = WienerFilterEnergy(position=m0, D=D, j=j)
(energy, convergence) = minimizer(energy)
m = energy.position
d_data = d.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html')
ss_data = ss.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=ss_data)], filename='ss.html')
sh_data = sh.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=sh_data)], filename='sh.html')
j_data = j.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=j_data)], filename='j.html')
jabs_data = np.abs(j.val.get_full_data())
jphase_data = np.angle(j.val.get_full_data())
if rank == 0:
pl.plot([go.Heatmap(z=jabs_data)], filename='j_abs.html')
pl.plot([go.Heatmap(z=jphase_data)], filename='j_phase.html')
m_data = m.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=m_data)], filename='map.html')
# grad_data = grad.val.get_full_data().real
# if rank == 0:
# pl.plot([go.Heatmap(z=grad_data)], filename='grad.html')
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