Skip to content

GitLab

Commit 41879da1 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

DiagonalOperator takes and returns diagonal without volume factors

parent 88872f64
Pipeline #19209 failed with stage
Hide whitespace changes
Inline Side-by-side
demos/critical_filtering.py
View file @ 41879da1
......@@ -72,7 +72,7 @@ if __name__ == "__main__":
R = AdjointFFTResponse(fft, Instrument)
noise = 1.
N = ift.DiagonalOperator(ift.Field(s_space,noise).weight(1))
N = ift.DiagonalOperator(ift.Field(s_space,noise))
n = ift.Field.from_random(domain=s_space,
random_type='normal',
std=np.sqrt(noise),
......
demos/log_normal_wiener_filter.py
View file @ 41879da1
......@@ -38,7 +38,7 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R])
# Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise)
N = ift.DiagonalOperator(ndiag)
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
......
demos/paper_demos/cartesian_wiener_filter.py
View file @ 41879da1
......@@ -62,7 +62,7 @@ if __name__ == "__main__":
diagonal = mock_power.power_synthesize_special(spaces=(0, 1))**2
diagonal = diagonal.real
S = ift.DiagonalOperator(diagonal)
S = ift.DiagonalOperator(diagonal.weight(-1))
np.random.seed(10)
......@@ -84,7 +84,7 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R])
# Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise)
N = ift.DiagonalOperator(ndiag)
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise),
......@@ -93,7 +93,7 @@ if __name__ == "__main__":
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(verbose=True, iteration_limit=100)
ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic, inverter=inverter)
......
demos/paper_demos/wiener_filter.py
View file @ 41879da1
......@@ -37,7 +37,7 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R])
# Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise)
N = ift.DiagonalOperator(ndiag)
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
......
demos/wiener_filter_via_curvature.py
View file @ 41879da1
......@@ -47,7 +47,7 @@ if __name__ == "__main__":
data_domain = R.target[0]
R_harmonic = ift.ComposedOperator([fft, R])
N = ift.DiagonalOperator(ift.Field(data_domain,mock_signal.var()/signal_to_noise).weight(1))
N = ift.DiagonalOperator(ift.Field(data_domain,mock_signal.var()/signal_to_noise))
noise = ift.Field.from_random(domain=data_domain,
random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise),
......
demos/wiener_filter_via_hamiltonian.py
View file @ 41879da1
......@@ -62,7 +62,7 @@ if __name__ == "__main__":
# Adding a harmonic transformation to the instrument
R = AdjointFFTResponse(fft, Instrument)
signal_to_noise = 1.
ndiag = ift.Field(s_space, ss.var()/signal_to_noise).weight(1)
ndiag = ift.Field(s_space, ss.var()/signal_to_noise)
N = ift.DiagonalOperator(ndiag)
n = ift.Field.from_random(domain=s_space,
random_type='normal',
......
nifty/field.py
View file @ 41879da1
......@@ -420,7 +420,7 @@ class Field(object):
res *= tmp
return res
def weight(self, power=1, inplace=False, spaces=None):
def weight(self, power=1, spaces=None, out=None):
""" Weights the pixels of `self` with their invidual pixel-volume.
Parameters
......@@ -428,20 +428,25 @@ class Field(object):
power : number
The pixels get weighted with the volume-factor**power.
inplace : boolean
If True, `self` will be weighted and returned. Otherwise, a copy
is made.
spaces : tuple of ints
Determines on which subspace the operation takes place.
out : Field or None
if not None, the result is returned in a new Field
otherwise the contents of "out" are overwritten with the result.
"out" may be identical to "self"!
Returns
-------
out : Field
The weighted field.
"""
new_field = self if inplace else self.copy()
if out is None:
out = self.copy()
else:
if out is not self:
out.copy_content_from(self)
if spaces is None:
spaces = range(len(self.domain))
......@@ -458,12 +463,12 @@ class Field(object):
new_shape[self.domain.axes[ind][0]:
self.domain.axes[ind][-1]+1] = wgt.shape
wgt = wgt.reshape(new_shape)
new_field *= wgt**power
out *= wgt**power
fct = fct**power
if fct != 1.:
new_field *= fct
out *= fct
return new_field
return out
def vdot(self, x=None, spaces=None):
""" Computes the volume-factor-aware dot product of 'self' with x.
......@@ -474,8 +479,8 @@ class Field(object):
The domain of x must contain `self.domain`
spaces : tuple of ints
If the domain of `self` and `x` are not the same, `spaces` specfies
the mapping.
If the domain of `self` and `x` are not the same, `spaces` defines
which domains of `x` are mapped to those of `self`.
Returns
-------
......@@ -487,9 +492,9 @@ class Field(object):
"the NIFTy field class")
# Compute the dot respecting the fact of discrete/continuous spaces
fct = 1.
tmp = self.scalar_weight(spaces)
if tmp is None:
fct = 1.
y = self.weight(power=1)
else:
y = self
......@@ -498,13 +503,14 @@ class Field(object):
if spaces is None:
return fct*dobj.vdot(y.val.ravel(), x.val.ravel())
else:
# create a diagonal operator which is capable of taking care of the
# axes-matching
from .operators.diagonal_operator import DiagonalOperator
diag = DiagonalOperator(y.conjugate(), self.domain,
spaces=spaces, copy=False)
dotted = diag(x)
return fct*dotted.sum(spaces=spaces)
spaces = utilities.cast_iseq_to_tuple(spaces)
active_axes = []
for i in spaces:
active_axes += self.domain.axes[i]
res = 0.
for sl in utilities.get_slice_list(self.shape, active_axes):
res += dobj.vdot(y.val, x.val[sl])
return res*fct
def norm(self):
""" Computes the L2-norm of the field values.
......@@ -515,7 +521,7 @@ class Field(object):
The L2-norm of the field values.
"""
return dobj.sqrt(dobj.abs(self.vdot(x=self)))
return np.sqrt(np.abs(self.vdot(x=self)))
def conjugate(self):
""" Returns the complex conjugate of the field.
......@@ -566,6 +572,10 @@ class Field(object):
def sum(self, spaces=None):
return self._contraction_helper('sum', spaces)
def integrate(self, spaces=None):
tmp = self.weight(1, spaces=spaces)
return tmp.sum(spaces)
def prod(self, spaces=None):
return self._contraction_helper('prod', spaces)
......
nifty/library/critical_filter/critical_power_energy.py
View file @ 41879da1
......@@ -88,7 +88,7 @@ class CriticalPowerEnergy(Energy):
@property
def gradient(self):
gradient = -self._theta.weight(-1)
gradient += (self._rho_prime).weight(-1)
gradient += self._rho_prime.weight(-1)
gradient += self._Tt
gradient = gradient.real
return gradient
......
nifty/operators/diagonal_operator/diagonal_operator.py
View file @ 41879da1
......@@ -37,8 +37,12 @@ class DiagonalOperator(EndomorphicOperator):
----------
diagonal : Field
The diagonal entries of the operator.
copy : boolean
Internal copy of the diagonal (default: True)
domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
The domain on which the Operator's input Field lives.
If None, use the domain of "diagonal".
spaces : tuple of int
The elements of "domain" on which the operator acts.
If None, it acts on all elements.
Attributes
----------
......@@ -60,7 +64,7 @@ class DiagonalOperator(EndomorphicOperator):
# ---Overwritten properties and methods---
def __init__(self, diagonal, domain=None, spaces=None, copy=True):
def __init__(self, diagonal, domain=None, spaces=None):
super(DiagonalOperator, self).__init__()
if not isinstance(diagonal, Field):
......@@ -87,7 +91,7 @@ class DiagonalOperator(EndomorphicOperator):
if diagonal.domain[i] != self._domain[j]:
raise ValueError("domain mismatch")
self._diagonal = diagonal if not copy else diagonal.copy()
self._diagonal = diagonal.weight(1)
self._self_adjoint = None
self._unitary = None
......@@ -103,21 +107,16 @@ class DiagonalOperator(EndomorphicOperator):
def _adjoint_inverse_times(self, x):
return self._times_helper(x, lambda z: z.conjugate().__rtruediv__)
def diagonal(self, copy=True):
def diagonal(self):
""" Returns the diagonal of the Operator.
Parameters
----------
copy : boolean
Whether the returned Field should be copied or not.
Returns
-------
out : Field
The diagonal of the Operator.
"""
return self._diagonal.copy() if copy else self._diagonal
return self._diagonal.weight(-1)
# ---Mandatory properties and methods---
......
nifty/operators/laplace_operator/laplace_operator.py
View file @ 41879da1
......@@ -111,7 +111,7 @@ class LaplaceOperator(EndomorphicOperator):
ret /= np.sqrt(dposc)
ret[prefix + (slice(None, 2),)] = 0.
ret[prefix + (-1,)] = 0.
return Field(self.domain, val=ret).weight(-0.5, spaces=(self._space,))
return Field(self.domain, val=ret)
def _adjoint_times(self, x):
axes = x.domain.axes[self._space]
......@@ -122,7 +122,7 @@ class LaplaceOperator(EndomorphicOperator):
sl_r = prefix + (slice(1, None),) # "right" slice
dpos = self._dpos.reshape((1,)*axis + (nval-1,))
dposc = self._dposc.reshape((1,)*axis + (nval,))
y = x.copy().weight(power=0.5, spaces=(self._space,)).val
y = x.val.copy()
y /= np.sqrt(dposc)
y[prefix + (slice(None, 2),)] = 0.
y[prefix + (-1,)] = 0.
......@@ -131,4 +131,4 @@ class LaplaceOperator(EndomorphicOperator):
ret[sl_l] = deriv
ret[prefix + (-1,)] = 0.
ret[sl_r] -= deriv
return Field(self.domain, val=ret).weight(-1, spaces=(self._space,))
return Field(self.domain, val=ret)
nifty/operators/response_operator/response_operator.py
View file @ 41879da1
......@@ -62,7 +62,7 @@ class ResponseOperator(LinearOperator):
space=spaces[x])
for x in range(nsigma)]
kernel_exposure = [DiagonalOperator(Field(self._domain[spaces[x]],
exposure[x]),
exposure[x]).weight(-1),
domain=self._domain,
spaces=(spaces[x],))
for x in range(nsigma)]
......
nifty/sugar.py
View file @ 41879da1
......@@ -26,11 +26,27 @@ from . import Space,\
FFTOperator,\
sqrt
__all__ = ['create_power_operator',
__all__ = ['create_power_field',
'create_power_operator',
'generate_posterior_sample',
'create_composed_fft_operator']
def create_power_field(domain, power_spectrum, dtype=None):
if not callable(power_spectrum):
raise TypeError("power_spectrum must be callable")
power_domain = PowerSpace(domain)
fp = Field(power_domain, val=power_spectrum(power_domain.k_lengths),
dtype=dtype)
f = fp.power_synthesize_special()
if not issubclass(fp.dtype.type, np.complexfloating):
f = f.real
f **= 2
return f
def create_power_operator(domain, power_spectrum, dtype=None):
""" Creates a diagonal operator with the given power spectrum.
......@@ -53,21 +69,7 @@ def create_power_operator(domain, power_spectrum, dtype=None):
DiagonalOperator : An operator that implements the given power spectrum.
"""
if not callable(power_spectrum):
raise TypeError("power_spectrum must be callable")
power_domain = PowerSpace(domain)
fp = Field(power_domain, val=power_spectrum(power_domain.k_lengths),
dtype=dtype)
f = fp.power_synthesize_special()
if not issubclass(fp.dtype.type, np.complexfloating):
f = f.real
f **= 2
return DiagonalOperator(Field(domain,f).weight(1))
return DiagonalOperator(create_power_field(domain, power_spectrum, dtype))
def generate_posterior_sample(mean, covariance):
""" Generates a posterior sample from a Gaussian distribution with given
......@@ -96,10 +98,10 @@ def generate_posterior_sample(mean, covariance):
R = covariance.R
N = covariance.N
power = sqrt(S.diagonal().power_analyze())
power = sqrt(S.diagonal().weight(1).power_analyze())
mock_signal = power.power_synthesize(real_signal=True)
noise = N.diagonal().weight(-1)
noise = N.diagonal()
mock_noise = Field.from_random(random_type="normal", domain=N.domain,
dtype=noise.dtype.type)
......
test/test_field.py
View file @ 41879da1
......@@ -32,7 +32,7 @@ from nifty2go import Field,\
DomainTuple,\
DiagonalOperator
from nifty2go.sugar import create_power_operator
from nifty2go.sugar import create_power_field
from test.common import expand
......@@ -80,8 +80,8 @@ class Test_Functionality(unittest.TestCase):
sk = fp.power_synthesize(spaces=(0, 1), real_signal=True)
sp = sk.power_analyze(spaces=(0, 1), keep_phase_information=False)
ps1 += sp.sum(spaces=1)/fp2.sum()
ps2 += sp.sum(spaces=0)/fp1.sum()
ps1 += sp.integrate(spaces=1)/fp2.sum()
ps2 += sp.integrate(spaces=0)/fp1.sum()
assert_allclose(ps1.val/samples, fp1.val, rtol=0.2)
assert_allclose(ps2.val/samples, fp2.val, rtol=0.2)
......@@ -104,12 +104,10 @@ class Test_Functionality(unittest.TestCase):
spec2 = lambda k: 42/(1+k)**3
fp2 = Field(p2, val=spec2(p2.k_lengths))
S_1 = create_power_operator(space1, lambda x: np.sqrt(spec1(x)))
S_2 = create_power_operator(space2, lambda x: np.sqrt(spec2(x)))
S_1 = DiagonalOperator(S_1.diagonal().weight(-1),domain=fulldomain,
spaces=0)
S_2 = DiagonalOperator(S_2.diagonal().weight(-1),domain=fulldomain,
spaces=1)
S_1 = create_power_field(space1, lambda x: np.sqrt(spec1(x)))
S_1 = DiagonalOperator(S_1.weight(-1), domain=fulldomain, spaces=0)
S_2 = create_power_field(space2, lambda x: np.sqrt(spec2(x)))
S_2 = DiagonalOperator(S_2.weight(-1), domain=fulldomain, spaces=1)
samples = 500
ps1 = 0.
......@@ -119,8 +117,8 @@ class Test_Functionality(unittest.TestCase):
rand_k = Field.from_random('normal', domain=fulldomain)
sk = S_1.times(S_2.times(rand_k))
sp = sk.power_analyze(spaces=(0, 1), keep_phase_information=False)
ps1 += sp.sum(spaces=1)/fp2.sum()
ps2 += sp.sum(spaces=0)/fp1.sum()
ps1 += sp.integrate(spaces=1)/fp2.sum()
ps2 += sp.integrate(spaces=0)/fp1.sum()
assert_allclose(ps1.val/samples, fp1.val, rtol=0.2)
assert_allclose(ps2.val/samples, fp2.val, rtol=0.2)
......
test/test_minimization/test_minimizers.py
View file @ 41879da1
......@@ -21,7 +21,7 @@ class Test_Minimizers(unittest.TestCase):
starting_point = ift.Field.from_random('normal', domain=space)*10
covariance_diagonal = ift.Field.from_random(
'uniform', domain=space) + 0.5
covariance = ift.DiagonalOperator(covariance_diagonal)
covariance = ift.DiagonalOperator(covariance_diagonal.weight(-1))
required_result = ift.Field(space, val=1.)
IC = ift.GradientNormController(tol_abs_gradnorm=1e-5)
......
test/test_operators/test_diagonal_operator.py
View file @ 41879da1
......@@ -17,8 +17,8 @@ from test.common import expand
class DiagonalOperator_Tests(unittest.TestCase):
spaces = generate_spaces()
@expand(product(spaces, [True, False]))
def test_property(self, space, copy):
@expand(product(spaces))
def test_property(self, space):
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(diag)
if D.domain[0] != space:
......@@ -28,59 +28,59 @@ class DiagonalOperator_Tests(unittest.TestCase):
if D.self_adjoint != True:
raise TypeError
@expand(product(spaces, [True, False]))
def test_times_adjoint(self, space, copy):
@expand(product(spaces))
def test_times_adjoint(self, space):
rand1 = Field.from_random('normal', domain=space)
rand2 = Field.from_random('normal', domain=space)
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(diag, copy=copy)
D = DiagonalOperator(diag)
tt1 = rand1.vdot(D.times(rand2))
tt2 = rand2.vdot(D.times(rand1))
assert_approx_equal(tt1, tt2)
@expand(product(spaces, [True, False]))
def test_times_inverse(self, space, copy):
@expand(product(spaces))
def test_times_inverse(self, space):
rand1 = Field.from_random('normal', domain=space)
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(diag, copy=copy)
D = DiagonalOperator(diag)
tt1 = D.times(D.inverse_times(rand1))
assert_allclose(rand1.val, tt1.val)
@expand(product(spaces, [True, False]))
def test_times(self, space, copy):
@expand(product(spaces))
def test_times(self, space):
rand1 = Field.from_random('normal', domain=space)
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(diag, copy=copy)
D = DiagonalOperator(diag)
tt = D.times(rand1)
assert_equal(tt.domain[0], space)
@expand(product(spaces, [True, False]))
def test_adjoint_times(self, space, copy):
@expand(product(spaces))
def test_adjoint_times(self, space):
rand1 = Field.from_random('normal', domain=space)
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(diag, copy=copy)
D = DiagonalOperator(diag)
tt = D.adjoint_times(rand1)
assert_equal(tt.domain[0], space)
@expand(product(spaces, [True, False]))
def test_inverse_times(self, space, copy):
@expand(product(spaces))
def test_inverse_times(self, space):
rand1 = Field.from_random('normal', domain=space)
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(diag, copy=copy)
D = DiagonalOperator(diag)
tt = D.inverse_times(rand1)
assert_equal(tt.domain[0], space)
@expand(product(spaces, [True, False]))
def test_adjoint_inverse_times(self, space, copy):
@expand(product(spaces))
def test_adjoint_inverse_times(self, space):
rand1 = Field.from_random('normal', domain=space)
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(diag, copy=copy)
D = DiagonalOperator(diag)
tt = D.adjoint_inverse_times(rand1)
assert_equal(tt.domain[0], space)
@expand(product(spaces, [True, False]))
def test_diagonal(self, space, copy):
@expand(product(spaces))
def test_diagonal(self, space):
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(diag, copy=copy)
D = DiagonalOperator(diag)
diag_op = D.diagonal()
assert_allclose(diag.val, diag_op.val)
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