Commit 4e45b6a2 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'NIFTy_5' into operator_documentation

parents 5244664a 52593191
......@@ -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')
......@@ -59,6 +59,8 @@ class data_object(object):
distaxis = -1
self._distaxis = distaxis
self._data = data
if local_shape(self._shape, self._distaxis) != self._data.shape:
raise ValueError("shape mismatch")
# def _sanity_checks(self):
# # check whether the distaxis is consistent
......
......@@ -79,7 +79,7 @@ def distaxis(arr):
def from_local_data(shape, arr, distaxis=-1):
if shape != arr.shape:
if tuple(shape) != arr.shape:
raise ValueError
return arr
......
......@@ -21,7 +21,17 @@ from .structured_domain import StructuredDomain
class DOFSpace(StructuredDomain):
"""Generic degree-of-freedom space."""
"""Generic degree-of-freedom space. It is defined as the domain of some
DOFDistributor.
Its entries represent the underlying degrees of freedom of some other
space, according to the dofdex.
Parameters
----------
dof_weights: 1-D numpy array
A numpy array containing the multiplicity of each individual degree of
freedom.
"""
_needed_for_hash = ["_dvol"]
......
......@@ -24,7 +24,6 @@ class SampledKullbachLeiblerDivergence(Energy):
return self.__class__(self._h.at(position), self._res_samples,
self._iteration_controller)
@property
@memo
def value(self):
......
......@@ -4,6 +4,6 @@ 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
from .smooth_sky import make_correlated_field, make_mf_correlated_field
......@@ -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
Method for construction of correlated fields
Parameters
----------
s_space : domain of sky model
s_space : Field domain
amplitude_model : model for correlation structure
'''
......@@ -22,24 +22,18 @@ 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,
def make_mf_correlated_field(s_space_spatial, s_space_energy,
amplitude_model_spatial, amplitude_model_energy):
'''
Method for construction of correlated sky model
Parameters
----------
s_space : domain of sky model
amplitude_model : model for correlation structure
Method for construction of correlated multi-frequency fields
'''
from .. import (DomainTuple, Field, MultiField,
PointwiseExponential, Variable)
......
......@@ -21,9 +21,9 @@ from ..sugar import makeOp
from .model import Model
def _joint_position(op1, op2):
a = op1.position._val
b = op2.position._val
def _joint_position(model1, model2):
a = model1.position._val
b = model2.position._val
# Note: In python >3.5 one could do {**a, **b}
ab = a.copy()
ab.update(b)
......@@ -31,57 +31,78 @@ def _joint_position(op1, op2):
class ScalarMul(Model):
def __init__(self, factor, op):
# TODO op -> model
super(ScalarMul, self).__init__(op.position)
"""Class representing a model multiplied by a scalar factor."""
def __init__(self, factor, model):
super(ScalarMul, self).__init__(model.position)
# TODO -> floating