Commit 0119d765 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

introduce a few numpy-like convenience functions for Field

parent ca7a12ac
Pipeline #19500 passed with stage
in 4 minutes and 47 seconds
...@@ -72,7 +72,7 @@ if __name__ == "__main__": ...@@ -72,7 +72,7 @@ if __name__ == "__main__":
R = AdjointFFTResponse(fft, Instrument) R = AdjointFFTResponse(fft, Instrument)
noise = 1. noise = 1.
N = ift.DiagonalOperator(ift.Field(s_space,noise)) N = ift.DiagonalOperator(ift.Field.full(s_space,noise))
n = ift.Field.from_random(domain=s_space, n = ift.Field.from_random(domain=s_space,
random_type='normal', random_type='normal',
std=np.sqrt(noise), std=np.sqrt(noise),
...@@ -96,7 +96,7 @@ if __name__ == "__main__": ...@@ -96,7 +96,7 @@ if __name__ == "__main__":
minimizer3 = ift.SteepestDescent(IC3) minimizer3 = ift.SteepestDescent(IC3)
# Set starting position # Set starting position
flat_power = ift.Field(p_space, val=1e-8) flat_power = ift.Field.full(p_space, 1e-8)
m0 = flat_power.power_synthesize(real_signal=True) m0 = flat_power.power_synthesize(real_signal=True)
def ps0(k): def ps0(k):
......
...@@ -28,7 +28,7 @@ if __name__ == "__main__": ...@@ -28,7 +28,7 @@ if __name__ == "__main__":
mock_signal = fft(mock_power.power_synthesize(real_signal=True)) mock_signal = fft(mock_power.power_synthesize(real_signal=True))
# Setting up an exemplary response # Setting up an exemplary response
mask = ift.Field(signal_space, val=1.) mask = ift.Field.ones(signal_space)
N10 = int(N_pixels/10) N10 = int(N_pixels/10)
#mask.val[N10*5:N10*9, N10*5:N10*9] = 0. #mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}| R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}|
...@@ -36,14 +36,14 @@ if __name__ == "__main__": ...@@ -36,14 +36,14 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R]) R_harmonic = ift.ComposedOperator([fft, R])
# Setting up the noise covariance and drawing a random noise realization # Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise) ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise)
N = ift.DiagonalOperator(ndiag) N = ift.DiagonalOperator(ndiag)
noise = ift.Field.from_random(domain=data_domain, random_type='normal', noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(ift.exp(mock_signal)) + noise #|\label{code:wf_mock_data}| data = R(ift.exp(mock_signal)) + noise #|\label{code:wf_mock_data}|
# Wiener filter # Wiener filter
m0 = ift.Field(harmonic_space, val=0.) m0 = ift.Field.zeros(harmonic_space)
ctrl = ift.GradientNormController(verbose=False,tol_abs_gradnorm=1) ctrl = ift.GradientNormController(verbose=False,tol_abs_gradnorm=1)
ctrl2 = ift.GradientNormController(verbose=True,tol_abs_gradnorm=0.1, name="outer") ctrl2 = ift.GradientNormController(verbose=True,tol_abs_gradnorm=0.1, name="outer")
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
......
...@@ -68,11 +68,11 @@ if __name__ == "__main__": ...@@ -68,11 +68,11 @@ if __name__ == "__main__":
# Setting up a exemplary response # Setting up a exemplary response
N1_10 = int(N_pixels_1/10) N1_10 = int(N_pixels_1/10)
mask_1 = ift.Field(signal_space_1, val=1.) mask_1 = ift.Field.ones(signal_space_1)
mask_1.val[N1_10*7:N1_10*9] = 0. mask_1.val[N1_10*7:N1_10*9] = 0.
N2_10 = int(N_pixels_2/10) N2_10 = int(N_pixels_2/10)
mask_2 = ift.Field(signal_space_2, val=1.) mask_2 = ift.Field.ones(signal_space_2)
mask_2.val[N2_10*7:N2_10*9] = 0. mask_2.val[N2_10*7:N2_10*9] = 0.
R = ift.ResponseOperator(signal_domain,spaces=(0,1), R = ift.ResponseOperator(signal_domain,spaces=(0,1),
...@@ -82,7 +82,7 @@ if __name__ == "__main__": ...@@ -82,7 +82,7 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R]) R_harmonic = ift.ComposedOperator([fft, R])
# Setting up the noise covariance and drawing a random noise realization # Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise) ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise)
N = ift.DiagonalOperator(ndiag) N = ift.DiagonalOperator(ndiag)
noise = ift.Field.from_random(domain=data_domain, random_type='normal', noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), std=mock_signal.std()/np.sqrt(signal_to_noise),
......
...@@ -27,7 +27,7 @@ if __name__ == "__main__": ...@@ -27,7 +27,7 @@ if __name__ == "__main__":
mock_signal = fft(mock_power.power_synthesize(real_signal=True)) mock_signal = fft(mock_power.power_synthesize(real_signal=True))
# Setting up an exemplary response # Setting up an exemplary response
mask = ift.Field(signal_space, val=1.) mask = ift.Field.ones(signal_space)
N10 = int(N_pixels/10) N10 = int(N_pixels/10)
mask.val[N10*5:N10*9, N10*5:N10*9] = 0. mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}| R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}|
...@@ -35,7 +35,7 @@ if __name__ == "__main__": ...@@ -35,7 +35,7 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R]) R_harmonic = ift.ComposedOperator([fft, R])
# Setting up the noise covariance and drawing a random noise realization # Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise) ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise)
N = ift.DiagonalOperator(ndiag) N = ift.DiagonalOperator(ndiag)
noise = ift.Field.from_random(domain=data_domain, random_type='normal', noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
......
...@@ -47,7 +47,7 @@ if __name__ == "__main__": ...@@ -47,7 +47,7 @@ if __name__ == "__main__":
data_domain = R.target[0] data_domain = R.target[0]
R_harmonic = ift.ComposedOperator([fft, R]) R_harmonic = ift.ComposedOperator([fft, R])
N = ift.DiagonalOperator(ift.Field(data_domain,mock_signal.var()/signal_to_noise)) N = ift.DiagonalOperator(ift.Field.full(data_domain,mock_signal.var()/signal_to_noise))
noise = ift.Field.from_random(domain=data_domain, noise = ift.Field.from_random(domain=data_domain,
random_type='normal', random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), std=mock_signal.std()/np.sqrt(signal_to_noise),
......
...@@ -56,13 +56,13 @@ if __name__ == "__main__": ...@@ -56,13 +56,13 @@ if __name__ == "__main__":
# Choosing the measurement instrument # Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.05) # Instrument = SmoothingOperator(s_space, sigma=0.05)
Instrument = ift.DiagonalOperator(ift.Field(s_space, 1.)) Instrument = ift.DiagonalOperator(ift.Field.ones(s_space))
# Instrument._diagonal.val[200:400, 200:400] = 0 # Instrument._diagonal.val[200:400, 200:400] = 0
# Adding a harmonic transformation to the instrument # Adding a harmonic transformation to the instrument
R = AdjointFFTResponse(fft, Instrument) R = AdjointFFTResponse(fft, Instrument)
signal_to_noise = 1. signal_to_noise = 1.
ndiag = ift.Field(s_space, ss.var()/signal_to_noise) ndiag = ift.Field.full(s_space, ss.var()/signal_to_noise)
N = ift.DiagonalOperator(ndiag) N = ift.DiagonalOperator(ndiag)
n = ift.Field.from_random(domain=s_space, n = ift.Field.from_random(domain=s_space,
random_type='normal', random_type='normal',
...@@ -78,7 +78,7 @@ if __name__ == "__main__": ...@@ -78,7 +78,7 @@ if __name__ == "__main__":
ctrl = ift.GradientNormController(verbose=True,tol_abs_gradnorm=0.1) ctrl = ift.GradientNormController(verbose=True,tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
# Setting starting position # Setting starting position
m0 = ift.Field(h_space, val=.0) m0 = ift.Field.zeros(h_space)
# Initializing the Wiener Filter energy # Initializing the Wiener Filter energy
energy = ift.library.WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S, energy = ift.library.WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S,
...@@ -88,8 +88,8 @@ if __name__ == "__main__": ...@@ -88,8 +88,8 @@ if __name__ == "__main__":
# Solving the problem analytically # Solving the problem analytically
m0 = D0.inverse_times(j) m0 = D0.inverse_times(j)
sample_variance = ift.Field(sh.domain, val=0.) sample_variance = ift.Field.zeros(sh.domain)
sample_mean = ift.Field(sh.domain, val=0.) sample_mean = ift.Field.zeros(sh.domain)
# sampling the uncertainty map # sampling the uncertainty map
n_samples = 50 n_samples = 50
......
...@@ -91,6 +91,54 @@ class Field(object): ...@@ -91,6 +91,54 @@ class Field(object):
else: else:
raise TypeError("unknown source type") raise TypeError("unknown source type")
@staticmethod
def full(domain, val, dtype=None):
if not np.isscalar(val):
raise TypeError("val must be a scalar")
return Field(DomainTuple.make(domain), val, dtype)
@staticmethod
def ones(domain, dtype=None):
return Field(DomainTuple.make(domain), 1., dtype)
@staticmethod
def zeros(domain, dtype=None):
return Field(DomainTuple.make(domain), 0., dtype)
@staticmethod
def empty(domain, dtype=None):
return Field(DomainTuple.make(domain), None, dtype)
@staticmethod
def full_like(field, val, dtype=None):
if not isinstance(field, Field):
raise TypeError("field must be of Field type")
return Field.full(field.domain, val, dtype)
@staticmethod
def zeros_like(field, dtype=None):
if not isinstance(field, Field):
raise TypeError("field must be of Field type")
if dtype is None:
dtype = field.dtype
return Field.zeros(field.domain, dtype)
@staticmethod
def ones_like(field, dtype=None):
if not isinstance(field, Field):
raise TypeError("field must be of Field type")
if dtype is None:
dtype = field.dtype
return Field.ones(field.domain, dtype)
@staticmethod
def empty_like(field, dtype=None):
if not isinstance(field, Field):
raise TypeError("field must be of Field type")
if dtype is None:
dtype = field.dtype
return Field.empty(field.domain, dtype)
@staticmethod @staticmethod
def _parse_domain(domain, val=None): def _parse_domain(domain, val=None):
if domain is None: if domain is None:
...@@ -225,7 +273,7 @@ class Field(object): ...@@ -225,7 +273,7 @@ class Field(object):
def _single_power_analyze(field, idx, binbounds): def _single_power_analyze(field, idx, binbounds):
from .operators.power_projection_operator import PowerProjectionOperator from .operators.power_projection_operator import PowerProjectionOperator
power_domain = PowerSpace(field.domain[idx], binbounds) power_domain = PowerSpace(field.domain[idx], binbounds)
ppo = PowerProjectionOperator(field.domain,idx,power_domain) ppo = PowerProjectionOperator(field.domain, power_domain, idx)
return ppo(field) return ppo(field)
def _compute_spec(self, spaces): def _compute_spec(self, spaces):
...@@ -242,7 +290,7 @@ class Field(object): ...@@ -242,7 +290,7 @@ class Field(object):
spec = sqrt(self) spec = sqrt(self)
for i in spaces: for i in spaces:
result_domain[i] = self.domain[i].harmonic_partner result_domain[i] = self.domain[i].harmonic_partner
ppo = PowerProjectionOperator(result_domain,i,self.domain[i]) ppo = PowerProjectionOperator(result_domain, self.domain[i], i)
spec = ppo.adjoint_times(spec) spec = ppo.adjoint_times(spec)
return spec return spec
......
...@@ -99,7 +99,7 @@ class DirectSmoothingOperator(EndomorphicOperator): ...@@ -99,7 +99,7 @@ class DirectSmoothingOperator(EndomorphicOperator):
distances = np.log(np.maximum(distances, 1e-15)) distances = np.log(np.maximum(distances, 1e-15))
ibegin, nval, wgt = self._precompute(distances) ibegin, nval, wgt = self._precompute(distances)
res = Field(x.domain, dtype=x.dtype) res = Field.empty_like(x)
for sl in utilities.get_slice_list(x.val.shape, (axis,)): for sl in utilities.get_slice_list(x.val.shape, (axis,)):
inp = x.val[sl] inp = x.val[sl]
out = np.zeros(inp.shape[0], dtype=inp.dtype) out = np.zeros(inp.shape[0], dtype=inp.dtype)
......
...@@ -50,7 +50,7 @@ class InvertibleOperatorMixin(object): ...@@ -50,7 +50,7 @@ class InvertibleOperatorMixin(object):
if self.__forward_x0 is not None: if self.__forward_x0 is not None:
x0 = self.__forward_x0 x0 = self.__forward_x0
else: else:
x0 = Field(self.target, val=0., dtype=x.dtype) x0 = Field.zeros(self.target, dtype=x.dtype)
(result, convergence) = self.__inverter(QuadraticEnergy( (result, convergence) = self.__inverter(QuadraticEnergy(
A=self.inverse_times, A=self.inverse_times,
...@@ -62,7 +62,7 @@ class InvertibleOperatorMixin(object): ...@@ -62,7 +62,7 @@ class InvertibleOperatorMixin(object):
if self.__backward_x0 is not None: if self.__backward_x0 is not None:
x0 = self.__backward_x0 x0 = self.__backward_x0
else: else:
x0 = Field(self.domain, val=0., dtype=x.dtype) x0 = Field.zeros(self.domain, dtype=x.dtype)
(result, convergence) = self.__inverter(QuadraticEnergy( (result, convergence) = self.__inverter(QuadraticEnergy(
A=self.adjoint_inverse_times, A=self.adjoint_inverse_times,
...@@ -74,7 +74,7 @@ class InvertibleOperatorMixin(object): ...@@ -74,7 +74,7 @@ class InvertibleOperatorMixin(object):
if self.__backward_x0 is not None: if self.__backward_x0 is not None:
x0 = self.__backward_x0 x0 = self.__backward_x0
else: else:
x0 = Field(self.domain, val=0., dtype=x.dtype) x0 = Field.zeros(self.domain, dtype=x.dtype)
(result, convergence) = self.__inverter(QuadraticEnergy( (result, convergence) = self.__inverter(QuadraticEnergy(
A=self.times, A=self.times,
...@@ -86,7 +86,7 @@ class InvertibleOperatorMixin(object): ...@@ -86,7 +86,7 @@ class InvertibleOperatorMixin(object):
if self.__forward_x0 is not None: if self.__forward_x0 is not None:
x0 = self.__forward_x0 x0 = self.__forward_x0
else: else:
x0 = Field(self.target, val=0., dtype=x.dtype) x0 = Field.zeros(self.target, dtype=x.dtype)
(result, convergence) = self.__inverter(QuadraticEnergy( (result, convergence) = self.__inverter(QuadraticEnergy(
A=self.adjoint_times, A=self.adjoint_times,
......
...@@ -23,7 +23,7 @@ from .. import dobj ...@@ -23,7 +23,7 @@ from .. import dobj
import numpy as np import numpy as np
class PowerProjectionOperator(LinearOperator): class PowerProjectionOperator(LinearOperator):
def __init__(self, domain, space=None, power_space=None): def __init__(self, domain, power_space=None, space=None):
super(PowerProjectionOperator, self).__init__() super(PowerProjectionOperator, self).__init__()
# Initialize domain and target # Initialize domain and target
......
...@@ -61,7 +61,7 @@ class SmoothnessOperator(EndomorphicOperator): ...@@ -61,7 +61,7 @@ class SmoothnessOperator(EndomorphicOperator):
result = self._laplace.adjoint_times(self._laplace(x)) result = self._laplace.adjoint_times(self._laplace(x))
result *= self._strength**2 result *= self._strength**2
else: else:
result = Field(x.domain, 0., x.dtype) result = Field.zeros_like(x)
return result return result
# ---Added properties and methods--- # ---Added properties and methods---
......
...@@ -142,7 +142,7 @@ class RGSpace(Space): ...@@ -142,7 +142,7 @@ class RGSpace(Space):
tmp[t2] = True tmp[t2] = True
return np.sqrt(np.nonzero(tmp)[0])*self.distances[0] return np.sqrt(np.nonzero(tmp)[0])*self.distances[0]
else: # do it the hard way else: # do it the hard way
tmp = self.get_k_length_array().unique() # expensive! tmp = np.unique(self.get_k_length_array()) # expensive!
tol = 1e-12*tmp[-1] tol = 1e-12*tmp[-1]
# remove all points that are closer than tol to their right # remove all points that are closer than tol to their right
# neighbors. # neighbors.
......
...@@ -22,7 +22,7 @@ class Test_Minimizers(unittest.TestCase): ...@@ -22,7 +22,7 @@ class Test_Minimizers(unittest.TestCase):
covariance_diagonal = ift.Field.from_random( covariance_diagonal = ift.Field.from_random(
'uniform', domain=space) + 0.5 'uniform', domain=space) + 0.5
covariance = ift.DiagonalOperator(covariance_diagonal.weight(-1)) covariance = ift.DiagonalOperator(covariance_diagonal.weight(-1))
required_result = ift.Field(space, val=1.) required_result = ift.Field.ones(space, dtype=np.float64)
IC = ift.GradientNormController(tol_abs_gradnorm=1e-5) IC = ift.GradientNormController(tol_abs_gradnorm=1e-5)
minimizer = minimizer_class(controller=IC) minimizer = minimizer_class(controller=IC)
......
...@@ -62,7 +62,7 @@ class SmoothingOperator_Tests(unittest.TestCase): ...@@ -62,7 +62,7 @@ class SmoothingOperator_Tests(unittest.TestCase):
@expand(product(spaces, [0., .5, 5.])) @expand(product(spaces, [0., .5, 5.]))
def test_times(self, space, sigma): def test_times(self, space, sigma):
op = FFTSmoothingOperator(space, sigma=sigma) op = FFTSmoothingOperator(space, sigma=sigma)
rand1 = Field(space, val=0.) rand1 = Field.zeros(space)
rand1.val[0] = 1. rand1.val[0] = 1.
tt1 = op.times(rand1) tt1 = op.times(rand1)
assert_allclose(1, tt1.sum()) assert_allclose(1, tt1.sum())
......
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