Skip to content
GitLab
Commit b6c358a0 authored by Philipp Arras's avatar Philipp Arras
Browse files

Merge branch 'new_energies' into 'NIFTy_5' 

New energies

See merge request ift/nifty-dev!2
parents a1485b1c e4337b58
Hide whitespace changes
Inline Side-by-side
nifty5/energies/__init__.py
View file @ b6c358a0
from .hamiltonian import Hamiltonian
from .kl import SampledKullbachLeiblerDivergence
from .hamiltonian import Hamiltonian
nifty5/energies/hamiltonian.py
View file @ b6c358a0
from nifty5 import Energy, InversionEnabler, SamplingEnabler, Variable, memo
from nifty5.library import UnitLogGauss
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):
......
nifty5/energies/kl.py
View file @ b6c358a0
from nifty5 import Energy, InversionEnabler, ScalingOperator, memo
from ..minimization.energy import Energy
from ..operators.inversion_enabler import InversionEnabler
from ..operators.scaling_operator import ScalingOperator
from ..utilities import memo
class SampledKullbachLeiblerDivergence(Energy):
......
nifty5/library/__init__.py
View file @ b6c358a0
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 .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
nifty5/library/noise_energy.py deleted 100644 → 0
View file @ a1485b1c
# 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
nifty5/library/nonlinear_power_energy.py deleted 100644 → 0
View file @ a1485b1c
# 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)
nifty5/library/nonlinear_wiener_filter_energy.py
View file @ b6c358a0
......@@ -16,53 +16,13 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from .wiener_filter_curvature import WienerFilterCurvature
from ..utilities import memo
from ..minimization.energy import Energy
from ..sugar import makeOp
from ..models.constant import Constant
from .unit_log_gauss import UnitLogGauss
from ..energies.hamiltonian import Hamiltonian
class NonlinearWienerFilterEnergy(Energy):
def __init__(self, position, d, Instrument, nonlinearity, ht, power, N, S,
iteration_controller=None,
iteration_controller_sampling=None):
super(NonlinearWienerFilterEnergy, self).__init__(position=position)
self.d = d.lock()
self.Instrument = Instrument
self.nonlinearity = nonlinearity
self.ht = ht
self.power = power
m = ht(power*position)
residual = d - Instrument(nonlinearity(m))
self.N = N
self.S = S
self._ic = iteration_controller
if iteration_controller_sampling is None:
iteration_controller_sampling = self._ic
self._ic_samp = iteration_controller_sampling
t1 = S.inverse_times(position)
t2 = N.inverse_times(residual)
self._value = 0.5 * (position.vdot(t1) + residual.vdot(t2)).real
self.R = (Instrument * makeOp(nonlinearity.derivative(m)) *
ht * makeOp(power))
self._gradient = (t1 - self.R.adjoint_times(t2)).lock()
def at(self, position):
return self.__class__(position, self.d, self.Instrument,
self.nonlinearity, self.ht, self.power, self.N,
self.S, self._ic, self._ic_samp)
@property
def value(self):
return self._value
@property
def gradient(self):
return self._gradient
@property
@memo
def curvature(self):
return WienerFilterCurvature(self.R, self.N, self.S, self._ic,
self._ic_samp)
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)
nifty5/library/poisson_energy.py deleted 100644 → 0
View file @ a1485b1c
# 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 ..minimization.energy import Energy
from ..operators.diagonal_operator import DiagonalOperator
from ..operators.sandwich_operator import SandwichOperator
from ..operators.inversion_enabler import InversionEnabler
from ..sugar import log
class PoissonEnergy(Energy):
def __init__(self, position, d, Instrument, nonlinearity, ht, Phi_h,
iteration_controller=None):
super(PoissonEnergy, self).__init__(position=position)
self._ic = iteration_controller
self._d = d
self._Instrument = Instrument
self._nonlinearity = nonlinearity
self._ht = ht
self._Phi_h = Phi_h
htpos = ht(position)
lam = Instrument(nonlinearity(htpos))
Rho = DiagonalOperator(nonlinearity.derivative(htpos))
eps = 1e-100 # to catch harmless 0/0 where data is zero
W = DiagonalOperator((d+eps)/(lam**2+eps))
phipos = Phi_h.inverse_adjoint_times(position)
prior_energy = 0.5*position.vdot(phipos)
likel_energy = lam.sum()-d.vdot(log(lam+eps))
self._value = prior_energy + likel_energy
R1 = Instrument*Rho*ht
self._grad = (phipos + R1.adjoint_times((lam-d)/(lam+eps))).lock()
self._curv = Phi_h.inverse + SandwichOperator.make(R1, W)
def at(self, position):
return self.__class__(position, self._d, self._Instrument,
self._nonlinearity, self._ht, self._Phi_h,
self._ic)
@property
def value(self):
return self._value
@property
def gradient(self):
return self._grad
@property
def curvature(self):
return InversionEnabler(self._curv, self._ic, self._Phi_h.inverse)
nifty5/library/poisson_log_likelihood.py
View file @ b6c358a0
......@@ -26,7 +26,7 @@ from ..sugar import log, makeOp
class PoissonLogLikelihood(Energy):
def __init__(self, lamb, d):
"""
s: Sky model object
lamb: Sky model object
value = 0.5 * s.vdot(s), i.e. a log-Gauss distribution with unit
covariance
......
test/test_energies/test_map.py
View file @ b6c358a0
......@@ -77,20 +77,22 @@ 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']
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')
s = ht(xi0 * A)
s = ht(ift.makeOp(A)(xi0_var))
R = ift.ScalingOperator(10., space)
N = ift.ScalingOperator(1., space)
d = R(f(s)) + n
sqrtN = ift.ScalingOperator(1., space)
d_model = R(ift.LocalModel(s, nonlinearity()))
d = d_model.value + n
S = ift.create_power_operator(hspace, power_spectrum=_flat_PS)
IC = ift.GradientNormController(iteration_limit=100,
tol_abs_gradnorm=1e-5)
energy = ift.library.NonlinearWienerFilterEnergy(
position=xi0, d=d, Instrument=R, nonlinearity=f, ht=ht, power=A,
N=N, S=S)
d, d_model, sqrtN, IC)
if isinstance(nonlinearity, ift.Linear):
ift.extra.check_value_gradient_curvature_consistency(
energy, ntries=10)
......
test/test_energies/test_noise.py deleted 100644 → 0
View file @ a1485b1c
# 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.
import unittest
import nifty5 as ift
import numpy as np
from itertools import product
from test.common import expand
from numpy.testing import assert_allclose
def _flat_PS(k):
return np.ones_like(k)
class Noise_Energy_Tests(unittest.TestCase):
@expand(product([ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[ift.Exponential, ift.Linear],
[23, 131, 32]))
def testNoise(self, space, nonlinearity, seed):
np.random.seed(seed)
f = nonlinearity()
dim = len(space.shape)
hspace = space.get_default_codomain()
ht = ift.HarmonicTransformOperator(hspace, target=space)
binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
pspace = ift.PowerSpace(hspace, binbounds=binbounds)
Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
xi = ift.Field.from_random(domain=hspace, random_type='normal')
def pspec(k): return 1 / (1 + k**2)**dim
tau = ift.PS_field(pspace, pspec)
A = Dist(ift.sqrt(tau))
var = 1.
n = ift.Field.from_random(
domain=space,
random_type='normal',
std=np.sqrt(var))
var = ift.Field.full(n.domain, var)
N = ift.DiagonalOperator(var)
eta0 = ift.log(var)
s = ht(xi * A)
R = ift.ScalingOperator(10., space)
d = R(f(s)) + n
alpha = ift.Field.full(d.domain, 2.)
q = ift.Field.full(d.domain, 1e-5)
IC = ift.GradientNormController(
iteration_limit=100,
tol_abs_gradnorm=1e-5)
S = ift.create_power_operator(hspace, power_spectrum=_flat_PS)
C = ift.library.NonlinearWienerFilterEnergy(
position=xi,
d=d,
Instrument=R,
nonlinearity=f,
ht=ht,
power=A,
N=N,
S=S,
iteration_controller=IC).curvature
res_sample_list = [d - R(f(ht(C.draw_sample(from_inverse=True) + xi)))
for _ in range(10)]
energy = ift.library.NoiseEnergy(eta0, alpha, q, res_sample_list)
ift.extra.check_value_gradient_consistency(energy, ntries=10)
test/test_energies/test_power.py deleted 100644 → 0
View file @ a1485b1c
# 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.
import unittest
import nifty5 as ift
import numpy as np
from itertools import product
from test.common import expand
def _flat_PS(k):
return np.ones_like(k)
class Energy_Tests(unittest.TestCase):
@expand(product([ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[ift.Exponential, ift.Linear],
[132, 42, 3]))
def testNonlinearPower(self, space, nonlinearity, seed):
np.random.seed(seed)
f = nonlinearity()
dim = len(space.shape)
hspace = space.get_default_codomain()
ht = ift.HarmonicTransformOperator(hspace, space)
binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=True)
pspace = ift.PowerSpace(hspace, binbounds=binbounds)
Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
xi = ift.Field.from_random(domain=hspace, random_type='normal')
def pspec(k): return 1 / (1 + k**2)**dim
tau0 = ift.PS_field(pspace, pspec)
A = Dist(ift.sqrt(tau0))
n = ift.Field.from_random(domain=space, random_type='normal')
s = ht(xi * A)
R = ift.ScalingOperator(10., space)
diag = ift.full(space, 1.)
N = ift.DiagonalOperator(diag)
d = R(f(s)) + n
IC = ift.GradientNormController(
iteration_limit=100,
tol_abs_gradnorm=1e-5)
S = ift.create_power_operator(hspace, power_spectrum=_flat_PS)
D = ift.library.NonlinearWienerFilterEnergy(
position=xi,
d=d,
Instrument=R,
nonlinearity=f,
power=A,
N=N,
S=S,
ht=ht,
iteration_controller=IC).curvature
energy = ift.library.NonlinearPowerEnergy(
position=tau0,
d=d,
xi=xi,
D=D,
Instrument=R,
Distributor=Dist,
nonlinearity=f,
ht=ht,
N=N,
samples=10)
ift.extra.check_value_gradient_consistency(energy, ntries=10)