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


parents 41d128bb a1485b1c
......@@ -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 .kl import SampledKullbachLeiblerDivergence
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)
......@@ -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.
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 .nonlinear_wiener_filter_energy import NonlinearWienerFilterEnergy
from .nonlinearities import Exponential, Linear, PositiveTanh, Tanh
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
from .hamiltonian import Hamiltonian
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,
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.
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)
......@@ -45,7 +45,7 @@ class UnitLogGauss(Energy):
def value(self):
return .5 * self._s.value.vdot(self._s.value).real
return .5 * self._s.value.squared_norm()
......@@ -45,7 +45,7 @@ def WienerFilterCurvature(R, N, S, iteration_controller=None,
The iteration controller to use for sampling.
M = SandwichOperator.make(R, N.inverse)
if iteration_controller is not None:
if iteration_controller_sampling is not None:
op = SamplingEnabler(M, S.inverse, iteration_controller_sampling,
......@@ -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
......@@ -50,10 +50,7 @@ class Model(NiftyMetaBase()):
if isinstance(other, Model):
from .binary_helpers import Add
return Add.make(self, other)
if isinstance(other, (Field, MultiField)):
from .constant import Constant
return self.__add__(Constant(self.position, other))
raise TypeError
return NotImplemented
def __sub__(self, other):
return self.__add__((-1) * other)
......@@ -67,12 +64,12 @@ class Model(NiftyMetaBase()):
return Mul.make(self, other)
if isinstance(other, (Field, MultiField)):
return makeOp(other)(self)
raise NotImplementedError
return NotImplemented
def __rmul__(self, other):
if isinstance(other, (float, int, Field)):
return self.__mul__(other)
raise NotImplementedError
return NotImplemented
def __str__(self):
s = ('----------------------------------------'
......@@ -151,6 +151,16 @@ class MultiField(object):
return np.sqrt(np.abs(self.vdot(x=self)))
def squared_norm(self):
""" Computes the square of the L2-norm of the field values.
The square of the L2-norm of the field values.
return abs(self.vdot(x=self))
def __neg__(self):
return MultiField({key: -val for key, val in self.items()})
......@@ -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):
Supports Markdown
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