Commit 7759d1f4 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

fixes for MPI; better type checks

parent 41bb082f
Pipeline #22460 passed with stage
in 4 minutes and 45 seconds
......@@ -66,19 +66,18 @@ if __name__ == "__main__":
# m3 = fft(me3[0].position)
# Plotting
ift.plotting.plot(mock_signal.real, name='mock_signal.png',
colormap="plasma", xlabel="Pixel Index",
ylabel="Pixel Index")
logdata = np.log(ift.dobj.to_global_data(data.val.real)).reshape(signal_space.shape)
ift.plotting.plot(mock_signal, name='mock_signal.png', colormap="plasma",
xlabel="Pixel Index", ylabel="Pixel Index")
logdata = np.log(ift.dobj.to_global_data(data.val)).reshape(signal_space.shape)
ift.plotting.plot(ift.Field(signal_space,
val=ift.dobj.from_global_data(logdata)),
name="log_of_data.png", colormap="plasma",
xlabel="Pixel Index", ylabel="Pixel Index")
# ift.plotting.plot(m1.real,name='m_LBFGS.png', colormap="plasma",
# ift.plotting.plot(m1,name='m_LBFGS.png', colormap="plasma",
# xlabel="Pixel Index", ylabel="Pixel Index")
ift.plotting.plot(m2.real, name='m_Newton.png', colormap="plasma",
ift.plotting.plot(m2, name='m_Newton.png', colormap="plasma",
xlabel="Pixel Index", ylabel="Pixel Index")
# ift.plotting.plot(m3.real, name='m_SteepestDescent.png',
# ift.plotting.plot(m3, name='m_SteepestDescent.png',
# colormap="plasma", xlabel="Pixel Index",
# ylabel="Pixel Index")
......
......@@ -46,7 +46,7 @@ if __name__ == "__main__":
# Instrument = SmoothingOperator(s_space, sigma=0.01)
mask = np.ones(s_space.shape)
mask[6000:8000] = 0.
mask = ift.Field(s_space, val=mask)
mask = ift.Field(s_space, val=ift.dobj.from_global_data(mask))
MaskOperator = ift.DiagonalOperator(mask)
InstrumentResponse = ift.ResponseOperator(s_space, sigma=[0.0])
......
......@@ -53,7 +53,6 @@ if __name__ == "__main__":
ift.dobj.to_global_data(mock_power_2.val))))
diagonal = ift.power_synthesize_special(mock_power, spaces=(0, 1))**2
diagonal = diagonal.real
S = ift.DiagonalOperator(diagonal)
......
......@@ -64,7 +64,7 @@ if __name__ == "__main__":
Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)
sp = ift.Field(p_space, val=pow_spec(p_space.k_lengths))
sp = ift.PS_field(p_space, pow_spec)
sh = ift.power_synthesize(sp, real_signal=True)
ss = fft.inverse_times(sh)
......
......@@ -47,7 +47,6 @@ if __name__ == "__main__":
mock_power = ift.PS_field(power_space, power_spectrum)
mock_harmonic = ift.power_synthesize(mock_power, real_signal=True)
mock_harmonic = mock_harmonic.real
mock_signal = fft(mock_harmonic)
exposure = 1.
......@@ -78,9 +77,9 @@ if __name__ == "__main__":
sspace2 = ift.RGSpace(shape, distances=L/N_pixels/nu.m)
ift.plotting.plot(ift.Field(sspace2, mock_signal.real.val)/nu.K,
ift.plotting.plot(ift.Field(sspace2, mock_signal.val)/nu.K,
name="mock_signal.png")
data = ift.dobj.to_global_data(data.val.real).reshape(sspace2.shape)/nu.K
data = ift.dobj.to_global_data(data.val).reshape(sspace2.shape)/nu.K
data = ift.Field(sspace2, val=ift.dobj.from_global_data(data))/nu.K
ift.plotting.plot(ift.Field(sspace2, val=data), name="data.png")
ift.plotting.plot(ift.Field(sspace2, m_s.real.val)/nu.K, name="map.png")
ift.plotting.plot(ift.Field(sspace2, m_s.val)/nu.K, name="map.png")
......@@ -23,7 +23,7 @@ import numpy as np
class Random(object):
@staticmethod
def pm1(dtype, shape):
if issubclass(dtype, (complex, np.complexfloating)):
if np.issubdtype(dtype, np.complexfloating):
x = np.array([1+0j, 0+1j, -1+0j, 0-1j], dtype=dtype)
x = x[np.random.randint(4, size=shape)]
else:
......@@ -32,7 +32,7 @@ class Random(object):
@staticmethod
def normal(dtype, shape, mean=0., std=1.):
if issubclass(dtype, (complex, np.complexfloating)):
if np.issubdtype(dtype, np.complexfloating):
x = np.empty(shape, dtype=dtype)
x.real = np.random.normal(mean.real, std*np.sqrt(0.5), shape)
x.imag = np.random.normal(mean.imag, std*np.sqrt(0.5), shape)
......@@ -42,12 +42,11 @@ class Random(object):
@staticmethod
def uniform(dtype, shape, low=0., high=1.):
if issubclass(dtype, (complex, np.complexfloating)):
if np.issubdtype(dtype, np.complexfloating):
x = np.empty(shape, dtype=dtype)
x.real = np.random.uniform(low, high, shape)
x.imag = np.random.uniform(low, high, shape)
elif dtype in [np.dtype('int8'), np.dtype('int16'), np.dtype('int32'),
np.dtype('int64')]:
elif np.issubdtype(dtype, np.integer):
x = np.random.random.randint(low, high+1, shape)
else:
x = np.random.uniform(low, high, shape)
......
......@@ -213,11 +213,15 @@ class Field(object):
@property
def real(self):
""" The real part of the field (data is not copied)."""
if not np.issubdtype(self.dtype, np.complexfloating):
raise ValueError(".real called on a non-complex Field")
return Field(self.domain, self.val.real)
@property
def imag(self):
""" The imaginary part of the field (data is not copied)."""
if not np.issubdtype(self.dtype, np.complexfloating):
raise ValueError(".imag called on a non-complex Field")
return Field(self.domain, self.val.imag)
def copy(self):
......
......@@ -89,7 +89,7 @@ class CriticalPowerEnergy(Energy):
gradient = -self._theta
gradient += self.alpha-0.5
gradient += Tt
self._gradient = gradient.real
self._gradient = gradient
def at(self, position):
return self.__class__(position, self.m, D=self.D, alpha=self.alpha,
......
......@@ -140,7 +140,7 @@ class DiagonalOperator(EndomorphicOperator):
@property
def self_adjoint(self):
if self._self_adjoint is None:
if not issubclass(self._diagonal.dtype.type, np.complexfloating):
if not np.issubdtype(self._diagonal.dtype, np.complexfloating):
self._self_adjoint = True
else:
self._self_adjoint = (self._diagonal.val.imag == 0).all()
......
......@@ -107,7 +107,7 @@ class FFTOperator(LinearOperator):
self._trafo = SphericalTransformation(pdom, hdom, space)
def _times_helper(self, x):
if issubclass(x.dtype.type, np.complexfloating):
if np.issubdtype(x.dtype, np.complexfloating):
res = (self._trafo.transform(x.real) +
1j * self._trafo.transform(x.imag))
else:
......
......@@ -45,9 +45,11 @@ def _find_closest(A, target):
def _makeplot(name):
import matplotlib.pyplot as plt
if dobj.rank != 0:
plt.close()
return
if name is None:
plt.show()
plt.close()
return
extension = os.path.splitext(name)[1]
if extension == ".pdf":
......
......@@ -94,10 +94,17 @@ def power_analyze(field, spaces=None, binbounds=None,
if len(spaces) == 0:
raise ValueError("No space for analysis specified.")
field_real = not np.issubdtype(field.dtype, np.complexfloating)
if (not field_real) and keep_phase_information:
raise ValueError("cannot keep phase from real-valued input Field")
if keep_phase_information:
parts = [field.real*field.real, field.imag*field.imag]
else:
parts = [field.real*field.real + field.imag*field.imag]
if field_real:
parts = [field**2]
else:
parts = [field.real*field.real + field.imag*field.imag]
parts = [part.weight(1, spaces) for part in parts]
for space_index in spaces:
......@@ -162,6 +169,9 @@ def power_synthesize(field, spaces=None, real_power=True, real_signal=True):
"""
spec = _compute_spec(field, spaces)
self_real = not np.issubdtype(spec.dtype, np.complexfloating)
if (not real_power) and self_real:
raise ValueError("can't draw complex realizations from real spectrum")
# create random samples: one or two, depending on whether the
# power spectrum is real or complex
......@@ -176,7 +186,7 @@ def power_synthesize(field, spaces=None, real_power=True, real_signal=True):
field.from_random('normal', mean=0., std=1.,
domain=spec.domain, dtype=np.float)
result[0] *= spec.real
result[0] *= spec if self_real else spec.real
if not real_power:
result[1] *= spec.imag
......@@ -190,7 +200,7 @@ def power_synthesize_special(field, spaces=None):
field.from_random('normal', mean=0., std=1.,
domain=spec.domain, dtype=np.complex)
return spec.real
return spec
def create_power_field(domain, power_spectrum, dtype=None):
......@@ -206,12 +216,8 @@ def create_power_field(domain, power_spectrum, dtype=None):
else:
power_domain = PowerSpace(domain)
fp = PS_field(power_domain, power_spectrum, dtype)
f = PowerProjectionOperator(domain, power_domain).adjoint_times(fp)
if not issubclass(fp.dtype.type, np.complexfloating):
f = f.real
return f
return PowerProjectionOperator(domain, power_domain).adjoint_times(fp)
def create_power_operator(domain, power_spectrum, space=None, dtype=None):
......
......@@ -129,7 +129,7 @@ def hartley(a, axes=None):
if axes is not None and \
not all(axis < len(a.shape) for axis in axes):
raise ValueError("Provided axes do not match array shape")
if issubclass(a.dtype.type, np.complexfloating):
if np.issubdtype(a.dtype, np.complexfloating):
raise TypeError("Hartley transform requires real-valued arrays.")
from pyfftw.interfaces.numpy_fft import rfftn
......@@ -171,7 +171,7 @@ def my_fftn_r2c(a, axes=None):
if axes is not None and \
not all(axis < len(a.shape) for axis in axes):
raise ValueError("Provided axes do not match array shape")
if issubclass(a.dtype.type, np.complexfloating):
if np.issubdtype(a.dtype, np.complexfloating):
raise TypeError("Transform requires real-valued input arrays.")
from pyfftw.interfaces.numpy_fft import rfftn
......
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