Commit e288370f authored by Philipp Arras's avatar Philipp Arras

Merge branch 'NIFTy_5' into demos

parents 6a5e113f b6c358a0
......@@ -45,8 +45,8 @@ test_python3:
stage: test
script:
- python3 setup.py 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
pages:
stage: release
......
......@@ -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",
"LogRGSpace"]
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))
@property
def harmonic(self):
return self._harmonic
@property
def shape(self):
return self._shape
def scalar_dvol(self):
return self._dvol
@property
def bindistances(self):
return np.array(self._bindistances)
@property
def size(self):
return np.prod(self._shape)
@property
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 .kl import SampledKullbachLeiblerDivergence
from .hamiltonian import Hamiltonian
from ..minimization.energy import Energy
from ..operators import InversionEnabler, SamplingEnabler
from ..models.variable import Variable
from ..utilities import memo
from ..library.unit_log_gauss import UnitLogGauss
class Hamiltonian(Energy):
def __init__(self, lh, iteration_controller,
iteration_controller_sampling=None):
"""
lh: Likelihood (energy object)
prior:
"""
super(Hamiltonian, self).__init__(lh.position)
self._lh = lh
self._ic = iteration_controller
if iteration_controller_sampling is None:
self._ic_samp = iteration_controller
else:
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._lh.at(position), self._ic, self._ic_samp)
@property
@memo
def value(self):
return self._lh.value + self._prior.value
@property
@memo
def gradient(self):
return self._lh.gradient + self._prior.gradient
@property
@memo
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 ..minimization.energy import Energy
from ..operators.inversion_enabler import InversionEnabler
from ..operators.scaling_operator import ScalingOperator
from ..utilities import 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 = h.at(self.position+ss)
self._energy_list.append(e)
def at(self, position):
return self.__class__(self._h.at(position), self._res_samples,
self._iteration_controller)
@property
@memo
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)
@property
@memo
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)
@property
@memo
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)
......@@ -487,6 +487,16 @@ class Field(object):
"""
return np.sqrt(abs(self.vdot(x=self)))
def squared_norm(self):
""" Computes the square of the L2-norm of the field values.
Returns
-------
float
The square of the L2-norm of the field values.
"""
return abs(self.vdot(x=self))
def conjugate(self):
""" Returns the complex conjugate of the field.
......
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 .wiener_filter_curvature import WienerFilterCurvature
from .wiener_filter_energy import WienerFilterEnergy
import numpy as np
from ..domains 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.
Parameters
----------
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,
SymmetrizingOperator)
from ..models import Variable, Constant, PointwiseExponential
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
smooth_spec = smooth_op(Variable(position)[keys[0]])
phi = Variable(position)[keys[1]] + Constant(position, norm_phi_mean)
linear_spec = slope(phi)
loglog_spec = smooth_spec + linear_spec
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.
Parameters
----------
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])
else:
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)
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from ..sugar import exp
from ..minimization.energy import Energy
from ..operators.diagonal_operator import DiagonalOperator
import numpy as np
class NoiseEnergy(Energy):
def __init__(self, position, alpha, q, res_sample_list):
super(NoiseEnergy, self).__init__(position)
self.N = DiagonalOperator(exp(position))
self.alpha = alpha
self.q = q
self.res_sample_list = res_sample_list
self._gradient = None
for s in res_sample_list:
lh = .5 * s.vdot(self.N.inverse_times(s))
grad = -.5 * self.N.inverse_times(s.conjugate()*s)
if self._gradient is None:
self._value = lh
self._gradient = grad
else:
self._value += lh
self._gradient += grad
expmpos = exp(-position)
self._value *= 1./len(res_sample_list)
possum = position.sum()
s1 = (alpha-1.)*possum if np.isscalar(alpha) \
else (alpha-1.).vdot(position)
s2 = q*expmpos.sum() if np.isscalar(q) else q.vdot(expmpos)
self._value += .5*possum + s1 + s2
self._gradient *= 1./len(res_sample_list)
self._gradient += (alpha-0.5) - q*expmpos
self._gradient.lock()
def at(self, position):
return self.__class__(position, self.alpha, self.q,
self.res_sample_list)
@property
def value(self):
return self._value
@property
def gradient(self):
return self._gradient
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from ..sugar import exp, makeOp
from ..minimization.energy import Energy
from ..operators.smoothness_operator import SmoothnessOperator
from ..operators.inversion_enabler import InversionEnabler
from ..utilities import memo
def _LinearizedPowerResponse(Instrument, nonlinearity, ht, Distributor, tau,
xi):
power = exp(0.5*tau)
position = ht(Distributor(power)*xi)
linearization = makeOp(nonlinearity.derivative(position))
return (0.5 * Instrument * linearization * ht * makeOp(xi) * Distributor *
makeOp(power))
class NonlinearPowerEnergy(Energy):
"""The Energy of the power spectrum according to the critical filter.
It describes the energy of the logarithmic amplitudes of the power spectrum
of a Gaussian random field under reconstruction uncertainty with smoothness
and inverse gamma prior. It is used to infer the correlation structure of a
correlated signal. The smoothness prior operates on logarithmic scale, i.e.
it prefers power-law-like power spectra.
Parameters
----------
position : Field
The current position of this energy.
xi : Field
The excitation field.
D : EndomorphicOperator
The curvature of the Gaussian encoding the posterior covariance.
If not specified, the map is assumed to be no reconstruction.
default : None
sigma : float
The parameter of the smoothness prior. Needs to be positive. A bigger
number means a stronger smoothness prior.
default : 0
samples : int
Number of samples used for the estimation of the uncertainty
corrections.
default : 3
"""
# MR FIXME: docstring incomplete and outdated
def __init__(self, position, d, N, xi, D, ht, Instrument, nonlinearity,
Distributor, sigma=0., samples=3, xi_sample_list=None,
iteration_controller=None):
super(NonlinearPowerEnergy, self).__init__(position)
self.xi = xi
self.D = D
self.d = d
self.N = N
self.T = SmoothnessOperator(domain=position.domain[0],
strength=sigma, logarithmic=True)
self.ht = ht
self.Instrument = Instrument
self.nonlinearity = nonlinearity
self.Distributor = Distributor
self.sigma = sigma
if xi_sample_list is None:
if samples is None or samples == 0:
xi_sample_list = [xi]
else:
xi_sample_list = [D.draw_sample(from_inverse=True) + xi
for _ in range(samples)]
self.xi_sample_list = xi_sample_list
self._ic = iteration_controller
A = Distributor(exp(.5 * position))
self._gradient = None
for xi_sample in self.xi_sample_list:
map_s = ht(A*xi_sample)
LinR = _LinearizedPowerResponse(Instrument, nonlinearity, ht,
Distributor, position, xi_sample)
residual = d - Instrument(nonlinearity(map_s))
tmp = N.inverse_times(residual)
lh = 0.5 * residual.vdot(tmp)
grad = LinR.adjoint_times(tmp)
if self._gradient is None:
self._value = lh
self._gradient = grad.copy()
else:
self._value += lh
self._gradient += grad
self._value *= 1. / len(self.xi_sample_list)
Tpos = self.T(position)
self._value += 0.5 * position.vdot(Tpos)
self._gradient *= -1. / len(self.xi_sample_list)
self._gradient += Tpos
self._gradient.lock()
def at(self, position):
return self.__class__(position, self.d, self.N, self.xi, self.D,
self.ht, self.Instrument, self.nonlinearity,
self.Distributor, sigma=self.sigma,
samples=len(self.xi_sample_list),
xi_sample_list=self.xi_sample_list,
iteration_controller=self._ic)
@property
def value(self):
return self._value
@property
def gradient(self):
return self._gradient
@property
@memo
def curvature(self):
result = None
for xi_sample in self.xi_sample_list:
LinearizedResponse = _LinearizedPowerResponse(
self.Instrument, self.nonlinearity, self.ht, self.Distributor,
self.position, xi_sample)
op = LinearizedResponse.adjoint*self.N.inverse*LinearizedResponse
result = op if result is None else result + op
result = result*(1./len(self.xi_sample_list)) + self.T
return InversionEnabler(result, self._ic)
......@@ -16,53 +16,13 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.