Commit ee930b45 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'demos' into 'NIFTy_5'

Demos

See merge request ift/nifty-dev!23
parents a1f06b10 459a2770
......@@ -15,6 +15,8 @@ build_docker_from_scratch:
- schedules
image: docker:stable
stage: build_docker
before_script:
- ls
script:
- docker login -u gitlab-ci-token -p $CI_BUILD_TOKEN gitlab-registry.mpcdf.mpg.de
- docker build -t $CONTAINER_TEST_IMAGE --no-cache .
......@@ -25,6 +27,8 @@ build_docker_from_cache:
- schedules
image: docker:stable
stage: build_docker
before_script:
- ls
script:
- docker login -u gitlab-ci-token -p $CI_BUILD_TOKEN gitlab-registry.mpcdf.mpg.de
- docker build -t $CONTAINER_TEST_IMAGE .
......@@ -33,7 +37,6 @@ build_docker_from_cache:
test_python2_with_coverage:
stage: test
script:
- python setup.py install --user -f
- mpiexec -n 2 --bind-to none nosetests -q 2> /dev/null
- nosetests -q --with-coverage --cover-package=nifty5 --cover-erase
- >
......@@ -44,12 +47,13 @@ test_python2_with_coverage:
test_python3:
stage: test
script:
- python3 setup.py install --user -f
- nosetests3 -q
- mpiexec -n 2 --bind-to none nosetests3 -q 2> /dev/null
pages:
stage: release
before_script:
- ls
script:
- python setup.py install --user -f
- sh docs/generate.sh
......@@ -62,13 +66,31 @@ pages:
before_script:
- export MPLBACKEND="agg"
- python setup.py install --user -f
- python3 setup.py install --user -f
run_ipynb:
stage: demo_runs
script:
- python setup.py install --user -f
- python3 setup.py install --user -f
- jupyter nbconvert --execute --ExecutePreprocessor.timeout=None demos/Wiener_Filter.ipynb
artifacts:
paths:
- '*.png'
run_getting_started_1:
stage: demo_runs
script:
- python demos/getting_started_1.py
- python3 demos/getting_started_1.py
run_getting_started_2:
stage: demo_runs
script:
- python demos/getting_started_2.py
- python3 demos/getting_started_2.py
run_getting_started_3:
stage: demo_runs
script:
- python demos/getting_started_3.py
- python3 demos/getting_started_3.py
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 import GaussianEnergy
if __name__ == '__main__':
# s_space = ift.RGSpace([1024])
s_space = ift.RGSpace([128,128])
# s_space = ift.HPSpace(64)
def make_chess_mask():
mask = np.ones(position_space.shape)
for i in range(4):
for j in range(4):
if (i+j) % 2 == 0:
mask[i*128//4:(i+1)*128//4, j*128//4:(j+1)*128//4] = 0
return mask
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)
def make_random_mask():
mask = ift.from_random('pm1', position_space)
mask = (mask+1)/2
return mask.to_global_data()
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)
if __name__ == '__main__':
# # description of the tutorial ###
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)
# Choose problem geometry and masking
# One dimensional regular grid
position_space = ift.RGSpace([1024])
mask = np.ones(position_space.shape)
# # Two dimensional regular grid with chess mask
# position_space = ift.RGSpace([128,128])
# mask = make_chess_mask()
MOCK_POSITION = ift.from_random('normal',total_domain)
data = Rs.at(MOCK_POSITION).value + N.draw_sample()
# # Sphere with half of its locations randomly masked
# position_space = ift.HPSpace(128)
# mask = make_random_mask()
NWR = ApplyData(data, ift.Field(d_space,val=noise), Rs)
harmonic_space = position_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space)
INITIAL_POSITION = ift.from_random('normal',total_domain)
likelihood = GaussianEnergy(INITIAL_POSITION, NWR)
# set correlation structure with a power spectrum and build
# prior correlation covariance
def power_spectrum(k):
return 100. / (20.+k**3)
power_space = ift.PowerSpace(harmonic_space)
PD = ift.PowerDistributor(harmonic_space, power_space)
prior_correlation_structure = PD(ift.PS_field(power_space, power_spectrum))
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)
S = ift.DiagonalOperator(prior_correlation_structure)
# build instrument response consisting of a discretization, mask
# and harmonic transformaion
GR = ift.GeometryRemover(position_space)
mask = ift.Field.from_global_data(position_space, mask)
Mask = ift.DiagonalOperator(mask)
R = GR * Mask * HT
H = Hamiltonian(likelihood, inverter)
H, convergence = minimizer(H)
result = s.at(H.position).value
data_space = GR.target
# setting the noise covariance
noise = 5.
N = ift.ScalingOperator(noise, data_space)
# creating mock data
MOCK_SIGNAL = S.draw_sample()
MOCK_NOISE = N.draw_sample()
data = R(MOCK_SIGNAL) + MOCK_NOISE
# building propagator D and information source j
j = R.adjoint_times(N.inverse_times(data))
D_inv = R.adjoint * N.inverse * R + S.inverse
# make it invertible
IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3)
D = ift.InversionEnabler(D_inv, IC, approximation=S.inverse).inverse
# WIENER FILTER
m = D(j)
# PLOTTING
# Truth, data, reconstruction, residuals
import nifty5 as ift
import numpy as np
def get_2D_exposure():
x_shape, y_shape = position_space.shape
exposure = np.ones(position_space.shape)
exposure[x_shape//3:x_shape//2, :] *= 2.
exposure[x_shape*4//5:x_shape, :] *= .1
exposure[x_shape//2:x_shape*3//2, :] *= 3.
exposure[:, x_shape//3:x_shape//2] *= 2.
exposure[:, x_shape*4//5:x_shape] *= .1
exposure[:, x_shape//2:x_shape*3//2] *= 3.
exposure = ift.Field.from_global_data(position_space, exposure)
return exposure
if __name__ == '__main__':
# ABOUT THIS CODE
np.random.seed(41)
# Set up the position space of the signal
#
# # One dimensional regular grid with uniform exposure
# position_space = ift.RGSpace(1024)
# exposure = np.ones(position_space.shape)
# Two-dimensional regular grid with inhomogeneous exposure
position_space = ift.RGSpace([512, 512])
exposure = get_2D_exposure()
# # Sphere with with uniform exposure
# position_space = ift.HPSpace(128)
# exposure = ift.Field.full(position_space, 1.)
# Defining harmonic space and transform
harmonic_space = position_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
domain = ift.MultiDomain.make({'xi': harmonic_space})
position = ift.from_random('normal', domain)
# Define power spectrum and amplitudes
def sqrtpspec(k):
return 1. / (20. + k**2)
p_space = ift.PowerSpace(harmonic_space)
pd = ift.PowerDistributor(harmonic_space, p_space)
a = ift.PS_field(p_space, sqrtpspec)
A = pd(a)
# Set up a sky model
xi = ift.Variable(position)['xi']
logsky_h = xi * A
logsky = HT(logsky_h)
sky = ift.PointwiseExponential(logsky)
M = ift.DiagonalOperator(exposure)
GR = ift.GeometryRemover(position_space)
# Set up instrumental response
R = GR * M
# Generate mock data
d_space = R.target[0]
lamb = R(sky)
mock_position = ift.from_random('normal', lamb.position.domain)
data = lamb.at(mock_position).value
data = np.random.poisson(data.to_global_data().astype(np.float64))
data = ift.Field.from_global_data(d_space, data)
# Compute likelihood and Hamiltonian
position = ift.from_random('normal', lamb.position.domain)
likelihood = ift.PoissonianEnergy(lamb, data)
ic_cg = ift.GradientNormController(iteration_limit=50)
ic_newton = ift.GradientNormController(name='Newton', iteration_limit=50,
tol_abs_gradnorm=1e-3)
minimizer = ift.RelaxedNewton(ic_newton)
# Minimize the Hamiltonian
H = ift.Hamiltonian(likelihood, ic_cg)
H, convergence = minimizer(H)
# Plot results
result_sky = sky.at(H.position).value
# PLOTTING
import nifty5 as ift
from nifty5.library.los_response import LOSResponse
from nifty5.library.amplitude_model import make_amplitude_model
from nifty5.library.smooth_sky import make_correlated_field
import numpy as np
from scipy.io import loadmat
def get_random_LOS(n_los):
starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
ends = list(np.random.uniform(0, 1, (n_los, 2)).T)
return starts, ends
if __name__ == '__main__':
# ## ABOUT THIS TUTORIAL
np.random.seed(42)
position_space = ift.RGSpace([128, 128])
# Setting up an amplitude model
A, amplitude_internals = make_amplitude_model(
position_space, 16, 1, 10, -4., 1, 0., 1.)
# Building the model for a correlated signal
harmonic_space = position_space.get_default_codomain()
ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
power_space = A.value.domain[0]
power_distributor = ift.PowerDistributor(harmonic_space, power_space)
position = {}
position['xi'] = ift.Field.from_random('normal', harmonic_space)
position = ift.MultiField(position)
xi = ift.Variable(position)['xi']
Amp = power_distributor(A)
correlated_field_h = Amp * xi
correlated_field = ht(correlated_field_h)
# # alternatively to the block above one can do:
# correlated_field, _ = make_correlated_field(position_space, A)
# apply some nonlinearity
signal = ift.PointwisePositiveTanh(correlated_field)
# Building the Line of Sight response
LOS_starts, LOS_ends = get_random_LOS(100)
R = LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends)
# build signal response model and model likelihood
signal_response = R(signal)
# specify noise
data_space = R.target
noise = .001
N = ift.ScalingOperator(noise, data_space)
# generate mock data
MOCK_POSITION = ift.from_random('normal', signal.position.domain)
data = signal_response.at(MOCK_POSITION).value + N.draw_sample()
# set up model likelihood
likelihood = ift.GaussianEnergy(signal_response, mean=data, covariance=N)
# set up minimization and inversion schemes
ic_cg = ift.GradientNormController(iteration_limit=10)
ic_sampling = ift.GradientNormController(iteration_limit=100)
ic_newton = ift.GradientNormController(name='Newton', iteration_limit=100)
minimizer = ift.RelaxedNewton(ic_newton)
# build model Hamiltonian
H = ift.Hamiltonian(likelihood, ic_cg,
iteration_controller_sampling=ic_sampling)
INITIAL_POSITION = ift.from_random('normal', H.position.domain)
position = INITIAL_POSITION
ift.plot(signal.at(MOCK_POSITION).value, name='truth.pdf')
ift.plot(R.adjoint_times(data), name='data.pdf')
ift.plot([A.at(MOCK_POSITION).value], name='power.pdf')
# number of samples used to estimate the KL
N_samples = 20
for i in range(5):
H = H.at(position)
samples = [H.curvature.draw_sample(from_inverse=True)
for _ in range(N_samples)]
KL = ift.SampledKullbachLeiblerDivergence(H, samples, ic_cg)
KL, convergence = minimizer(KL)
position = KL.position
ift.plot(signal.at(position).value, name='reconstruction.pdf')
ift.plot([A.at(position).value, A.at(MOCK_POSITION).value],
name='power.pdf')
avrg = 0.
va = 0.
powers = []
for sample in samples:
sam = signal.at(sample + position).value
powers.append(A.at(sample+position).value)
avrg += sam
va += sam**2
avrg /= len(samples)
va /= len(samples)
va -= avrg**2
std = ift.sqrt(va)
ift.plot(avrg, name='avrg.pdf')
ift.plot(std, name='std.pdf')
ift.plot([A.at(position).value, A.at(MOCK_POSITION).value]+powers,
name='power.pdf')
......@@ -4,6 +4,5 @@ from .gaussian_energy import GaussianEnergy
from .los_response import LOSResponse
from .point_sources import PointSources
from .poissonian_energy import PoissonianEnergy
from .smooth_sky import make_smooth_mf_sky_model, make_smooth_sky_model
from .wiener_filter_curvature import WienerFilterCurvature
from .wiener_filter_energy import WienerFilterEnergy
......@@ -22,7 +22,7 @@ def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
Npixdof : #pix in dof_space
ceps_a, ceps_k0 : Smoothnessparameters in ceps_kernel
ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel
eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
a = ceps_a, k0 = ceps_k0
......@@ -48,7 +48,7 @@ def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
phi_sig = np.array([sv, iv])
slope = SlopeOperator(param_space, logk_space, phi_sig)
norm_phi_mean = Field(param_space, val=phi_mean/phi_sig)
norm_phi_mean = Field.from_global_data(param_space, phi_mean/phi_sig)
fields = {keys[0]: Field.from_random('normal', dof_space),
keys[1]: Field.from_random('normal', param_space)}
......@@ -93,7 +93,7 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
# Prepare q_array
q_array = np.zeros((dim,) + shape)
if dim == 1:
ks = domain.get_k_length_array().val
ks = domain.get_k_length_array().to_global_data()
q_array = np.array([ks])
else:
for i in range(dim):
......@@ -120,4 +120,4 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
# Do summation
cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
return Field(domain, val=cepstrum_field)
return Field.from_global_data(domain, cepstrum_field)
......@@ -21,40 +21,40 @@ class PointSources(Model):
@property
@memo
def value(self):
points = self.position['points'].to_global_data()
points = self.position['points'].local_data
points = np.clip(points, None, 8.2)
points = Field(self.position['points'].domain, points)
points = Field.from_local_data(self.position['points'].domain, points)
return self.IG(points, self._alpha, self._q)
@property
@memo
def gradient(self):
u = self.position['points']
inner = norm.pdf(u.val)
outer_inv = invgamma.pdf(invgamma.ppf(norm.cdf(u.val),
u = self.position['points'].local_data
inner = norm.pdf(u)
outer_inv = invgamma.pdf(invgamma.ppf(norm.cdf(u),
self._alpha,
scale=self._q),
self._alpha, scale=self._q)
# FIXME
outer_inv = np.clip(outer_inv, 1e-20, None)
outer = 1/outer_inv
grad = Field(u.domain, val=inner*outer)
grad = Field.from_local_data(u.domain, inner*outer)
grad = makeOp(MultiField({'points': grad}))
return SelectionOperator(grad.target, 'points')*grad
@staticmethod
def IG(field, alpha, q):
foo = invgamma.ppf(norm.cdf(field.val), alpha, scale=q)
return Field(field.domain, val=foo)
foo = invgamma.ppf(norm.cdf(field.local_data), alpha, scale=q)
return Field.from_local_data(field.domain, foo)
@staticmethod
def IG_prime(field, alpha, q):
inner = norm.pdf(field.val)
outer = invgamma.pdf(invgamma.ppf(norm.cdf(field.val), alpha, scale=q), alpha, scale=q)
inner = norm.pdf(field.local_data)
outer = invgamma.pdf(invgamma.ppf(norm.cdf(field.local_data), alpha, scale=q), alpha, scale=q)
# # FIXME
# outer = np.clip(outer, 1e-20, None)
outer = 1/outer
return Field(field.domain, val=inner*outer)
return Field.from_local_data(field.domain, inner*outer)
@staticmethod
def inverseIG(u, alpha, q):
......
def make_smooth_sky_model(s_space, amplitude_model):
def make_correlated_field(s_space, amplitude_model):
'''
Method for construction of correlated sky model
......@@ -22,12 +22,12 @@ def make_smooth_sky_model(s_space, amplitude_model):
xi = Variable(position)['xi']
A = power_distributor(amplitude_model)
logsky_h = A * xi
logsky = ht(logsky_h)
internals = {'logsky_h': logsky_h,
correlated_field_h = A * xi
correlated_field = ht(correlated_field_h)
internals = {'correlated_field_h': correlated_field_h,
'power_distributor': power_distributor,
'ht': ht}
return PointwiseExponential(logsky), internals
return correlated_field, internals
def make_smooth_mf_sky_model(s_space_spatial, s_space_energy,
......
......@@ -67,7 +67,7 @@ class MultiField(object):
dtype = MultiField.build_dtype(dtype, domain)
return MultiField({key: Field.from_random(random_type, domain[key],
dtype[key], **kwargs)
for key in domain.keys()})
for key in sorted(domain.keys())})
def fill(self, fill_value):
"""Fill `self` uniformly with `fill_value`
......
......@@ -67,7 +67,7 @@ class ExpTransform(LinearOperator):
def apply(self, x, mode):
self._check_input(x, mode)
x = x.val
x = x.to_global_data()
ndim = len(self.target.shape)
idx = ()
for d in range(ndim):
......@@ -90,7 +90,7 @@ class ExpTransform(LinearOperator):
x = xnew
idx = (slice(None),) + idx
return Field(self._tgt(mode), val=x)
return Field.from_global_data(self._tgt(mode), x)
@property
def capability(self):
......
from ..domain_tuple import DomainTuple
from ..field import Field
from .. import dobj
from ..utilities import hartley
from .linear_operator import LinearOperator
......@@ -34,11 +35,13 @@ class QHTOperator(LinearOperator):
x = x.val * self.domain[0].scalar_dvol()
n = len(self.domain[0].shape)
rng = range(n) if mode == self.TIMES else reversed(range(n))
# MR FIXME: this needs to be fixed properly for MPI
x = dobj.to_global_data(x)
for i in rng:
sl = (slice(None),)*i + (slice(1, None),)
x[sl] = hartley(x[sl], axes=(i,))
return Field(self._tgt(mode), val=x)
return Field.from_global_data(self._tgt(mode), x)
@property
def capability(self):
......
......@@ -18,7 +18,7 @@ class SlopeOperator(LinearOperator):
self.pos = np.zeros((self.ndim,) + self.target[0].shape)
if self.ndim == 1:
self.pos[0] = self.target[0].get_k_length_array().val
self.pos[0] = self.target[0].get_k_length_array().to_global_data()
else:
shape = self.target[0].shape
for i in range(self.ndim):
......@@ -46,18 +46,19 @@ class SlopeOperator(LinearOperator):
# Times
if mode == self.TIMES:
inp = x.val
inp = x.to_global_data()
res = self.sigmas[-1] * inp[-1]
for i in range(self.ndim):
res += self.sigmas[i] * inp[i] * self.pos[i]
return Field(self.target, val=res)
return Field.from_global_data(self.target, res)
# Adjoint times
res = np.zeros(self.domain[0].shape)
res[-1] = np.sum(x.val) * self.sigmas[-1]
xglob = x.to_global_data()
res[-1] = np.sum(xglob) * self.sigmas[-1]
for i in range(self.ndim):
res[i] = np.sum(self.pos[i] * x.val) * self.sigmas[i]
return Field(self.domain, val=res)
res[i] = np.sum(self.pos[i] * xglob) * self.sigmas[i]
return Field.from_global_data(self.domain, res)
@property
def capability(self):
......
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