Commit ae697a00 authored by Jakob Knollmueller's avatar Jakob Knollmueller

remove all old demos, GN prototype in

parent 67e96903
import nifty5 as ift
from nifty5.library.nonlinearities import Linear, Tanh, Exponential
import numpy as np
np.random.seed(42)
def adjust_zero_mode(m0, t0):
mtmp = m0.to_global_data().copy()
zero_position = len(m0.shape)*(0,)
zero_mode = mtmp[zero_position]
mtmp[zero_position] = zero_mode / abs(zero_mode)
ttmp = t0.to_global_data().copy()
ttmp[0] += 2 * np.log(abs(zero_mode))
return (ift.Field.from_global_data(m0.domain, mtmp),
ift.Field.from_global_data(t0.domain, ttmp))
if __name__ == "__main__":
noise_level = 1.
p_spec = (lambda k: (.5 / (k + 1) ** 3))
nonlinearity = Linear()
# Set up position space
s_space = ift.RGSpace((128, 128))
h_space = s_space.get_default_codomain()
# Define harmonic transformation and associated harmonic space
HT = ift.HarmonicTransformOperator(h_space, target=s_space)
# Setting up power space
p_space = ift.PowerSpace(h_space,
binbounds=ift.PowerSpace.useful_binbounds(
h_space, logarithmic=True))
# Choosing the prior correlation structure and defining
# correlation operator
p = ift.PS_field(p_space, p_spec)
log_p = ift.log(p)
S = ift.create_power_operator(h_space, power_spectrum=lambda k: 1e-5)
# Drawing a sample sh from the prior distribution in harmonic space
sh = S.draw_sample()
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01)
mask = np.ones(s_space.shape)
# mask[6000:8000] = 0.
mask[30:70, 30:70] = 0.
mask = ift.Field.from_global_data(s_space, mask)
MaskOperator = ift.DiagonalOperator(mask)
R = ift.GeometryRemover(s_space)
R = R*MaskOperator
# R = R*HT
# R = R * ift.create_harmonic_smoothing_operator((harmonic_space,), 0,
# response_sigma)
MeasurementOperator = R
d_space = MeasurementOperator.target
Distributor = ift.PowerDistributor(target=h_space, power_space=p_space)
power = Distributor(ift.exp(0.5*log_p))
# Creating the mock data
true_sky = nonlinearity(HT(power*sh))
noiseless_data = MeasurementOperator(true_sky)
noise_amplitude = noiseless_data.val.std()*noise_level
N = ift.ScalingOperator(noise_amplitude**2, d_space)
n = N.draw_sample()
# Creating the mock data
d = noiseless_data + n
m0 = ift.full(h_space, 1e-7)
t0 = ift.full(p_space, -4.)
power0 = Distributor.times(ift.exp(0.5 * t0))
plotdict = {"colormap": "Planck-like"}
zmin = true_sky.min()
zmax = true_sky.max()
ift.plot(true_sky, title="True sky", name="true_sky.png", **plotdict)
ift.plot(MeasurementOperator.adjoint_times(d), title="Data",
name="data.png", **plotdict)
IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
tol_abs_gradnorm=1e-3)
LS = ift.LineSearchStrongWolfe(c2=0.02)
minimizer = ift.RelaxedNewton(IC1, line_searcher=LS)
IC = ift.GradientNormController(iteration_limit=500,
tol_abs_gradnorm=1e-3)
for i in range(20):
power0 = Distributor(ift.exp(0.5*t0))
map0_energy = ift.library.NonlinearWienerFilterEnergy(
m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S, IC)
# Minimization with chosen minimizer
map0_energy, convergence = minimizer(map0_energy)
m0 = map0_energy.position
# Updating parameters for correlation structure reconstruction
D0 = map0_energy.curvature
# Initializing the power energy with updated parameters
power0_energy = ift.library.NonlinearPowerEnergy(
position=t0, d=d, N=N, xi=m0, D=D0, ht=HT,
Instrument=MeasurementOperator, nonlinearity=nonlinearity,
Distributor=Distributor, sigma=1., samples=2,
iteration_controller=IC)
power0_energy = minimizer(power0_energy)[0]
# Setting new power spectrum
t0 = power0_energy.position
# break degeneracy between amplitude and excitation by setting
# excitation monopole to 1
m0, t0 = adjust_zero_mode(m0, t0)
ift.plot(nonlinearity(HT(power0*m0)), title="Reconstructed sky",
name="reconstructed_sky.png", zmin=zmin, zmax=zmax, **plotdict)
ymin = np.min(p.to_global_data())
ift.plot([ift.exp(t0), p], title="Power spectra",
name="ps.png", ymin=ymin, **plotdict)
import nifty5 as ift
import numpy as np
from global_newton.models_other.apply_data import ApplyData
from global_newton.models_energy.hamiltonian import Hamiltonian
from nifty5.library.unit_log_gauss import UnitLogGauss
if __name__ == '__main__':
# s_space = ift.RGSpace([1024])
s_space = ift.RGSpace([128,128])
# s_space = ift.HPSpace(64)
h_space = s_space.get_default_codomain()
total_domain = ift.MultiDomain.make({'xi': h_space})
HT = ift.HarmonicTransformOperator(h_space, s_space)
def sqrtpspec(k):
return 16. / (20.+k**2)
GR = ift.GeometryRemover(s_space)
d_space = GR.target
B = ift.FFTSmoothingOperator(s_space,0.1)
mask = np.ones(s_space.shape)
mask[64:89,76:100] = 0.
mask = ift.Field(s_space,val=mask)
Mask = ift.DiagonalOperator(mask)
R = GR * Mask * B
noise = 1.
N = ift.ScalingOperator(noise, d_space)
p_space = ift.PowerSpace(h_space)
pd = ift.PowerDistributor(h_space, p_space)
position = ift.from_random('normal', total_domain)
xi = ift.Variable(position)['xi']
a = ift.Constant(position, ift.PS_field(p_space, sqrtpspec))
A = pd(a)
s_h = A * xi
s = HT(s_h)
Rs = R(s)
MOCK_POSITION = ift.from_random('normal',total_domain)
data = Rs.at(MOCK_POSITION).value + N.draw_sample()
NWR = ApplyData(data, ift.Field(d_space,val=noise), Rs)
INITIAL_POSITION = ift.from_random('normal',total_domain)
likelihood = UnitLogGauss(INITIAL_POSITION, NWR)
IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3)
inverter = ift.ConjugateGradient(controller=IC)
IC2 = ift.GradientNormController(name='Newton', iteration_limit=15)
minimizer = ift.RelaxedNewton(IC2)
H = Hamiltonian(likelihood, inverter)
H, convergence = minimizer(H)
result = s.at(H.position).value
import nifty5 as ift
from nifty5.library.nonlinearities import Exponential
import numpy as np
np.random.seed(42)
def adjust_zero_mode(m0, t0):
mtmp = m0.to_global_data().copy()
zero_position = len(m0.shape)*(0,)
zero_mode = mtmp[zero_position]
mtmp[zero_position] = zero_mode / abs(zero_mode)
ttmp = t0.to_global_data().copy()
ttmp[0] += 2 * np.log(abs(zero_mode))
return (ift.Field.from_global_data(m0.domain, mtmp),
ift.Field.from_global_data(t0.domain, ttmp))
if __name__ == "__main__":
noise_level = 1.
p_spec = (lambda k: (1. / (k + 1) ** 2))
# nonlinearity = Linear()
nonlinearity = Exponential()
# Set up position space
# s_space = ift.RGSpace([1024])
s_space = ift.HPSpace(32)
h_space = s_space.get_default_codomain()
# Define harmonic transformation and associated harmonic space
HT = ift.HarmonicTransformOperator(h_space, target=s_space)
# Setting up power space
p_space = ift.PowerSpace(h_space,
binbounds=ift.PowerSpace.useful_binbounds(
h_space, logarithmic=True))
# Choosing the prior correlation structure and defining
# correlation operator
p = ift.PS_field(p_space, p_spec)
log_p = ift.log(p)
S = ift.create_power_operator(h_space, lambda k: 1.)
# Drawing a sample sh from the prior distribution in harmonic space
sh = S.draw_sample()
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01)
mask = np.ones(s_space.shape)
mask[6000:8000] = 0.
mask = ift.Field.from_global_data(s_space, mask)
MaskOperator = ift.DiagonalOperator(mask)
R = ift.GeometryRemover(s_space)
R = R*MaskOperator
# R = R*HT
# R = R * ift.create_harmonic_smoothing_operator((harmonic_space,), 0,
# response_sigma)
MeasurementOperator = R
d_space = MeasurementOperator.target
Distributor = ift.PowerDistributor(target=h_space, power_space=p_space)
power = Distributor(ift.exp(0.5*log_p))
# Creating the mock data
true_sky = nonlinearity(HT(power*sh))
noiseless_data = MeasurementOperator(true_sky)
noise_amplitude = noiseless_data.val.std()*noise_level
N = ift.ScalingOperator(noise_amplitude**2, d_space)
n = N.draw_sample()
# Creating the mock data
d = noiseless_data + n
m0 = ift.full(h_space, 1e-7)
t0 = ift.full(p_space, -4.)
power0 = Distributor.times(ift.exp(0.5 * t0))
IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
tol_abs_gradnorm=1e-3)
LS = ift.LineSearchStrongWolfe(c2=0.02)
minimizer = ift.RelaxedNewton(IC1, line_searcher=LS)
IC = ift.GradientNormController(iteration_limit=500,
tol_abs_gradnorm=1e-3)
for i in range(20):
power0 = Distributor(ift.exp(0.5*t0))
map0_energy = ift.library.NonlinearWienerFilterEnergy(
m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S, IC)
# Minimization with chosen minimizer
map0_energy, convergence = minimizer(map0_energy)
m0 = map0_energy.position
# Updating parameters for correlation structure reconstruction
D0 = map0_energy.curvature
# Initializing the power energy with updated parameters
power0_energy = ift.library.NonlinearPowerEnergy(
position=t0, d=d, N=N, xi=m0, D=D0, ht=HT,
Instrument=MeasurementOperator, nonlinearity=nonlinearity,
Distributor=Distributor, sigma=1., samples=2, iteration_controller=IC)
power0_energy = minimizer(power0_energy)[0]
# Setting new power spectrum
t0 = power0_energy.position
# break degeneracy between amplitude and excitation by setting
# excitation monopole to 1
m0, t0 = adjust_zero_mode(m0, t0)
plotdict = {"colormap": "Planck-like"}
ift.plot(true_sky, name="true_sky.png", **plotdict)
ift.plot(nonlinearity(HT(power0*m0)),
name="reconstructed_sky.png", **plotdict)
ift.plot(MeasurementOperator.adjoint_times(d), name="data.png", **plotdict)
ift.plot([ift.exp(t0), p], name="ps.png")
import nifty5 as ift
from nifty5.library.nonlinearities import Linear, Exponential, Tanh
import numpy as np
np.random.seed(20)
if __name__ == "__main__":
noise_level = 0.3
correlation_length = 0.1
p_spec = lambda k: (1. / (k*correlation_length + 1) ** 4)
nonlinearity = Tanh()
# nonlinearity = Linear()
# nonlinearity = Exponential()
# Set up position space
s_space = ift.RGSpace(1024)
h_space = s_space.get_default_codomain()
# Define harmonic transformation and associated harmonic space
HT = ift.HarmonicTransformOperator(h_space, target=s_space)
S = ift.ScalingOperator(1., h_space)
# Drawing a sample sh from the prior distribution in harmonic space
sh = S.draw_sample()
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01)
mask = np.ones(s_space.shape)
mask[600:800] = 0.
mask = ift.Field.from_global_data(s_space, mask)
R = ift.GeometryRemover(s_space) * ift.DiagonalOperator(mask)
d_space = R.target
p_op = ift.create_power_operator(h_space, p_spec)
power = ift.sqrt(p_op(ift.full(h_space, 1.)))
# Creating the mock data
true_sky = nonlinearity(HT(power*sh))
noiseless_data = R(true_sky)
noise_amplitude = noiseless_data.val.std()*noise_level
N = ift.ScalingOperator(noise_amplitude**2, d_space)
n = N.draw_sample()
# Creating the mock data
d = noiseless_data + n
IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
tol_abs_gradnorm=1e-4)
LS = ift.LineSearchStrongWolfe(c2=0.02)
minimizer = ift.RelaxedNewton(IC1, line_searcher=LS)
IC = ift.GradientNormController(iteration_limit=2000,
tol_abs_gradnorm=1e-3)
# initial guess
m = ift.full(h_space, 1e-7)
map_energy = ift.library.NonlinearWienerFilterEnergy(
m, d, R, nonlinearity, HT, power, N, S, IC)
# Minimization with chosen minimizer
map_energy, convergence = minimizer(map_energy)
m = map_energy.position
recsky = nonlinearity(HT(power*m))
data = R.adjoint_times(d)
lo = np.min([true_sky.min(), recsky.min(), data.min()])
hi = np.max([true_sky.max(), recsky.max(), data.max()])
plotdict = {"colormap": "Planck-like", "ymin": lo, "ymax": hi}
ift.plot(true_sky, name="true_sky.png", **plotdict)
ift.plot(recsky, name="reconstructed_sky.png", **plotdict)
ift.plot(data, name="data.png", **plotdict)
import numpy as np
import nifty5 as ift
if __name__ == "__main__":
signal_to_noise = 0.5 # The signal to noise ratio
# Setting up parameters
L_1 = 2. # Total side-length of the domain
N_pixels_1 = 512 # Grid resolution (pixels per axis)
L_2 = 2. # Total side-length of the domain
N_pixels_2 = 512 # Grid resolution (pixels per axis)
correlation_length_1 = 1.
field_variance_1 = 2. # Variance of field in position space
response_sigma_1 = 0.05 # Smoothing length of response
correlation_length_2 = 1.
field_variance_2 = 2. # Variance of field in position space
response_sigma_2 = 0.01 # Smoothing length of response
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.
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
signal_space_1 = ift.RGSpace([N_pixels_1], distances=L_1/N_pixels_1)
harmonic_space_1 = signal_space_1.get_default_codomain()
signal_space_2 = ift.RGSpace([N_pixels_2], distances=L_2/N_pixels_2)
harmonic_space_2 = signal_space_2.get_default_codomain()
signal_domain = ift.DomainTuple.make((signal_space_1, signal_space_2))
harmonic_domain = ift.DomainTuple.make((harmonic_space_1,
harmonic_space_2))
ht_1 = ift.HarmonicTransformOperator(harmonic_domain, space=0)
ht_2 = ift.HarmonicTransformOperator(ht_1.target, space=1)
ht = ht_2*ht_1
S = (ift.create_power_operator(harmonic_domain, power_spectrum_1, 0) *
ift.create_power_operator(harmonic_domain, power_spectrum_2, 1))
np.random.seed(10)
mock_signal = S.draw_sample()
# Setting up a exemplary response
N1_10 = int(N_pixels_1/10)
mask_1 = np.ones(signal_space_1.shape)
mask_1[N1_10*7:N1_10*9] = 0.
mask_1 = ift.Field.from_global_data(signal_space_1, mask_1)
N2_10 = int(N_pixels_2/10)
mask_2 = np.ones(signal_space_2.shape)
mask_2[N2_10*7:N2_10*9] = 0.
mask_2 = ift.Field.from_global_data(signal_space_2, mask_2)
R = ift.GeometryRemover(signal_domain)
R = R*ift.DiagonalOperator(mask_1, signal_domain, spaces=0)
R = R*ift.DiagonalOperator(mask_2, signal_domain, spaces=1)
R = R*ht
R = R * ift.create_harmonic_smoothing_operator(harmonic_domain, 0,
response_sigma_1)
R = R * ift.create_harmonic_smoothing_operator(harmonic_domain, 1,
response_sigma_2)
data_domain = R.target
noiseless_data = R(mock_signal)
noise_amplitude = noiseless_data.val.std()/signal_to_noise
# Setting up the noise covariance and drawing a random noise realization
N = ift.ScalingOperator(noise_amplitude**2, data_domain)
noise = N.draw_sample()
data = noiseless_data + noise
# Wiener filter
j = R.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1)
sampling_ctrl = ift.GradientNormController(name="sampling",
tol_abs_gradnorm=1e2)
wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R, iteration_controller=ctrl,
iteration_controller_sampling=sampling_ctrl)
m_k = wiener_curvature.inverse_times(j)
m = ht(m_k)
plotdict = {"colormap": "Planck-like"}
plot_space = ift.RGSpace((N_pixels_1, N_pixels_2))
ift.plot(ht(mock_signal).cast_domain(plot_space),
name='mock_signal.png', **plotdict)
ift.plot(data.cast_domain(plot_space), name='data.png', **plotdict)
ift.plot(m.cast_domain(plot_space), name='map.png', **plotdict)
# sampling the uncertainty map
mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 50)
ift.plot(ift.sqrt(variance).cast_domain(plot_space),
name="uncertainty.png", **plotdict)
ift.plot((mean+m).cast_domain(plot_space),
name="posterior_mean.png", **plotdict)
import nifty5 as ift
import numpy as np
if __name__ == "__main__":
# Setting up parameters
L = 2. # Total side-length of the domain
N_pixels = 512 # Grid resolution (pixels per axis)
correlation_length_scale = .2 # 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
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
signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
harmonic_space = signal_space.get_default_codomain()
ht = ift.HarmonicTransformOperator(harmonic_space, target=signal_space)
# Creating the mock signal
S = ift.create_power_operator(harmonic_space,
power_spectrum=power_spectrum)
mock_signal = S.draw_sample()
# Setting up an exemplary response
mask = np.ones(signal_space.shape)
N10 = int(N_pixels/10)
mask[N10*5:N10*9, N10*5:N10*9] = 0.
mask = ift.Field.from_global_data(signal_space, mask).lock()
R = ift.GeometryRemover(signal_space)
R = R*ift.DiagonalOperator(mask)
R = R*ht
R = R * ift.create_harmonic_smoothing_operator((harmonic_space,), 0,
response_sigma)
data_domain = R.target[0]
noiseless_data = R(mock_signal)
noise_amplitude = noiseless_data.val.std()/signal_to_noise
# Setting up the noise covariance and drawing a random noise realization
N = ift.ScalingOperator(noise_amplitude**2, data_domain)
noise = N.draw_sample()
data = noiseless_data + noise
# Wiener filter
j = R.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=1e-2)
sampling_ctrl = ift.GradientNormController(name="sampling",
tol_abs_gradnorm=2e1)
wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R, iteration_controller=ctrl,
iteration_controller_sampling=sampling_ctrl)
m_k = wiener_curvature.inverse_times(j)
m = ht(m_k)
plotdict = {"colormap": "Planck-like"}
ift.plot(ht(mock_signal), name="mock_signal.png", **plotdict)
ift.plot(data.cast_domain(signal_space), name="data.png", **plotdict)
ift.plot(m, name="map.png", **plotdict)
# sampling the uncertainty map
mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 50)
ift.plot(ift.sqrt(variance), name="uncertainty.png", **plotdict)
ift.plot(mean+m, name="posterior_mean.png", **plotdict)
# Program to generate figures of article "Information theory for fields"
# by Torsten Ensslin, Annalen der Physik, submitted to special edition
# "Physics of Information" in April 2018, arXiv:1804.03350
# to get exact figure of paper comment out lines marked by "COMMENT OUT"
import numpy as np
import nifty5 as ift
import matplotlib.pyplot as plt
class Exp3(object):
def __call__(self, x):
return ift.exp(3*x)
def derivative(self, x):
return 3*ift.exp(3*x)
if __name__ == "__main__":
np.random.seed(20)
# Set up physical constants
nu = 1. # excitation field level
kappa = 10. # diffusion constant
eps = 1e-8 # small number to tame zero mode
sigma_n = 0.2 # noise level
sigma_n2 = sigma_n**2
L = 1. # Total length of interval or volume the field lives on
nprobes = 1000 # Number of probes for uncertainty quantification used in paper
nprobes = 10 # COMMENT OUT TO REPRODUCE PAPER FIGURE EXACTLY
# Define resolution (pixels per dimension)
N_pixels = 1024
# Define data gaps
N1a = int(0.6*N_pixels)
N1b = int(0.64*N_pixels)
N2a = int(0.67*N_pixels)
N2b = int(0.8*N_pixels)
# Set up derived constants
amp = nu/(2*kappa) # spectral normalization
pow_spec = lambda k: amp / (eps + k**2)
lambda2 = 2*kappa*sigma_n2/nu # resulting correlation length squared
lambda1 = np.sqrt(lambda2)
pixel_width = L/N_pixels
x = np.arange(0, 1, pixel_width)
# Set up the geometry
s_domain = ift.RGSpace([N_pixels], distances=pixel_width)
h_domain = s_domain.get_default_codomain()
HT = ift.HarmonicTransformOperator(h_domain, s_domain)
aHT = HT.adjoint
# Create mock signal
Phi_h = ift.create_power_operator(h_domain, power_spectrum=pow_spec)
phi_h = Phi_h.draw_sample()
# remove zero mode
glob = phi_h.to_global_data()
glob[0] = 0.
phi_h = ift.Field.from_global_data(phi_h.domain, glob)
phi = HT(phi_h)
# Setting up an exemplary response
GeoRem = ift.GeometryRemover(s_domain)
d_domain = GeoRem.target[0]
mask = np.ones(d_domain.shape)
mask[N1a:N1b] = 0.
mask[N2a:N2b] = 0.
fmask = ift.Field.from_global_data(d_domain, mask)
Mask = ift.DiagonalOperator(fmask)
R0 = Mask*GeoRem