Skip to content

GitLab

Commit 3818a258 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

byebye bare

parent f2b5860d
Pipeline #18528 failed with stage
in 3 minutes and 55 seconds
Hide whitespace changes
Inline Side-by-side
demos/D2O.conf 0 → 100644
View file @ 3818a258
[DEFAULT]
mpi_module = MPI
mpi_init_checks = True
default_distribution_strategy = equal
default_comm = COMM_WORLD
demos/critical_filtering.py
View file @ 3818a258
......@@ -5,7 +5,7 @@ np.random.seed(42)
#np.seterr(all="raise",under="ignore")
def plot_parameters(m, t, p, p_d):
x = ift.log(t.domain[0].kindex)
x = np.log(t.domain[0].kindex)
m = fft.adjoint_times(m)
t = t.val.real
p = p.val.real
......@@ -72,10 +72,10 @@ if __name__ == "__main__":
R = AdjointFFTResponse(fft, Instrument)
noise = 1.
N = ift.DiagonalOperator(s_space, diagonal=noise, bare=True)
N = ift.DiagonalOperator(s_space, ift.Field(s_space,noise).weight(1))
n = ift.Field.from_random(domain=s_space,
random_type='normal',
std=ift.sqrt(noise),
std=np.sqrt(noise),
mean=0)
# Create mock data
......@@ -101,7 +101,7 @@ if __name__ == "__main__":
def ps0(k):
return (1./(1.+k)**2)
t0 = ift.Field(p_space, val=ift.log(1./(1+p_space.kindex)**2))
t0 = ift.Field(p_space, val=np.log(1./(1+p_space.kindex)**2))
for i in range(500):
S0 = ift.create_power_operator(h_space, power_spectrum=ps0)
......
demos/log_normal_wiener_filter.py
View file @ 3818a258
......@@ -38,7 +38,8 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R], default_spaces=[0, 0])
# Setting up the noise covariance and drawing a random noise realization
N = ift.DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise, bare=True)
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
N = ift.DiagonalOperator(data_domain, ndiag)
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(ift.exp(mock_signal)) + noise #|\label{code:wf_mock_data}|
......
demos/paper_demos/cartesian_wiener_filter.py
View file @ 3818a258
......@@ -81,8 +81,8 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R], default_spaces=(0, 1, 0, 1))
# Setting up the noise covariance and drawing a random noise realization
N = ift.DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise,
bare=True)
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
N = ift.DiagonalOperator(data_domain, 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/wiener_filter.py
View file @ 3818a258
......@@ -21,7 +21,7 @@ if __name__ == "__main__":
signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
harmonic_space = signal_space.get_default_codomain()
fft = ift.FFTOperator(harmonic_space, target=signal_space)
power_space = ift.PowerSpace(harmonic_space)
power_space = ift.PowerSpace(harmonic_space,binbounds=ift.PowerSpace.useful_binbounds(harmonic_space,logarithmic=True))
# Creating the mock signal |\label{code:wf_mock_signal}|
S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum)
......@@ -37,14 +37,15 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R], default_spaces=[0, 0])
# Setting up the noise covariance and drawing a random noise realization
N = ift.DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise, bare=True)
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
N = ift.DiagonalOperator(data_domain, ndiag)
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(mock_signal) + noise #|\label{code:wf_mock_data}|
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.DefaultIterationController(verbose=False,tol_abs_gradnorm=0.1)
ctrl = ift.DefaultIterationController(verbose=True,tol_abs_gradnorm=0.1,iteration_limit=10)
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic,inverter=inverter)
m_k = wiener_curvature.inverse_times(j) #|\label{code:wf_wiener_filter}|
......
demos/wiener_filter_via_curvature.py
View file @ 3818a258
......@@ -48,8 +48,7 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R], default_spaces=[0, 0])
N = ift.DiagonalOperator(data_domain,
diagonal=mock_signal.var()/signal_to_noise,
bare=True)
diagonal=ift.Field(data_domain,mock_signal.var()/signal_to_noise).weight(1))
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 @ 3818a258
......@@ -62,7 +62,8 @@ if __name__ == "__main__":
# Adding a harmonic transformation to the instrument
R = AdjointFFTResponse(fft, Instrument)
signal_to_noise = 1.
N = ift.DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
ndiag = ift.Field(s_space, ss.var()/signal_to_noise).weight(1)
N = ift.DiagonalOperator(s_space, ndiag)
n = ift.Field.from_random(domain=s_space,
random_type='normal',
std=ss.std()/np.sqrt(signal_to_noise),
......
nifty/__init__.py
View file @ 3818a258
......@@ -47,3 +47,5 @@ from .sugar import *
from . import plotting
from . import library
from . import dobj
nifty/field.py
View file @ 3818a258
......@@ -473,7 +473,7 @@ class Field(object):
return new_field
def vdot(self, x=None, spaces=None, bare=False):
def vdot(self, x=None, spaces=None):
""" Computes the volume-factor-aware dot product of 'self' with x.
Parameters
......@@ -485,9 +485,6 @@ class Field(object):
If the domain of `self` and `x` are not the same, `spaces` specfies
the mapping.
bare : boolean
If true, no volume factors will be included in the computation.
Returns
-------
out : float, complex
......@@ -499,15 +496,12 @@ class Field(object):
# Compute the dot respecting the fact of discrete/continuous spaces
fct = 1.
if bare:
y = self
tmp = self.scalar_weight(spaces)
if tmp is None:
y = self.weight(power=1)
else:
tmp = self.scalar_weight(spaces)
if tmp is None:
y = self.weight(power=1)
else:
y = self
fct = tmp
y = self
fct = tmp
if spaces is None:
return fct*np.vdot(y.val.ravel(), x.val.ravel())
......
nifty/library/critical_filter/critical_power_energy.py
View file @ 3818a258
......@@ -81,7 +81,7 @@ class CriticalPowerEnergy(Energy):
@property
def value(self):
energy = self._theta.sum()
energy += self.position.vdot(self._rho_prime, bare=True)
energy += self.position.weight(-1).vdot(self._rho_prime)
energy += 0.5 * self.position.vdot(self._Tt)
return energy.real
......
nifty/operators/diagonal_operator/diagonal_operator.py
View file @ 3818a258
......@@ -39,9 +39,6 @@ class DiagonalOperator(EndomorphicOperator):
The domain on which the Operator's input Field lives.
diagonal : {scalar, list, array, Field}
The diagonal entries of the operator.
bare : boolean
Indicates whether the input for the diagonal is bare or not
(default: False).
copy : boolean
Internal copy of the diagonal (default: True)
default_spaces : tuple of ints *optional*
......@@ -63,16 +60,6 @@ class DiagonalOperator(EndomorphicOperator):
Raises
------
Notes
-----
The ambiguity of bare or non-bare diagonal entries is based on the choice
of a matrix representation of the operator in question. The naive choice
of absorbing the volume weights into the matrix leads to a matrix-vector
calculus with the non-bare entries which seems intuitive, though.
The choice of keeping matrix entries and volume weights separate
deals with the bare entries that allow for correct interpretation
of the matrix entries; e.g., as variance in case of an covariance operator.
See Also
--------
EndomorphicOperator
......@@ -81,7 +68,7 @@ class DiagonalOperator(EndomorphicOperator):
# ---Overwritten properties and methods---
def __init__(self, domain=(), diagonal=None, bare=False, copy=True,
def __init__(self, domain=(), diagonal=None, copy=True,
default_spaces=None):
super(DiagonalOperator, self).__init__(default_spaces)
......@@ -89,7 +76,7 @@ class DiagonalOperator(EndomorphicOperator):
self._self_adjoint = None
self._unitary = None
self.set_diagonal(diagonal=diagonal, bare=bare, copy=copy)
self.set_diagonal(diagonal=diagonal, copy=copy)
def _times(self, x, spaces):
return self._times_helper(x, spaces, operation=lambda z: z.__mul__)
......@@ -107,13 +94,11 @@ class DiagonalOperator(EndomorphicOperator):
operation=lambda z:
z.conjugate().__rtruediv__)
def diagonal(self, bare=False, copy=True):
def diagonal(self, copy=True):
""" Returns the diagonal of the Operator.
Parameters
----------
bare : boolean
Whether the returned Field values should be bare or not.
copy : boolean
Whether the returned Field should be copied or not.
......@@ -123,28 +108,18 @@ class DiagonalOperator(EndomorphicOperator):
The diagonal of the Operator.
"""
if bare:
return self._diagonal.weight(power=-1)
elif copy:
return self._diagonal.copy()
else:
return self._diagonal
return self._diagonal.copy() if copy else self._diagonal
def inverse_diagonal(self, bare=False):
def inverse_diagonal(self):
""" Returns the inverse-diagonal of the operator.
Parameters
----------
bare : boolean
Whether the returned Field values should be bare or not.
Returns
-------
out : Field
The inverse of the diagonal of the Operator.
"""
return 1./self.diagonal(bare=bare, copy=False)
return 1./self.diagonal(copy=False)
# ---Mandatory properties and methods---
......@@ -169,16 +144,13 @@ class DiagonalOperator(EndomorphicOperator):
# ---Added properties and methods---
def set_diagonal(self, diagonal, bare=False, copy=True):
def set_diagonal(self, diagonal, copy=True):
""" Sets the diagonal of the Operator.
Parameters
----------
diagonal : {scalar, list, array, Field}
The diagonal entries of the operator.
bare : boolean
Indicates whether the input for the diagonal is bare or not
(default: False).
copy : boolean
Specifies if a copy of the input shall be made (default: True).
......@@ -187,13 +159,6 @@ class DiagonalOperator(EndomorphicOperator):
# use the casting functionality from Field to process `diagonal`
f = Field(domain=self.domain, val=diagonal, copy=copy)
# weight if the given values were `bare` is True
# do inverse weightening if the other way around
if bare:
# If `copy` is True, we won't change external data by weightening
# Otherwise, inplace weighting would change the external field
f.weight(inplace=copy)
# Reset the self_adjoint property:
self._self_adjoint = None
......
nifty/probing/mixin_classes/trace_prober_mixin.py
View file @ 3818a258
......@@ -37,9 +37,9 @@ class TraceProberMixin(object):
def finish_probe(self, probe, pre_result):
if self.__evaluate_probe_in_signal_space:
fft = create_composed_fft_operator(self._domain, all_to='position')
result = fft(probe[1]).vdot(fft(pre_result), bare=True)
result = fft(probe[1]).weight(-1).vdot(fft(pre_result))
else:
result = probe[1].vdot(pre_result, bare=True)
result = probe[1].weight(-1).vdot(pre_result)
self.__sum_of_probings += result
if self.compute_variance:
......
nifty/sugar.py
View file @ 3818a258
......@@ -66,7 +66,7 @@ def create_power_operator(domain, power_spectrum, dtype=None):
f = f.real
f **= 2
return DiagonalOperator(domain, diagonal=f, bare=True)
return DiagonalOperator(domain, diagonal=Field(domain,f).weight(1))
def generate_posterior_sample(mean, covariance):
......@@ -99,10 +99,10 @@ def generate_posterior_sample(mean, covariance):
power = sqrt(S.diagonal().power_analyze())
mock_signal = power.power_synthesize(real_signal=True)
noise = N.diagonal(bare=True)
noise = N.diagonal().weight(-1)
mock_noise = Field.from_random(random_type="normal", domain=N.domain,
dtype=noise.dtype)
dtype=noise.dtype.type)
mock_noise *= sqrt(noise)
mock_data = R(mock_signal) + mock_noise
......
test/test_operators/test_diagonal_operator.py
View file @ 3818a258
......@@ -17,8 +17,8 @@ from test.common import expand
class DiagonalOperator_Tests(unittest.TestCase):
spaces = generate_spaces()
@expand(product(spaces, [True, False], [True, False]))
def test_property(self, space, bare, copy):
@expand(product(spaces, [True, False]))
def test_property(self, space, copy):
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag)
if D.domain[0] != space:
......@@ -28,53 +28,53 @@ class DiagonalOperator_Tests(unittest.TestCase):
if D.self_adjoint != True:
raise TypeError
@expand(product(spaces, [True, False], [True, False]))
def test_times_adjoint(self, space, bare, copy):
@expand(product(spaces, [True, False]))
def test_times_adjoint(self, space, copy):
rand1 = Field.from_random('normal', domain=space)
rand2 = Field.from_random('normal', domain=space)
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
D = DiagonalOperator(space, diagonal=diag, copy=copy)
tt1 = rand1.vdot(D.times(rand2))
tt2 = rand2.vdot(D.times(rand1))
assert_approx_equal(tt1, tt2)
@expand(product(spaces, [True, False], [True, False]))
def test_times_inverse(self, space, bare, copy):
@expand(product(spaces, [True, False]))
def test_times_inverse(self, space, copy):
rand1 = Field.from_random('normal', domain=space)
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
D = DiagonalOperator(space, diagonal=diag, copy=copy)
tt1 = D.times(D.inverse_times(rand1))
assert_allclose(rand1.val, tt1.val)
@expand(product(spaces, [True, False], [True, False]))
def test_times(self, space, bare, copy):
@expand(product(spaces, [True, False]))
def test_times(self, space, copy):
rand1 = Field.from_random('normal', domain=space)
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
D = DiagonalOperator(space, diagonal=diag, copy=copy)
tt = D.times(rand1)
assert_equal(tt.domain[0], space)
@expand(product(spaces, [True, False], [True, False]))
def test_adjoint_times(self, space, bare, copy):
@expand(product(spaces, [True, False]))
def test_adjoint_times(self, space, copy):
rand1 = Field.from_random('normal', domain=space)
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
D = DiagonalOperator(space, diagonal=diag, copy=copy)
tt = D.adjoint_times(rand1)
assert_equal(tt.domain[0], space)
@expand(product(spaces, [True, False], [True, False]))
def test_inverse_times(self, space, bare, copy):
@expand(product(spaces, [True, False]))
def test_inverse_times(self, space, copy):
rand1 = Field.from_random('normal', domain=space)
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
D = DiagonalOperator(space, diagonal=diag, copy=copy)
tt = D.inverse_times(rand1)
assert_equal(tt.domain[0], space)
@expand(product(spaces, [True, False], [True, False]))
def test_adjoint_inverse_times(self, space, bare, copy):
@expand(product(spaces, [True, False]))
def test_adjoint_inverse_times(self, space, copy):
rand1 = Field.from_random('normal', domain=space)
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
D = DiagonalOperator(space, diagonal=diag, copy=copy)
tt = D.adjoint_inverse_times(rand1)
assert_equal(tt.domain[0], space)
......
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