Commit ab072f88 authored by Philipp Arras's avatar Philipp Arras
Browse files

Merge branch 'energy_tests' into addUnits

parents 8d604d2c 5dbe3b0d
......@@ -6,7 +6,7 @@ import numericalunits as nu
if __name__ == "__main__":
# In MPI mode, the random seed for numericalunits must be set by hand
#nu.reset_units(43)
dimensionality = 1
dimensionality = 2
np.random.seed(43)
# Setting up variable parameters
......@@ -48,14 +48,14 @@ if __name__ == "__main__":
np.random.seed(43)
mock_power = ift.PS_field(power_space, power_spectrum)
mock_harmonic = ift.power_synthesize(mock_power, real_signal=True)
mock_signal = fft(mock_harmonic)
mock_signal = ift.power_synthesize(mock_power, real_signal=True)
sensitivity = (1./nu.m)**dimensionality/nu.K
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
sensitivity=(sensitivity,))
R = ift.GeometryRemover(signal_space)
R = R*ift.ScalingOperator(sensitivity, signal_space)
R = R*fft
R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma)
data_domain = R.target[0]
R_harmonic = R*fft
noise_amplitude = 1./signal_to_noise*field_sigma*sensitivity*((L/N_pixels)**dimensionality)
print "noise amplitude:", noise_amplitude
......@@ -67,22 +67,22 @@ if __name__ == "__main__":
data = R(mock_signal) + noise
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data))
j = R.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(
verbose=True, tol_abs_gradnorm=1e-5/(nu.K*(nu.m**dimensionality)))
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R_harmonic, inverter=inverter)
S=S, N=N, R=R, inverter=inverter)
m = wiener_curvature.inverse_times(j)
m_s = fft(m)
sspace2 = ift.RGSpace(shape, distances=L/N_pixels/nu.m)
ift.plot(ift.Field(sspace2, mock_signal.val)/nu.K, title="mock_signal.png")
ift.plot(ift.Field(sspace2, fft(mock_signal).val)/nu.K, name="mock_signal.png")
#data = ift.dobj.to_global_data(data.val).reshape(sspace2.shape)
#data = ift.Field(sspace2, val=ift.dobj.from_global_data(data))
ift.plot(ift.Field(sspace2, val=R.adjoint_times(data).val), title="data.png")
print "msig",np.min(mock_signal.val)/nu.K, np.max(mock_signal.val)/nu.K
ift.plot(ift.Field(sspace2, val=fft(R.adjoint_times(data)).val), name="data.png")
print "msig",np.min(fft(mock_signal).val)/nu.K, np.max(fft(mock_signal).val)/nu.K
print "map",np.min(m_s.val)/nu.K, np.max(m_s.val)/nu.K
ift.plot(ift.Field(sspace2, m_s.val)/nu.K, title="map.png")
ift.plot(ift.Field(sspace2, m_s.val)/nu.K, name="map.png")
......@@ -24,6 +24,7 @@ from .operators.fft_operator import FFTOperator
from .operators.fft_smoothing_operator import FFTSmoothingOperator
from .operators.direct_smoothing_operator import DirectSmoothingOperator
from .operators.response_operator import ResponseOperator
from .operators.response_operator import GeometryRemover
from .operators.laplace_operator import LaplaceOperator
from .operators.power_projection_operator import PowerProjectionOperator
from .operators.inversion_enabler import InversionEnabler
......
......@@ -35,7 +35,7 @@ class CriticalPowerEnergy(Energy):
Parameters
----------
position : Field,
The current position of this energy.
The current position of this energy. (Logarithm of power spectrum)
m : Field,
The map whose power spectrum has to be inferred
D : EndomorphicOperator,
......@@ -101,8 +101,8 @@ class CriticalPowerEnergy(Energy):
self._theta = exp(-self.position) * (self.q + self._w*0.5)
Tt = self.T(self.position)
energy = self._theta.integrate()
energy += self.position.integrate()*(self.alpha-0.5)
energy = self._theta.vdot(Field.ones_like(self._theta))
energy += self.position.vdot(Field.ones_like(self.position)) *(self.alpha-0.5)
energy += 0.5*self.position.vdot(Tt)
self._value = energy.real
......
......@@ -43,7 +43,7 @@ class LogNormalWienerFilterEnergy(Energy):
The prior signal covariance in harmonic space.
"""
def __init__(self, position, d, R, N, S, inverter, fft=None):
def __init__(self, position, d, R, N, S, inverter, ht=None):
super(LogNormalWienerFilterEnergy, self).__init__(position=position)
self.d = d
self.R = R
......@@ -51,21 +51,21 @@ class LogNormalWienerFilterEnergy(Energy):
self.S = S
self._inverter = inverter
if fft is None:
self._fft = create_composed_fft_operator(self.S.domain,
all_to='position')
if ht is None:
self._ht = create_composed_fft_operator(self.S.domain,
all_to='position')
else:
self._fft = fft
self._ht = ht
self._expp_sspace = exp(self._fft(self.position))
self._expp_sspace = exp(self._ht(self.position))
Sp = self.S.inverse_times(self.position)
expp = self._fft.adjoint_times(self._expp_sspace)
expp = self._ht.adjoint_times(self._expp_sspace)
Rexppd = self.R(expp) - self.d
NRexppd = self.N.inverse_times(Rexppd)
self._value = 0.5*(self.position.vdot(Sp) + Rexppd.vdot(NRexppd))
exppRNRexppd = self._fft.adjoint_times(
self._expp_sspace * self._fft(self.R.adjoint_times(NRexppd)))
exppRNRexppd = self._ht.adjoint_times(
self._expp_sspace * self._ht(self.R.adjoint_times(NRexppd)))
self._gradient = Sp + exppRNRexppd
def at(self, position):
......
......@@ -20,9 +20,11 @@ from .. import Field, exp
from ..operators.diagonal_operator import DiagonalOperator
from ..minimization.energy import Energy
# TODO Take only residual_sample_list as argument
class NoiseEnergy(Energy):
def __init__(self, position, d, m, D, t, HarmonicTransform, Instrument,
def __init__(self, position, d, m, D, t, ht, Instrument,
nonlinearity, alpha, q, Projection, munit=1., sunit=1.,
dunit=1., samples=3, sample_list=None, inverter=None):
super(NoiseEnergy, self).__init__(position=position)
......@@ -32,7 +34,7 @@ class NoiseEnergy(Energy):
self.N = DiagonalOperator(diagonal=dunit**2 * exp(self.position))
self.t = t
self.samples = samples
self.ht = HarmonicTransform
self.ht = ht
self.Instrument = Instrument
self.nonlinearity = nonlinearity
self.munit = munit
......@@ -73,7 +75,7 @@ class NoiseEnergy(Energy):
self._gradient += grad
self._value *= 1. / len(self.sample_list)
self._value += .5 * self.position.integrate()
self._value += .5 * self.position.sum()
self._value += (self.alpha - 1.).vdot(self.position) + \
self.q.vdot(exp(-self.position))
......
......@@ -20,15 +20,13 @@ from ..operators.inversion_enabler import InversionEnabler
from .response_operators import LinearizedPowerResponse
def NonlinearPowerCurvature(position, HarmonicTransform, Instrument,
nonlinearity, Projection, N, T, sample_list,
inverter, munit=1., sunit=1.):
def NonlinearPowerCurvature(position, ht, Instrument, nonlinearity,
Projection, N, T, sample_list, inverter, munit=1., sunit=1.):
result = None
for sample in sample_list:
LinR = LinearizedPowerResponse(Instrument, nonlinearity,
HarmonicTransform, Projection, position,
sample, munit, sunit)
op = LinR.adjoint * N.inverse * LinR
LinearizedResponse = LinearizedPowerResponse(
Instrument, nonlinearity, ht, Projection, position, sample, munit, sunit)
op = LinearizedResponse.adjoint*N.inverse*LinearizedResponse
result = op if result is None else result + op
result = result * (1. / len(sample_list)) + T
return InversionEnabler(result, inverter)
......@@ -51,11 +51,17 @@ class NonlinearPowerEnergy(Energy):
default : 3
"""
def __init__(self, position, d, N, m, D, HarmonicTransform, Instrument,
nonlinearity, Projection, sigma=0., samples=3,
sample_list=None, munit=1., sunit=1., inverter=None):
def __init__(self, position, d, N, m, D, ht, Instrument, nonlinearity,
Projection, sigma=0., samples=3, sample_list=None,
inverter=None, munit=1., sunit=1.):
super(NonlinearPowerEnergy, self).__init__(position)
self.d, self.N, self.m, self.D, self.ht = d, N, m, D, HarmonicTransform
self.m = m
self.D = D
self.d = d
self.N = N
self.T = SmoothnessOperator(domain=self.position.domain[0],
strength=sigma, logarithmic=True)
self.ht = ht
self.Instrument = Instrument
self.nonlinearity = nonlinearity
self.Projection = Projection
......@@ -71,9 +77,6 @@ class NonlinearPowerEnergy(Energy):
self.sample_list = sample_list
self.inverter = inverter
self.T = SmoothnessOperator(domain=self.position.domain[0],
strength=sigma, logarithmic=True)
A = Projection.adjoint_times(munit * exp(.5 * position)) # unit: munit
map_s = self.ht(A * m)
Tpos = self.T(position)
......@@ -92,8 +95,12 @@ class NonlinearPowerEnergy(Energy):
sunit)
residual = self.d - \
self.Instrument(sunit * self.nonlinearity(map_s))
self.Instrument(sunit * self.nonlinearity(
self.ht(self.power*sample)))
lh = 0.5 * residual.vdot(self.N.inverse_times(residual))
LinR = LinearizedPowerResponse(
self.Instrument, self.nonlinearity, self.ht, self.Projection,
self.position, sample)
grad = LinR.adjoint_times(self.N.inverse_times(residual))
if self._gradient is None:
......
......@@ -23,19 +23,20 @@ from .response_operators import LinearizedSignalResponse
class NonlinearWienerFilterEnergy(Energy):
def __init__(self, position, d, Instrument, nonlinearity, HarmonicTransform,
power, N, S, sunit=1., inverter=None):
def __init__(self, position, d, Instrument, nonlinearity, ht, power, N, S,
inverter=None, sunit=1.):
super(NonlinearWienerFilterEnergy, self).__init__(position=position)
self.d = d
self.sunit = sunit
self.Instrument = Instrument
self.nonlinearity = nonlinearity
self.ht = HarmonicTransform
self.ht = ht
self.power = power
position_map = self.ht(self.power * self.position)
self.LinearizedResponse = \
LinearizedSignalResponse(Instrument, nonlinearity, self.ht, power,
position_map, sunit)
self.LinearizedResponse = LinearizedSignalResponse(Instrument, nonlinearity, ht, power,
position_map, sunit)
position_map = ht(self.power * self.position)
residual = d - Instrument(sunit * nonlinearity(position_map))
self.N = N
self.S = S
......
......@@ -19,16 +19,13 @@
from ..field import exp
def LinearizedSignalResponse(Instrument, nonlinearity, HarmonicTransform, power,
s, sunit):
return sunit * (Instrument * nonlinearity.derivative(s) *
HarmonicTransform * power)
def LinearizedSignalResponse(Instrument, nonlinearity, ht, power, s, sunit):
return sunit * (Instrument * nonlinearity.derivative(s) * ht * power)
def LinearizedPowerResponse(Instrument, nonlinearity, HarmonicTransform,
Projection, t, m, munit, sunit):
power = exp(0.5 * t) * munit
position = HarmonicTransform(Projection.adjoint_times(power) * m)
def LinearizedPowerResponse(Instrument, nonlinearity, ht, Projection, t, m, munit, sunit):
power = exp(0.5*t) * munit
position = ht(Projection.adjoint_times(power) * m)
linearization = nonlinearity.derivative(position)
return sunit * (0.5 * Instrument * linearization * HarmonicTransform * m *
return sunit * (0.5 * Instrument * linearization * ht * m *
Projection.adjoint * power)
......@@ -32,7 +32,8 @@ __all__ = ['PS_field',
'power_synthesize_nonrandom',
'create_power_field',
'create_power_operator',
'create_composed_fft_operator']
'create_composed_fft_operator',
'create_harmonic_smoothing_operator']
def PS_field(pspace, func, dtype=None):
......@@ -256,3 +257,9 @@ def create_composed_fft_operator(domain, codomain=None, all_to='other'):
if res is None:
raise ValueError("empty operator")
return res
def create_harmonic_smoothing_operator(domain, space, sigma):
kfunc = domain[space].get_fft_smoothing_kernel_function(sigma)
return DiagonalOperator(kfunc(domain[space].get_k_length_array()), domain,
space)
# 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-2017 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 nifty4 as ift
import numpy as np
from itertools import product
from test.common import expand
from numpy.testing import assert_allclose
# TODO Add also other space types
# TODO Set tolerances and eps to reasonable values
class Map_Energy_Tests(unittest.TestCase):
@expand(product([ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[4, 78, 23]))
def testLinearMap(self, space, seed):
np.random.seed(seed)
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)
P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
def pspec(k): return 1 / (1 + k**2)**dim
pspec = ift.PS_field(pspace, pspec)
A = P.adjoint_times(ift.sqrt(pspec))
n = ift.Field.from_random(domain=space, random_type='normal')
s0 = xi0 * A
diag = ift.Field.ones(space) * 10
Instrument = ift.DiagonalOperator(diag)
R = Instrument * ht
diag = ift.Field.ones(space)
N = ift.DiagonalOperator(diag)
d = R(s0) + n
direction = ift.Field.from_random('normal', hspace)
direction /= np.sqrt(direction.var())
eps = 1e-7
s1 = s0 + eps * direction
IC = ift.GradientNormController(
name='IC',
verbose=False,
iteration_limit=100,
tol_abs_gradnorm=1e-5)
inverter = ift.ConjugateGradient(IC)
S = ift.create_power_operator(hspace, power_spectrum=lambda k: 1.)
energy0 = ift.library.WienerFilterEnergy(
position=s0, d=d, R=R, N=N, S=S, inverter=inverter)
energy1 = ift.library.WienerFilterEnergy(
position=s1, d=d, R=R, N=N, S=S, inverter=inverter)
a = (energy1.value - energy0.value) / eps
b = energy0.gradient.vdot(direction)
tol = 1e-5
assert_allclose(a, b, rtol=tol, atol=tol)
@expand(product([ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[4, 78, 23]))
def testLognormalMap(self, space, seed):
np.random.seed(seed)
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)
P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
def pspec(k): return 1 / (1 + k**2)**dim
pspec = ift.PS_field(pspace, pspec)
A = P.adjoint_times(ift.sqrt(pspec))
n = ift.Field.from_random(domain=space, random_type='normal')
sh0 = xi0 * A
s = ht(sh0)
diag = ift.Field.ones(space) * 10
Instrument = ift.DiagonalOperator(diag)
R = Instrument * ht
diag = ift.Field.ones(space)
N = ift.DiagonalOperator(diag)
d = Instrument(ift.exp(s)) + n
direction = ift.Field.from_random('normal', hspace)
direction /= np.sqrt(direction.var())
eps = 1e-6
sh1 = sh0 + eps * direction
IC = ift.GradientNormController(
name='IC',
verbose=False,
iteration_limit=100,
tol_abs_gradnorm=1e-5)
inverter = ift.ConjugateGradient(IC)
S = ift.create_power_operator(hspace, power_spectrum=lambda k: 1.)
energy0 = ift.library.LogNormalWienerFilterEnergy(
position=sh0, d=d, R=R, N=N, S=S, inverter=inverter)
energy1 = ift.library.LogNormalWienerFilterEnergy(
position=sh1, d=d, R=R, N=N, S=S, inverter=inverter)
a = (energy1.value - energy0.value) / eps
b = energy0.gradient.vdot(direction)
tol = 1e-3
assert_allclose(a, b, rtol=tol, atol=tol)
@expand(product([ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[ift.library.Exponential, ift.library.Linear],
[4, 78, 23]))
def testNonlinearMap(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)
P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
def pspec(k): return 1 / (1 + k**2)**dim
pspec = ift.PS_field(pspace, pspec)
A = P.adjoint_times(ift.sqrt(pspec))
n = ift.Field.from_random(domain=space, random_type='normal')
s = ht(xi0 * A)
diag = ift.Field.ones(space) * 10
R = ift.DiagonalOperator(diag)
diag = ift.Field.ones(space)
N = ift.DiagonalOperator(diag)
d = R(f(s)) + n
direction = ift.Field.from_random('normal', hspace)
direction /= np.sqrt(direction.var())
eps = 1e-7
xi1 = xi0 + eps * direction
S = ift.create_power_operator(hspace, power_spectrum=lambda k: 1.)
energy0 = ift.library.NonlinearWienerFilterEnergy(
position=xi0, d=d, Instrument=R, nonlinearity=f, ht=ht, power=A, N=N, S=S)
energy1 = ift.library.NonlinearWienerFilterEnergy(
position=xi1, d=d, Instrument=R, nonlinearity=f, ht=ht, power=A, N=N, S=S)
a = (energy1.value - energy0.value) / eps
b = energy0.gradient.vdot(direction)
tol = 1e-4
assert_allclose(a, b, rtol=tol, atol=tol)
# 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-2017 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 nifty4 as ift
import numpy as np
from itertools import product
from test.common import expand
from numpy.testing import assert_allclose
# TODO Add also other space types
class Noise_Energy_Tests(unittest.TestCase):
@expand(product([ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[ift.library.Exponential, ift.library.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)
P = ift.PowerProjectionOperator(domain=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 = P.adjoint_times(ift.sqrt(tau))
var = 1.
n = ift.Field.from_random(domain=space, random_type='normal', std=np.sqrt(var))
var = ift.Field(n.domain, val=var)
N = ift.DiagonalOperator(var)
eta0 = ift.log(var)
s = ht(xi * A)
diag = ift.Field.ones(space) * 10
R = ift.DiagonalOperator(diag)
diag = ift.Field.ones(space)
d = R(f(s)) + n
alpha = ift.Field(d.domain, val=2.)
q = ift.Field(d.domain, val=1e-5)
direction = ift.Field.from_random('normal', d.domain)
direction /= np.sqrt(direction.var())
eps = 1e-8
eta1 = eta0 + eps * direction
IC = ift.GradientNormController(
name='IC',
verbose=False,
iteration_limit=100,
tol_abs_gradnorm=1e-5)
inverter = ift.ConjugateGradient(IC)
S = ift.create_power_operator(hspace, power_spectrum=lambda k: 1.)
D = ift.library.NonlinearWienerFilterEnergy(
position=xi,
d=d,
Instrument=R,
nonlinearity=f,
ht=ht,
power=A,
N=N,
S=S,
inverter=inverter).curvature
Nsamples = 10
sample_list = [D.generate_posterior_sample() + xi for i in range(Nsamples)]
energy0 = ift.library.NoiseEnergy(
position=eta0, d=d, m=xi, D=D, t=tau, Instrument=R,
alpha=alpha, q=q, Projection=P, nonlinearity=f,
ht=ht, sample_list=sample_list)
energy1 = ift.library.NoiseEnergy(
position=eta1, d=d, m=xi, D=D, t=tau, Instrument=R,
alpha=alpha, q=q, Projection=P, nonlinearity=f,
ht=ht, sample_list=sample_list)
a = (energy1.value - energy0.value) / eps
b = energy0.gradient.vdot(direction)
tol = 1e-5
assert_allclose(a, b, rtol=tol, atol=tol)
# 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