Commit 9fed8b01 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

stage 1/n (broken)

parent c6aa285b
Pipeline #23672 failed with stage
in 4 minutes and 1 second
...@@ -4,8 +4,8 @@ import numericalunits as nu ...@@ -4,8 +4,8 @@ import numericalunits as nu
if __name__ == "__main__": if __name__ == "__main__":
# In MPI mode, the random seed for numericalunits must be set by hand # In MPI mode, the random seed for numericalunits must be set by hand
nu.reset_units(43) #nu.reset_units(43)
dimensionality = 2 dimensionality = 1
np.random.seed(43) np.random.seed(43)
# Setting up variable parameters # Setting up variable parameters
...@@ -49,7 +49,7 @@ if __name__ == "__main__": ...@@ -49,7 +49,7 @@ if __name__ == "__main__":
mock_harmonic = ift.power_synthesize(mock_power, real_signal=True) mock_harmonic = ift.power_synthesize(mock_power, real_signal=True)
mock_signal = fft(mock_harmonic) mock_signal = fft(mock_harmonic)
exposure = 1. exposure = 1./nu.K
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
exposure=(exposure,)) exposure=(exposure,))
data_domain = R.target[0] data_domain = R.target[0]
...@@ -57,17 +57,18 @@ if __name__ == "__main__": ...@@ -57,17 +57,18 @@ if __name__ == "__main__":
N = ift.DiagonalOperator( N = ift.DiagonalOperator(
ift.Field.full(data_domain, ift.Field.full(data_domain,
mock_signal.var()/signal_to_noise).weight(1)) mock_signal.var()/signal_to_noise))
noise = ift.Field.from_random( noise = ift.Field.from_random(
domain=data_domain, random_type='normal', 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(mock_signal) + noise data = R(mock_signal) #+ noise
print data.val[5]
# Wiener filter print mock_signal.val[5]/nu.K
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data)) j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController( ctrl = ift.GradientNormController(
verbose=True, tol_abs_gradnorm=1e-4/nu.K/(nu.m**(0.5*dimensionality))) verbose=True, tol_abs_gradnorm=1e-4/nu.K)
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature( wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R_harmonic, inverter=inverter) S=S, N=N, R=R_harmonic, inverter=inverter)
...@@ -79,7 +80,7 @@ if __name__ == "__main__": ...@@ -79,7 +80,7 @@ if __name__ == "__main__":
ift.plotting.plot(ift.Field(sspace2, mock_signal.val)/nu.K, ift.plotting.plot(ift.Field(sspace2, mock_signal.val)/nu.K,
name="mock_signal.png") name="mock_signal.png")
data = ift.dobj.to_global_data(data.val).reshape(sspace2.shape)/nu.K data = ift.dobj.to_global_data(data.val).reshape(sspace2.shape)
data = ift.Field(sspace2, val=ift.dobj.from_global_data(data))/nu.K data = ift.Field(sspace2, val=ift.dobj.from_global_data(data))
ift.plotting.plot(ift.Field(sspace2, val=data), name="data.png") ift.plotting.plot(ift.Field(sspace2, val=data), name="data.png")
ift.plotting.plot(ift.Field(sspace2, m_s.val)/nu.K, name="map.png") ift.plotting.plot(ift.Field(sspace2, m_s.val)/nu.K, name="map.png")
...@@ -331,16 +331,8 @@ class Field(object): ...@@ -331,16 +331,8 @@ class Field(object):
spaces = utilities.parse_spaces(spaces, ndom) spaces = utilities.parse_spaces(spaces, ndom)
if len(spaces) == ndom: if len(spaces) == ndom:
tmp = self.scalar_weight(spaces) return dobj.vdot(self.val, x.val)
if tmp is None: raise NotImplementedError
fct = 1.
y = self.weight(power=1)
else:
y = self
fct = tmp
return fct*dobj.vdot(y.val, x.val)
# If we arrive here, we have to do a partial dot product. # If we arrive here, we have to do a partial dot product.
# For the moment, do this the explicit, non-optimized way # For the moment, do this the explicit, non-optimized way
return (self.conjugate()*x).integrate(spaces=spaces) return (self.conjugate()*x).integrate(spaces=spaces)
......
...@@ -28,8 +28,8 @@ class GeometryRemover(LinearOperator): ...@@ -28,8 +28,8 @@ class GeometryRemover(LinearOperator):
def apply(self, x, mode): def apply(self, x, mode):
self._check_input(x, mode) self._check_input(x, mode)
if mode == self.TIMES: if mode == self.TIMES:
return Field(self._target, val=x.val) return Field(self._target, val=x.weight(1).val)
return Field(self._domain, val=x.val).weight(power=-1) return Field(self._domain, val=x.val).weight(1)
def ResponseOperator(domain, sigma, exposure): def ResponseOperator(domain, sigma, exposure):
......
...@@ -173,8 +173,8 @@ def power_synthesize(field, spaces=None, real_power=True, real_signal=True): ...@@ -173,8 +173,8 @@ def power_synthesize(field, spaces=None, real_power=True, real_signal=True):
# power spectrum is real or complex # power spectrum is real or complex
result = [field.from_random('normal', mean=0., std=1., result = [field.from_random('normal', mean=0., std=1.,
domain=spec.domain, domain=spec.domain,
dtype=np.float if real_signal dtype=np.float64 if real_signal
else np.complex) else np.complex128)
for x in range(1 if real_power else 2)] for x in range(1 if real_power else 2)]
result[0] *= spec if self_real else spec.real result[0] *= spec if self_real else spec.real
...@@ -231,7 +231,7 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None): ...@@ -231,7 +231,7 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None):
space = 0 space = 0
space = int(space) space = int(space)
return DiagonalOperator( return DiagonalOperator(
create_power_field(domain[space], power_spectrum, dtype).weight(1), create_power_field(domain[space], power_spectrum, dtype),
domain=domain, domain=domain,
spaces=space) spaces=space)
...@@ -239,7 +239,6 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None): ...@@ -239,7 +239,6 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None):
def create_composed_fft_operator(domain, codomain=None, all_to='other'): def create_composed_fft_operator(domain, codomain=None, all_to='other'):
if codomain is None: if codomain is None:
codomain = [None]*len(domain) codomain = [None]*len(domain)
tdom = list(domain)
res = None res = None
for i, space in enumerate(domain): for i, space in enumerate(domain):
if not isinstance(space, Space): if not isinstance(space, Space):
...@@ -247,13 +246,9 @@ def create_composed_fft_operator(domain, codomain=None, all_to='other'): ...@@ -247,13 +246,9 @@ def create_composed_fft_operator(domain, codomain=None, all_to='other'):
if (all_to == 'other' or if (all_to == 'other' or
(all_to == 'position' and space.harmonic) or (all_to == 'position' and space.harmonic) or
(all_to == 'harmonic' and not space.harmonic)): (all_to == 'harmonic' and not space.harmonic)):
if codomain[i] is None: tdom = domain if res is None else res.target
tgt = domain[i].get_default_codomain() op = FFTOperator(domain=tdom, target=codomain[i], space=i)
else:
tgt = codomain[i]
op = FFTOperator(domain=tdom, target=tgt, space=i)
res = op if res is None else op*res res = op if res is None else op*res
tdom[i] = tgt
if res is None: if res is None:
raise ValueError("empty operator") raise ValueError("empty operator")
return res return res
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