Commit e53f2d33 authored by Philipp Arras's avatar Philipp Arras

Merge branch 'migrate_from_gn' into 'NIFTy_5'

Migrate from gn

See merge request ift/nifty-dev!3
parents 52e8355f bb223a89
......@@ -45,8 +45,8 @@ test_python3:
stage: test
- python3 install --user -f
- mpiexec -n 2 --bind-to none nosetests3 -q 2> /dev/null
- nosetests3 -q
- mpiexec -n 2 --bind-to none nosetests3 -q 2> /dev/null
stage: release
import nifty5 as ift
from nifty5.library.nonlinearities import Linear, Tanh, Exponential
import numpy as np
from nifty5 import Exponential, Linear, Tanh
import nifty5 as ift
from nifty5.library.nonlinearities import Exponential
from nifty5 import Exponential
import numpy as np
......@@ -4,6 +4,9 @@ from . import dobj
from .domains import *
from .domain_tuple import DomainTuple
from .field import Field
from .nonlinearities import Exponential, Linear, PositiveTanh, Tanh
from .models import *
from .operators import *
from .probing.utils import probe_with_posterior_samples, probe_diagonal, \
......@@ -12,6 +15,7 @@ from .minimization import *
from .sugar import *
from .plotting.plot import plot
from . import library
from . import extra
......@@ -21,5 +25,7 @@ from .logger import logger
from .multi import *
from .energies import *
# We deliberately don't set __all__ here, because we don't want people to do a
# "from nifty5 import *"; that would swamp the global namespace.
......@@ -7,6 +7,8 @@ from .hp_space import HPSpace
from .gl_space import GLSpace
from .dof_space import DOFSpace
from .power_space import PowerSpace
from .log_rg_space import LogRGSpace
__all__ = ["Domain", "UnstructuredDomain", "StructuredDomain", "RGSpace",
"LMSpace", "HPSpace", "GLSpace", "DOFSpace", "PowerSpace"]
"LMSpace", "HPSpace", "GLSpace", "DOFSpace", "PowerSpace",
from functools import reduce
from ..sugar import exp
import numpy as np
from ..dobj import ibegin
from ..field import Field
from .structured_domain import StructuredDomain
class LogRGSpace(StructuredDomain):
_needed_for_hash = ['_shape', '_bindistances', '_t_0', '_harmonic']
def __init__(self, shape, bindistances, t_0, harmonic=False):
super(LogRGSpace, self).__init__()
self._harmonic = bool(harmonic)
if np.isscalar(shape):
shape = (shape,)
self._shape = tuple(int(i) for i in shape)
self._bindistances = tuple(bindistances)
self._t_0 = tuple(t_0)
self._dim = int(reduce(lambda x, y: x * y, self._shape))
self._dvol = float(reduce(lambda x, y: x * y, self._bindistances))
def harmonic(self):
return self._harmonic
def shape(self):
return self._shape
def scalar_dvol(self):
return self._dvol
def bindistances(self):
return np.array(self._bindistances)
def size(self):
def t_0(self):
return np.array(self._t_0)
def __repr__(self):
return ("LogRGSpace(shape=%r, harmonic=%r)"
% (self.shape, self.harmonic))
def get_default_codomain(self):
if self._harmonic:
raise ValueError("only supported for nonharmonic space")
codomain_bindistances = 1. / (self.bindistances * self.shape)
return LogRGSpace(self.shape, codomain_bindistances,
np.zeros(len(self.shape)), True)
def get_k_length_array(self):
out = Field(self, dtype=np.float64)
oloc = out.local_data
ib = ibegin(out.val)
res = np.arange(oloc.shape[0], dtype=np.float64) + ib[0]
res = np.minimum(res, self.shape[0]-res)*self.bindistances[0]
if len(self.shape) == 1:
oloc[()] = res
return out
res *= res
for i in range(1, len(self.shape)):
tmp = np.arange(oloc.shape[i], dtype=np.float64) + ib[i]
tmp = np.minimum(tmp, self.shape[i]-tmp)*self.bindistances[i]
tmp *= tmp
res = np.add.outer(res, tmp)
oloc[()] = np.sqrt(res)
return out
def get_expk_length_array(self):
# FIXME This is a hack! Only for plotting. Seems not to be the final version.
out = exp(self.get_k_length_array())
out.val[1:] = out.val[:-1]
out.val[0] = 0
return out
from .hamiltonian import Hamiltonian
from .kl import SampledKullbachLeiblerDivergence
from nifty5 import Energy, InversionEnabler, SamplingEnabler, Variable, memo
from nifty5.library import UnitLogGauss
class Hamiltonian(Energy):
def __init__(self, lh, iteration_controller,
lh: Likelihood (energy object)
super(Hamiltonian, self).__init__(lh.position)
self._lh = lh
self._ic = iteration_controller
if iteration_controller_sampling is None:
self._ic_samp = iteration_controller
self._ic_samp = iteration_controller_sampling
self._prior = UnitLogGauss(Variable(self.position))
self._precond = self._prior.curvature
def at(self, position):
return self.__class__(, self._ic, self._ic_samp)
def value(self):
return self._lh.value + self._prior.value
def gradient(self):
return self._lh.gradient + self._prior.gradient
def curvature(self):
prior_curv = self._prior.curvature
c = SamplingEnabler(self._lh.curvature, prior_curv.inverse,
self._ic_samp, prior_curv.inverse)
return InversionEnabler(c, self._ic, self._precond)
def __str__(self):
res = 'Likelihood:\t{:.2E}\n'.format(self._lh.value)
res += 'Prior:\t\t{:.2E}'.format(self._prior.value)
return res
from nifty5 import Energy, InversionEnabler, ScalingOperator, memo
class SampledKullbachLeiblerDivergence(Energy):
def __init__(self, h, res_samples, iteration_controller):
h: Hamiltonian
N: Number of samples to be used
super(SampledKullbachLeiblerDivergence, self).__init__(h.position)
self._h = h
self._res_samples = res_samples
self._iteration_controller = iteration_controller
self._energy_list = []
for ss in res_samples:
e =
def at(self, position):
return self.__class__(, self._res_samples,
def value(self):
v = self._energy_list[0].value
for energy in self._energy_list[1:]:
v += energy.value
return v / len(self._energy_list)
def gradient(self):
g = self._energy_list[0].gradient
for energy in self._energy_list[1:]:
g += energy.gradient
return g / len(self._energy_list)
def curvature(self):
# MR FIXME: This looks a bit strange...
approx = self._energy_list[-1]._prior.curvature
curvature_list = [e.curvature for e in self._energy_list]
op = curvature_list[0]
for curv in curvature_list[1:]:
op = op + curv
op = op * ScalingOperator(1./len(curvature_list), op.domain)
return InversionEnabler(op, self._iteration_controller, approx)
from .amplitude_model import make_amplitude_model
from .apply_data import ApplyData
from .los_response import LOSResponse
from .noise_energy import NoiseEnergy
from .nonlinear_power_energy import NonlinearPowerEnergy
from .nonlinear_wiener_filter_energy import NonlinearWienerFilterEnergy
from .nonlinearities import Exponential, Linear, PositiveTanh, Tanh
from .poisson_energy import PoissonEnergy
from .unit_log_gauss import UnitLogGauss
from .point_sources import PointSources
from .poisson_log_likelihood import PoissonLogLikelihood
from .smooth_sky import make_smooth_mf_sky_model, make_smooth_sky_model
from .unit_log_gauss import UnitLogGauss
from .wiener_filter_curvature import WienerFilterCurvature
from .wiener_filter_energy import WienerFilterEnergy
import numpy as np
from import PowerSpace, UnstructuredDomain
from ..field import Field
from ..multi import MultiField
from ..sugar import makeOp, sqrt
def _ceps_kernel(dof_space, k, a, k0):
return a**2/(1+(k/(k0*dof_space.bindistances[0]))**2)**2
def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
keys=['tau', 'phi']):
Method for construction of amplitude model
Computes a smooth power spectrum.
Output lives in PowerSpace.
Npixdof : #pix in dof_space
ceps_a, ceps_k0 : Smoothnessparameters in ceps_kernel
eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
a = ceps_a, k0 = ceps_k0
sm, sv : slope_mean = expected exponent of powerlaw (e.g. -4),
slope_variance (default=1)
im, iv : y-intercept_mean, y-intercept_variance of power_slope
from ..operators import (ExpTransform, QHTOperator, SlopeOperator,
h_space = s_space.get_default_codomain()
p_space = PowerSpace(h_space)
exp_transform = ExpTransform(p_space, Npixdof)
logk_space = exp_transform.domain[0]
dof_space = logk_space.get_default_codomain()
param_space = UnstructuredDomain(2)
qht = QHTOperator(dof_space, logk_space)
sym = SymmetrizingOperator(logk_space)
phi_mean = np.array([sm, im])
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)
fields = {keys[0]: Field.from_random('normal', dof_space),
keys[1]: Field.from_random('normal', param_space)}
position = MultiField(fields)
dof_space = position[keys[0]].domain[0]
kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
cepstrum = create_cepstrum_amplitude_field(dof_space, kern)
ceps = makeOp(sqrt(cepstrum))
smooth_op = sym * qht * ceps
from ..models import Variable
smooth_spec = smooth_op(Variable(position)[keys[0]])
phi = Variable(position)[keys[1]] + norm_phi_mean
linear_spec = slope(phi)
loglog_spec = smooth_spec + linear_spec
from ..models import PointwiseExponential
xlog_ampl = PointwiseExponential(0.5*loglog_spec)
internals = {'loglog_spec': loglog_spec,
'qht': qht,
'ceps': ceps,
'norm_phi_mean': norm_phi_mean}
return exp_transform(xlog_ampl), internals
def create_cepstrum_amplitude_field(domain, cepstrum):
"""Creates a ...
Writes the sum of all modes into the zero-mode.
domain: ???
cepstrum: Callable
dim = len(domain.shape)
dist = domain.bindistances
shape = domain.shape
# Prepare q_array
q_array = np.zeros((dim,) + shape)
if dim == 1:
ks = domain.get_k_length_array().val
q_array = np.array([ks])
for i in range(dim):
ks = np.minimum(shape[i] - np.arange(shape[i]) +
1, np.arange(shape[i])) * dist[i]
fst_dims = (1,) * i
lst_dims = (1,) * (dim - i - 1)
q_array[i] += ks.reshape(fst_dims + (shape[i],) + lst_dims)
# Fill cepstrum field (all non-zero modes)
no_zero_modes = (slice(1, None),) * dim
ks = q_array[(slice(None),) + no_zero_modes]
cepstrum_field = np.zeros(shape)
cepstrum_field[no_zero_modes] = cepstrum(ks)
# Fill cepstrum field (zero-mode subspaces)
for i in range(dim):
# Prepare indices
fst_dims = (slice(None),) * i
lst_dims = (slice(None),) * (dim - i - 1)
sl = fst_dims + (slice(1, None),) + lst_dims
sl2 = fst_dims + (0,) + lst_dims
# Do summation
cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
return Field(domain, val=cepstrum_field)
def ApplyData(data, var, model_data):
from .. import DiagonalOperator, Constant, sqrt
sqrt_n = DiagonalOperator(sqrt(var))
data = Constant(model_data.position, data)
return sqrt_n.inverse(model_data - data)
import numpy as np
from scipy.stats import invgamma, norm
from ..field import Field
from ..sugar import makeOp
from ..multi import MultiField
from ..models import Model
from ..operators import SelectionOperator
from ..utilities import memo
class PointSources(Model):
def __init__(self, position, alpha, q):
super(PointSources, self).__init__(position)
self._alpha = alpha
self._q = q
def at(self, position):
return self.__class__(position, self._alpha, self._q)
def value(self):
points = self.position['points'].to_global_data()
points = np.clip(points, None, 8.2)
points = Field(self.position['points'].domain, points)
return self.IG(points, self._alpha, self._q)
def gradient(self):
u = self.position['points']
inner = norm.pdf(u.val)
outer_inv = invgamma.pdf(invgamma.ppf(norm.cdf(u.val),
self._alpha, scale=self._q)
outer_inv = np.clip(outer_inv, 1e-20, None)
outer = 1/outer_inv
grad = Field(u.domain, val=inner*outer)
grad = makeOp(MultiField({'points': grad}))
return SelectionOperator(, 'points')*grad
def IG(field, alpha, q):
foo = invgamma.ppf(norm.cdf(field.val), alpha, scale=q)
return Field(field.domain, val=foo)
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)
# outer = np.clip(outer, 1e-20, None)
outer = 1/outer
return Field(field.domain, val=inner*outer)
def inverseIG(u, alpha, q):
res = norm.ppf(invgamma.cdf(u, alpha, scale=q))
# res = np.clip(res, 0, None)
return res
def make_smooth_sky_model(s_space, amplitude_model):
Method for construction of correlated sky model
s_space : domain of sky model
amplitude_model : model for correlation structure
from .. import (FFTOperator, Field, MultiField, PointwiseExponential,
PowerDistributor, Variable)
h_space = s_space.get_default_codomain()
ht = FFTOperator(h_space, s_space)
p_space = amplitude_model.value.domain[0]
power_distributor = PowerDistributor(h_space, p_space)
position = {}
position['xi'] = Field.from_random('normal', h_space)
position['tau'] = amplitude_model.position['tau']
position['phi'] = amplitude_model.position['phi']
position = MultiField(position)
xi = Variable(position)['xi']
A = power_distributor(amplitude_model)
logsky_h = A * xi
logsky = ht(logsky_h)
internals = {'logsky_h': logsky_h,
'power_distributor': power_distributor,
'ht': ht}
return PointwiseExponential(logsky), internals
def make_smooth_mf_sky_model(s_space_spatial, s_space_energy,
amplitude_model_spatial, amplitude_model_energy):
Method for construction of correlated sky model
s_space : domain of sky model
amplitude_model : model for correlation structure
from .. import (DomainTuple, Field, MultiField,
PointwiseExponential, Variable)
from ..operators import (DomainDistributor, PowerDistributor,
h_space_spatial = s_space_spatial.get_default_codomain()
h_space_energy = s_space_energy.get_default_codomain()
h_space = DomainTuple.make((h_space_spatial, h_space_energy))
ht1 = HarmonicTransformOperator(h_space, space=0)
ht2 = HarmonicTransformOperator(, space=1)
ht = ht2*ht1
p_space_spatial = amplitude_model_spatial.value.domain[0]
p_space_energy = amplitude_model_energy.value.domain[0]
pd_spatial = PowerDistributor(h_space, p_space_spatial, 0)
pd_energy = PowerDistributor(pd_spatial.domain, p_space_energy, 1)
pd = pd_spatial*pd_energy
dom_distr_0 = DomainDistributor(pd.domain, 0)
dom_distr_1 = DomainDistributor(pd.domain, 1)
a_spatial = dom_distr_1(amplitude_model_spatial)
a_energy = dom_distr_0(amplitude_model_energy)
a = a_spatial*a_energy
A = pd(a)
position = MultiField({'xi': Field.from_random('normal', h_space)})
xi = Variable(position)['xi']
logsky_h = A*xi
logsky = ht(logsky_h)
return PointwiseExponential(logsky)
......@@ -16,7 +16,7 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from ..library.nonlinearities import Exponential, PositiveTanh, Tanh
from ..nonlinearities import Exponential, PositiveTanh, Tanh
from ..sugar import makeOp
from .model import Model
......@@ -16,7 +16,7 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from ..sugar import full, exp, tanh
from .sugar import full, exp, tanh
class Linear(object):
from .diagonal_operator import DiagonalOperator
from .dof_distributor import DOFDistributor
from .domain_distributor import DomainDistributor
from .endomorphic_operator import EndomorphicOperator
from .exp_transform import ExpTransform
from .fft_operator import FFTOperator
from .fft_smoothing_operator import FFTSmoothingOperator
from .geometry_remover import GeometryRemover
......@@ -10,15 +12,20 @@ from .laplace_operator import LaplaceOperator
from .linear_operator import LinearOperator
from .multi_adaptor import MultiAdaptor
from .power_distributor import PowerDistributor
from .qht_operator import QHTOperator
from .sampling_enabler import SamplingEnabler
from .sandwich_operator import SandwichOperator
from .scaling_operator import ScalingOperator
from .selection_operator import SelectionOperator
from .slope_operator import SlopeOperator
from .smoothness_operator import SmoothnessOperator
from .symmetrizing_operator import SymmetrizingOperator
__all__ = ["LinearOperator", "EndomorphicOperator", "ScalingOperator",
"DiagonalOperator", "HarmonicTransformOperator", "FFTOperator",
"FFTSmoothingOperator", "GeometryRemover",
"LaplaceOperator", "SmoothnessOperator", "PowerDistributor",
"InversionEnabler", "SandwichOperator", "SamplingEnabler",
"DOFDistributor", "SelectionOperator", "MultiAdaptor"]
"DOFDistributor", "SelectionOperator", "MultiAdaptor",
"ExpTransform", "SymmetrizingOperator", "QHTOperator",
"SlopeOperator", "DomainDistributor"]
import numpy as np
from ..field import Field
from .. import dobj
from ..domain_tuple import DomainTuple
from .linear_operator import LinearOperator
class DomainDistributor(LinearOperator):
def __init__(self, target, axis):
# TODO Replace this by a DiagonalOperator
if dobj.ntask > 1:
raise NotImplementedError('UpProj class does not support MPI.')
assert len(target) == 2
assert axis in [0, 1]
if axis == 0:
domain = target[1]
self._size = target[0].size
domain = target[0]
self._size = target[1].size
self._axis = axis
self._domain = DomainTuple.make(domain)
self._target = DomainTuple.make(target)
def domain(self):
return self._domain
def target(self):
return self._target
def apply(<