Commit d2fb3aee authored by Jakob Knollmueller's avatar Jakob Knollmueller
Browse files

merge master

parents f73bffdc 7585873e
......@@ -172,7 +172,7 @@
" tol_abs_gradnorm=0.1)\n",
" # WienerFilterCurvature is (R.adjoint*N.inverse*R + Sh.inverse) plus some handy\n",
" # helper methods.\n",
" return ift.library.WienerFilterCurvature(R,N,Sh,iteration_controller=IC,iteration_controller_sampling=IC)"
" return ift.WienerFilterCurvature(R,N,Sh,iteration_controller=IC,iteration_controller_sampling=IC)"
]
},
{
......
......@@ -34,7 +34,7 @@ if __name__ == '__main__':
# position_space = ift.HPSpace(128)
# mask = make_random_mask()
# set up corresponding harmonic transform and space
harmonic_space = position_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space)
......
......@@ -16,7 +16,7 @@ from .minimization import *
from .sugar import *
from .plotting.plot import plot
from . import library
from .library import *
from . import extra
from .utilities import memo
......
......@@ -34,7 +34,7 @@ class GLSpace(StructuredDomain):
pixelization.
nlon : int, optional
Number of longitudinal bins that are used for this pixelization.
Default value is 2*nlat + 1.
Default value is 2*nlat - 1.
"""
_needed_for_hash = ["_nlat", "_nlon"]
......
# 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 ..library.gaussian_energy import GaussianEnergy
from ..minimization.energy import Energy
from ..operators import InversionEnabler, SamplingEnabler
from ..models.variable import Variable
from ..operators import InversionEnabler, SamplingEnabler
from ..utilities import memo
from ..library.unit_log_gauss import UnitLogGauss
class Hamiltonian(Energy):
......@@ -15,11 +33,8 @@ class Hamiltonian(Energy):
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._ic_samp = iteration_controller_sampling
self._prior = GaussianEnergy(Variable(self.position))
self._precond = self._prior.curvature
def at(self, position):
......@@ -39,8 +54,11 @@ class Hamiltonian(Energy):
@memo
def curvature(self):
prior_curv = self._prior.curvature
c = SamplingEnabler(self._lh.curvature, prior_curv.inverse,
self._ic_samp, prior_curv.inverse)
if self._ic_samp is None:
c = self._lh.curvature + prior_curv
else:
c = SamplingEnabler(self._lh.curvature, prior_curv.inverse,
self._ic_samp, prior_curv.inverse)
return InversionEnabler(c, self._ic, self._precond)
def __str__(self):
......
from .amplitude_model import make_amplitude_model
from .apply_data import ApplyData
from .gaussian_energy import GaussianEnergy
from .los_response import LOSResponse
from .nonlinear_wiener_filter_energy import NonlinearWienerFilterEnergy
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_correlated_field
from .wiener_filter_curvature import WienerFilterCurvature
from .wiener_filter_energy import WienerFilterEnergy
def ApplyData(data, var, model_data):
from .. import DiagonalOperator, Constant, sqrt
# TODO This is rather confusing. Delete that eventually.
from ..operators.diagonal_operator import DiagonalOperator
from ..models.constant import Constant
from ..sugar import sqrt
sqrt_n = DiagonalOperator(sqrt(var))
data = Constant(model_data.position, data)
return sqrt_n.inverse(model_data - data)
......@@ -17,45 +17,50 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from ..minimization.energy import Energy
from ..operators.inversion_enabler import InversionEnabler
from ..operators.sandwich_operator import SandwichOperator
from ..utilities import memo
class UnitLogGauss(Energy):
def __init__(self, s, inverter=None):
class GaussianEnergy(Energy):
def __init__(self, inp, mean=None, covariance=None):
"""
s: Sky model object
inp: Model object
value = 0.5 * s.vdot(s), i.e. a log-Gauss distribution with unit
covariance
"""
super(UnitLogGauss, self).__init__(s.position)
self._s = s
self._inverter = inverter
super(GaussianEnergy, self).__init__(inp.position)
self._inp = inp
self._mean = mean
self._cov = covariance
def at(self, position):
return self.__class__(self._s.at(position), self._inverter)
return self.__class__(self._inp.at(position), self._mean, self._cov)
@property
@memo
def _gradient_helper(self):
return self._s.gradient
def residual(self):
if self._mean is not None:
return self._inp.value - self._mean
return self._inp.value
@property
@memo
def value(self):
return .5 * self._s.value.squared_norm()
if self._cov is None:
return .5 * self.residual.vdot(self.residual).real
return .5 * self.residual.vdot(self._cov.inverse(self.residual)).real
@property
@memo
def gradient(self):
return self._gradient_helper.adjoint(self._s.value)
if self._cov is None:
return self._inp.gradient.adjoint(self.residual)
return self._inp.gradient.adjoint(self._cov.inverse(self.residual))
@property
@memo
def curvature(self):
c = SandwichOperator.make(self._gradient_helper)
if self._inverter is None:
return c
return InversionEnabler(c, self._inverter)
if self._cov is None:
return SandwichOperator.make(self._inp.gradient, None)
return SandwichOperator.make(self._inp.gradient, self._cov.inverse)
......@@ -23,7 +23,7 @@ from ..operators.sandwich_operator import SandwichOperator
from ..sugar import log, makeOp
class PoissonLogLikelihood(Energy):
class PoissonianEnergy(Energy):
def __init__(self, lamb, d):
"""
lamb: Sky model object
......@@ -31,7 +31,7 @@ class PoissonLogLikelihood(Energy):
value = 0.5 * s.vdot(s), i.e. a log-Gauss distribution with unit
covariance
"""
super(PoissonLogLikelihood, self).__init__(lamb.position)
super(PoissonianEnergy, self).__init__(lamb.position)
self._lamb = lamb
self._d = d
......
......@@ -164,10 +164,25 @@ class MultiField(object):
def __neg__(self):
return MultiField({key: -val for key, val in self.items()})
def __abs__(self):
return MultiField({key: abs(val) for key, val in self.items()})
def conjugate(self):
return MultiField({key: sub_field.conjugate()
for key, sub_field in self.items()})
def all(self):
for v in self.values():
if not v.all():
return False
return True
def any(self):
for v in self.values():
if v.any():
return True
return False
def isEquivalentTo(self, other):
"""Determines (as quickly as possible) whether `self`'s content is
identical to `other`'s content."""
......
......@@ -10,6 +10,7 @@ from .harmonic_transform_operator import HarmonicTransformOperator
from .inversion_enabler import InversionEnabler
from .laplace_operator import LaplaceOperator
from .linear_operator import LinearOperator
from .mask_operator import MaskOperator
from .multi_adaptor import MultiAdaptor
from .power_distributor import PowerDistributor
from .qht_operator import QHTOperator
......@@ -23,7 +24,7 @@ from .symmetrizing_operator import SymmetrizingOperator
__all__ = ["LinearOperator", "EndomorphicOperator", "ScalingOperator",
"DiagonalOperator", "HarmonicTransformOperator", "FFTOperator",
"FFTSmoothingOperator", "GeometryRemover",
"FFTSmoothingOperator", "GeometryRemover", "MaskOperator",
"LaplaceOperator", "SmoothnessOperator", "PowerDistributor",
"InversionEnabler", "SandwichOperator", "SamplingEnabler",
"DOFDistributor", "SelectionOperator", "MultiAdaptor",
......
......@@ -16,13 +16,48 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from ..models.constant import Constant
from .unit_log_gauss import UnitLogGauss
from ..energies.hamiltonian import Hamiltonian
import numpy as np
from ..domain_tuple import DomainTuple
from ..domains.unstructured_domain import UnstructuredDomain
from ..field import Field
from .linear_operator import LinearOperator
def NonlinearWienerFilterEnergy(measured_data, data_model, sqrtN, iteration_controller):
d = measured_data.lock()
residual = Constant(data_model.position, d) - data_model
lh = UnitLogGauss(sqrtN.inverse(residual))
return Hamiltonian(lh, iteration_controller)
class MaskOperator(LinearOperator):
def __init__(self, mask):
if not isinstance(mask, Field):
raise TypeError
self._domain = DomainTuple.make(mask.domain)
self._mask = np.logical_not(mask.to_global_data())
self._target = DomainTuple.make(UnstructuredDomain(self._mask.sum()))
def data_indices(self):
if len(self.domain.shape) == 1:
return np.arange(self.domain.shape[0])[self._mask]
if len(self.domain.shape) == 2:
return np.indices(self.domain.shape).transpose((1, 2, 0))[self._mask]
def apply(self, x, mode):
self._check_input(x, mode)
if mode == self.TIMES:
res = x.to_global_data()[self._mask]
return Field.from_global_data(self.target, res)
x = x.to_global_data()
res = np.empty(self.domain.shape, x.dtype)
res[self._mask] = x
res[~self._mask] = 0
return Field.from_global_data(self.domain, res)
@property
def capability(self):
return self.TIMES | self.ADJOINT_TIMES
@property
def domain(self):
return self._domain
@property
def target(self):
return self._target
......@@ -24,16 +24,17 @@ from .endomorphic_operator import EndomorphicOperator
class SamplingEnabler(EndomorphicOperator):
"""Class which augments the capability of another operator object via
numerical inversion.
"""Class which acts as the operator object (likelihood + prior)
and enables sampling from its inverse even if the operator object
itself does not support it.
Parameters
----------
op : :class:`EndomorphicOperator`
The operator to be enhanced.
The InversionEnabler object will support the same operation modes as
`op`, and additionally the inverse set. The newly-added modes will
be computed by iterative inversion.
likelihood : :class:`EndomorphicOperator`
Curvature of the likelihood
prior : :class:`EndomorphicOperator`
Inverse curvature of the prior
iteration_controller : :class:`IterationController`
The iteration controller to use for the iterative numerical inversion
done by a :class:`ConjugateGradient` object.
......
......@@ -57,7 +57,7 @@ class Energy_Tests(unittest.TestCase):
tol_abs_gradnorm=1e-5)
S = ift.create_power_operator(hspace, power_spectrum=_flat_PS)
energy = ift.library.WienerFilterEnergy(
energy = ift.WienerFilterEnergy(
position=s0, d=d, R=R, N=N, S=S, iteration_controller=IC)
ift.extra.check_value_gradient_curvature_consistency(
energy, ntries=10)
......@@ -66,10 +66,10 @@ class Energy_Tests(unittest.TestCase):
ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[ift.Tanh, ift.Exponential, ift.Linear],
[1, 1e-2, 1e2],
[4, 78, 23]))
def testNonlinearMap(self, space, nonlinearity, seed):
def testGaussianEnergy(self, space, nonlinearity, noise, seed):
np.random.seed(seed)
f = nonlinearity()
dim = len(space.shape)
hspace = space.get_default_codomain()
ht = ift.HarmonicTransformOperator(hspace, target=space)
......@@ -77,23 +77,23 @@ class Energy_Tests(unittest.TestCase):
pspace = ift.PowerSpace(hspace, binbounds=binbounds)
Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
xi0_var = ift.Variable(ift.MultiField({'xi':xi0}))['xi']
xi0_var = ift.Variable(ift.MultiField({'xi': xi0}))['xi']
def pspec(k): return 1 / (1 + k**2)**dim
pspec = ift.PS_field(pspace, pspec)
A = Dist(ift.sqrt(pspec))
n = ift.Field.from_random(domain=space, random_type='normal')
N = ift.ScalingOperator(noise, space)
n = N.draw_sample()
s = ht(ift.makeOp(A)(xi0_var))
R = ift.ScalingOperator(10., space)
sqrtN = ift.ScalingOperator(1., space)
d_model = R(ift.LocalModel(s, nonlinearity()))
d = d_model.value + n
IC = ift.GradientNormController(iteration_limit=100,
tol_abs_gradnorm=1e-5)
energy = ift.library.NonlinearWienerFilterEnergy(
d, d_model, sqrtN, IC)
if isinstance(nonlinearity, ift.Linear):
if noise == 1:
N = None
energy = ift.GaussianEnergy(d_model, d, N)
if isinstance(nonlinearity(), ift.Linear):
ift.extra.check_value_gradient_curvature_consistency(
energy, ntries=10)
else:
......
......@@ -61,6 +61,17 @@ class Consistency_Tests(unittest.TestCase):
op = ift.HarmonicTransformOperator(sp)
ift.extra.consistency_check(op, dtype, dtype)
@expand(product(_p_spaces, [np.float64, np.complex128]))
def testMask(self, sp, dtype):
# Create mask
f = ift.from_random('normal', sp).to_global_data()
mask = np.zeros_like(f)
mask[f > 0] = 1
mask = ift.Field.from_global_data(sp, mask)
# Test MaskOperator
op = ift.MaskOperator(mask)
ift.extra.consistency_check(op, dtype, dtype)
@expand(product(_h_spaces+_p_spaces, [np.float64, np.complex128]))
def testDiagonal(self, sp, dtype):
op = ift.DiagonalOperator(ift.Field.from_random("normal", sp,
......
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