Commit ab8ac790 authored by Martin Reinecke's avatar Martin Reinecke

PEP8 and various hacks to make new functionality work with MPI (this neeeds further work)

parent 64ea7dd9
...@@ -6,19 +6,19 @@ def make_chess_mask(): ...@@ -6,19 +6,19 @@ def make_chess_mask():
mask = np.ones(position_space.shape) mask = np.ones(position_space.shape)
for i in range(4): for i in range(4):
for j in range(4): for j in range(4):
if (i+j)%2 == 0: if (i+j) % 2 == 0:
mask[i*128//4:(i+1)*128//4, j*128//4:(j+1)*128//4] = 0 mask[i*128//4:(i+1)*128//4, j*128//4:(j+1)*128//4] = 0
return mask return mask
def make_random_mask(): def make_random_mask():
mask = ift.from_random('pm1',position_space) mask = ift.from_random('pm1', position_space)
mask = (mask+1)/2 mask = (mask+1)/2
return mask.val return mask.val
if __name__ == '__main__': if __name__ == '__main__':
## describtion of the tutorial ### # # description of the tutorial ###
# Choose problem geometry and masking # Choose problem geometry and masking
...@@ -34,11 +34,11 @@ if __name__ == '__main__': ...@@ -34,11 +34,11 @@ if __name__ == '__main__':
# position_space = ift.HPSpace(128) # position_space = ift.HPSpace(128)
# mask = make_random_mask() # mask = make_random_mask()
harmonic_space = position_space.get_default_codomain() harmonic_space = position_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space) HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space)
# set correlation structure with a power spectrum and build prior correlation covariance # set correlation structure with a power spectrum and build
# prior correlation covariance
def power_spectrum(k): def power_spectrum(k):
return 100. / (20.+k**3) return 100. / (20.+k**3)
power_space = ift.PowerSpace(harmonic_space) power_space = ift.PowerSpace(harmonic_space)
...@@ -47,9 +47,10 @@ if __name__ == '__main__': ...@@ -47,9 +47,10 @@ if __name__ == '__main__':
S = ift.DiagonalOperator(prior_correlation_structure) S = ift.DiagonalOperator(prior_correlation_structure)
# build instrument response consisting of a discretization, mask and harmonic transformaion # build instrument response consisting of a discretization, mask
# and harmonic transformaion
GR = ift.GeometryRemover(position_space) GR = ift.GeometryRemover(position_space)
mask = ift.Field(position_space,val=mask) mask = ift.Field.from_global_data(position_space, mask)
Mask = ift.DiagonalOperator(mask) Mask = ift.DiagonalOperator(mask)
R = GR * Mask * HT R = GR * Mask * HT
...@@ -69,10 +70,10 @@ if __name__ == '__main__': ...@@ -69,10 +70,10 @@ if __name__ == '__main__':
D_inv = R.adjoint * N.inverse * R + S.inverse D_inv = R.adjoint * N.inverse * R + S.inverse
# make it invertible # make it invertible
IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3) IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3)
D = ift.InversionEnabler(D_inv,IC,approximation=S.inverse).inverse D = ift.InversionEnabler(D_inv, IC, approximation=S.inverse).inverse
# WIENER FILTER # WIENER FILTER
m = D(j) m = D(j)
##PLOTTING # PLOTTING
#Truth, data, reconstruction, residuals # Truth, data, reconstruction, residuals
...@@ -6,14 +6,14 @@ def get_2D_exposure(): ...@@ -6,14 +6,14 @@ def get_2D_exposure():
x_shape, y_shape = position_space.shape x_shape, y_shape = position_space.shape
exposure = np.ones(position_space.shape) exposure = np.ones(position_space.shape)
exposure[x_shape/3:x_shape/2,:] *= 2. exposure[x_shape/3:x_shape/2, :] *= 2.
exposure[x_shape*4/5:x_shape,:] *= .1 exposure[x_shape*4/5:x_shape, :] *= .1
exposure[x_shape/2:x_shape*3/2,:] *=3. exposure[x_shape/2:x_shape*3/2, :] *= 3.
exposure[:,x_shape / 3:x_shape / 2] *= 2. exposure[:, x_shape/3:x_shape/2] *= 2.
exposure[:,x_shape * 4 / 5:x_shape] *= .1 exposure[:, x_shape*4/5:x_shape] *= .1
exposure[:,x_shape / 2:x_shape * 3 / 2] *= 3. exposure[:, x_shape/2:x_shape*3/2] *= 3.
exposure = ift.Field(position_space, val=exposure) exposure = ift.Field.from_global_data(position_space, exposure)
return exposure return exposure
...@@ -27,14 +27,13 @@ if __name__ == '__main__': ...@@ -27,14 +27,13 @@ if __name__ == '__main__':
# position_space = ift.RGSpace(1024) # position_space = ift.RGSpace(1024)
# exposure = np.ones(position_space.shape) # exposure = np.ones(position_space.shape)
# Two-dimensional regular grid with inhomogeneous exposure
# Two dimensional regular grid with inhomogeneous exposure
position_space = ift.RGSpace([512, 512]) position_space = ift.RGSpace([512, 512])
exposure = get_2D_exposure() exposure = get_2D_exposure()
# # Sphere with with uniform exposure # # Sphere with with uniform exposure
# position_space = ift.HPSpace(128) # position_space = ift.HPSpace(128)
# exposure = np.ones(position_space.shape) # exposure = ift.Field.full(position_space, 1.)
# Defining harmonic space and transform # Defining harmonic space and transform
harmonic_space = position_space.get_default_codomain() harmonic_space = position_space.get_default_codomain()
...@@ -58,7 +57,6 @@ if __name__ == '__main__': ...@@ -58,7 +57,6 @@ if __name__ == '__main__':
logsky = HT(logsky_h) logsky = HT(logsky_h)
sky = ift.PointwiseExponential(logsky) sky = ift.PointwiseExponential(logsky)
exposure = ift.Field(position_space, val=exposure)
M = ift.DiagonalOperator(exposure) M = ift.DiagonalOperator(exposure)
GR = ift.GeometryRemover(position_space) GR = ift.GeometryRemover(position_space)
# Set up instrumental response # Set up instrumental response
...@@ -68,14 +66,16 @@ if __name__ == '__main__': ...@@ -68,14 +66,16 @@ if __name__ == '__main__':
d_space = R.target[0] d_space = R.target[0]
lamb = R(sky) lamb = R(sky)
mock_position = ift.from_random('normal', lamb.position.domain) mock_position = ift.from_random('normal', lamb.position.domain)
data = np.random.poisson(lamb.at(mock_position).value.val.astype(np.float64)) data = lamb.at(mock_position).value
data = ift.Field.from_local_data(d_space, data) data = np.random.poisson(data.to_global_data().astype(np.float64))
data = ift.Field.from_global_data(d_space, data)
# Compute likelihood and Hamiltonian # Compute likelihood and Hamiltonian
position = ift.from_random('normal', lamb.position.domain) position = ift.from_random('normal', lamb.position.domain)
likelihood = ift.PoissonianEnergy(lamb, data) likelihood = ift.PoissonianEnergy(lamb, data)
ic_cg = ift.GradientNormController(iteration_limit=50) ic_cg = ift.GradientNormController(iteration_limit=50)
ic_newton = ift.GradientNormController(name='Newton',iteration_limit=50, tol_abs_gradnorm=1e-3) ic_newton = ift.GradientNormController(name='Newton', iteration_limit=50,
tol_abs_gradnorm=1e-3)
minimizer = ift.RelaxedNewton(ic_newton) minimizer = ift.RelaxedNewton(ic_newton)
# Minimize the Hamiltonian # Minimize the Hamiltonian
...@@ -84,6 +84,4 @@ if __name__ == '__main__': ...@@ -84,6 +84,4 @@ if __name__ == '__main__':
# Plot results # Plot results
result_sky = sky.at(H.position).value result_sky = sky.at(H.position).value
##PLOTTING # PLOTTING
...@@ -7,18 +7,20 @@ from scipy.io import loadmat ...@@ -7,18 +7,20 @@ from scipy.io import loadmat
def get_random_LOS(n_los): def get_random_LOS(n_los):
starts = list(np.random.uniform(0,1,(n_los,2)).T) starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
ends = list(np.random.uniform(0,1,(n_los,2)).T) ends = list(np.random.uniform(0, 1, (n_los, 2)).T)
return starts, ends return starts, ends
if __name__ == '__main__': if __name__ == '__main__':
### ABOUT THIS TUTORIAL # ## ABOUT THIS TUTORIAL
np.random.seed(42) np.random.seed(42)
position_space = ift.RGSpace([128,128]) position_space = ift.RGSpace([128, 128])
# Setting up an amplitude model # Setting up an amplitude model
A, amplitude_internals = make_amplitude_model(position_space,16, 1, 10, -4., 1, 0., 1.) A, amplitude_internals = make_amplitude_model(
position_space, 16, 1, 10, -4., 1, 0., 1.)
# Building the model for a correlated signal # Building the model for a correlated signal
harmonic_space = position_space.get_default_codomain() harmonic_space = position_space.get_default_codomain()
...@@ -46,8 +48,7 @@ if __name__ == '__main__': ...@@ -46,8 +48,7 @@ if __name__ == '__main__':
# specify noise # specify noise
data_space = R.target data_space = R.target
noise = .001 noise = .001
N = ift.ScalingOperator(noise,data_space) N = ift.ScalingOperator(noise, data_space)
# generate mock data # generate mock data
MOCK_POSITION = ift.from_random('normal', signal.position.domain) MOCK_POSITION = ift.from_random('normal', signal.position.domain)
...@@ -63,28 +64,31 @@ if __name__ == '__main__': ...@@ -63,28 +64,31 @@ if __name__ == '__main__':
minimizer = ift.RelaxedNewton(ic_newton) minimizer = ift.RelaxedNewton(ic_newton)
# build model Hamiltonian # build model Hamiltonian
H = ift.Hamiltonian(likelihood,ic_cg,iteration_controller_sampling=ic_sampling) H = ift.Hamiltonian(likelihood, ic_cg,
iteration_controller_sampling=ic_sampling)
INITIAL_POSITION = ift.from_random('normal',H.position.domain) INITIAL_POSITION = ift.from_random('normal', H.position.domain)
position = INITIAL_POSITION position = INITIAL_POSITION
ift.plot(signal.at(MOCK_POSITION).value,name='truth.pdf') ift.plot(signal.at(MOCK_POSITION).value, name='truth.pdf')
ift.plot(R.adjoint_times(data),name='data.pdf') ift.plot(R.adjoint_times(data), name='data.pdf')
ift.plot([ A.at(MOCK_POSITION).value], name='power.pdf') ift.plot([A.at(MOCK_POSITION).value], name='power.pdf')
# number of samples used to estimate the KL # number of samples used to estimate the KL
N_samples = 20 N_samples = 20
for i in range(5): for i in range(5):
H = H.at(position) H = H.at(position)
samples = [H.curvature.draw_sample(from_inverse=True) for _ in range(N_samples)] samples = [H.curvature.draw_sample(from_inverse=True)
for _ in range(N_samples)]
KL = ift.SampledKullbachLeiblerDivergence(H, samples, ic_cg) KL = ift.SampledKullbachLeiblerDivergence(H, samples, ic_cg)
KL, convergence = minimizer(KL) KL, convergence = minimizer(KL)
position = KL.position position = KL.position
ift.plot(signal.at(position).value,name='reconstruction.pdf') ift.plot(signal.at(position).value, name='reconstruction.pdf')
ift.plot([A.at(position).value, A.at(MOCK_POSITION).value],name='power.pdf') ift.plot([A.at(position).value, A.at(MOCK_POSITION).value],
name='power.pdf')
avrg = 0. avrg = 0.
va = 0. va = 0.
...@@ -101,15 +105,5 @@ if __name__ == '__main__': ...@@ -101,15 +105,5 @@ if __name__ == '__main__':
std = ift.sqrt(va) std = ift.sqrt(va)
ift.plot(avrg, name='avrg.pdf') ift.plot(avrg, name='avrg.pdf')
ift.plot(std, name='std.pdf') ift.plot(std, name='std.pdf')
ift.plot([A.at(position).value, A.at(MOCK_POSITION).value]+powers, name='power.pdf') ift.plot([A.at(position).value, A.at(MOCK_POSITION).value]+powers,
name='power.pdf')
...@@ -22,7 +22,7 @@ def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv, ...@@ -22,7 +22,7 @@ def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
Npixdof : #pix in dof_space Npixdof : #pix in dof_space
ceps_a, ceps_k0 : Smoothnessparameters in ceps_kernel ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel
eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
a = ceps_a, k0 = ceps_k0 a = ceps_a, k0 = ceps_k0
...@@ -48,7 +48,7 @@ def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv, ...@@ -48,7 +48,7 @@ def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
phi_sig = np.array([sv, iv]) phi_sig = np.array([sv, iv])
slope = SlopeOperator(param_space, logk_space, phi_sig) slope = SlopeOperator(param_space, logk_space, phi_sig)
norm_phi_mean = Field(param_space, val=phi_mean/phi_sig) norm_phi_mean = Field.from_global_data(param_space, phi_mean/phi_sig)
fields = {keys[0]: Field.from_random('normal', dof_space), fields = {keys[0]: Field.from_random('normal', dof_space),
keys[1]: Field.from_random('normal', param_space)} keys[1]: Field.from_random('normal', param_space)}
...@@ -93,7 +93,7 @@ def create_cepstrum_amplitude_field(domain, cepstrum): ...@@ -93,7 +93,7 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
# Prepare q_array # Prepare q_array
q_array = np.zeros((dim,) + shape) q_array = np.zeros((dim,) + shape)
if dim == 1: if dim == 1:
ks = domain.get_k_length_array().val ks = domain.get_k_length_array().to_global_data()
q_array = np.array([ks]) q_array = np.array([ks])
else: else:
for i in range(dim): for i in range(dim):
...@@ -120,4 +120,4 @@ def create_cepstrum_amplitude_field(domain, cepstrum): ...@@ -120,4 +120,4 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
# Do summation # Do summation
cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i) cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
return Field(domain, val=cepstrum_field) return Field.from_global_data(domain, cepstrum_field)
...@@ -67,7 +67,7 @@ class MultiField(object): ...@@ -67,7 +67,7 @@ class MultiField(object):
dtype = MultiField.build_dtype(dtype, domain) dtype = MultiField.build_dtype(dtype, domain)
return MultiField({key: Field.from_random(random_type, domain[key], return MultiField({key: Field.from_random(random_type, domain[key],
dtype[key], **kwargs) dtype[key], **kwargs)
for key in domain.keys()}) for key in sorted(domain.keys())})
def fill(self, fill_value): def fill(self, fill_value):
"""Fill `self` uniformly with `fill_value` """Fill `self` uniformly with `fill_value`
......
...@@ -67,7 +67,7 @@ class ExpTransform(LinearOperator): ...@@ -67,7 +67,7 @@ class ExpTransform(LinearOperator):
def apply(self, x, mode): def apply(self, x, mode):
self._check_input(x, mode) self._check_input(x, mode)
x = x.val x = x.to_global_data()
ndim = len(self.target.shape) ndim = len(self.target.shape)
idx = () idx = ()
for d in range(ndim): for d in range(ndim):
...@@ -90,7 +90,7 @@ class ExpTransform(LinearOperator): ...@@ -90,7 +90,7 @@ class ExpTransform(LinearOperator):
x = xnew x = xnew
idx = (slice(None),) + idx idx = (slice(None),) + idx
return Field(self._tgt(mode), val=x) return Field.from_global_data(self._tgt(mode), x)
@property @property
def capability(self): def capability(self):
......
from ..domain_tuple import DomainTuple from ..domain_tuple import DomainTuple
from ..field import Field from ..field import Field
from .. import dobj
from ..utilities import hartley from ..utilities import hartley
from .linear_operator import LinearOperator from .linear_operator import LinearOperator
...@@ -34,11 +35,13 @@ class QHTOperator(LinearOperator): ...@@ -34,11 +35,13 @@ class QHTOperator(LinearOperator):
x = x.val * self.domain[0].scalar_dvol() x = x.val * self.domain[0].scalar_dvol()
n = len(self.domain[0].shape) n = len(self.domain[0].shape)
rng = range(n) if mode == self.TIMES else reversed(range(n)) rng = range(n) if mode == self.TIMES else reversed(range(n))
# MR FIXME: this needs to be fixed properly for MPI
x = dobj.to_global_data(x)
for i in rng: for i in rng:
sl = (slice(None),)*i + (slice(1, None),) sl = (slice(None),)*i + (slice(1, None),)
x[sl] = hartley(x[sl], axes=(i,)) x[sl] = hartley(x[sl], axes=(i,))
return Field(self._tgt(mode), val=x) return Field.from_global_data(self._tgt(mode), x)
@property @property
def capability(self): def capability(self):
......
...@@ -18,7 +18,7 @@ class SlopeOperator(LinearOperator): ...@@ -18,7 +18,7 @@ class SlopeOperator(LinearOperator):
self.pos = np.zeros((self.ndim,) + self.target[0].shape) self.pos = np.zeros((self.ndim,) + self.target[0].shape)
if self.ndim == 1: if self.ndim == 1:
self.pos[0] = self.target[0].get_k_length_array().val self.pos[0] = self.target[0].get_k_length_array().to_global_data()
else: else:
shape = self.target[0].shape shape = self.target[0].shape
for i in range(self.ndim): for i in range(self.ndim):
...@@ -46,18 +46,19 @@ class SlopeOperator(LinearOperator): ...@@ -46,18 +46,19 @@ class SlopeOperator(LinearOperator):
# Times # Times
if mode == self.TIMES: if mode == self.TIMES:
inp = x.val inp = x.to_global_data()
res = self.sigmas[-1] * inp[-1] res = self.sigmas[-1] * inp[-1]
for i in range(self.ndim): for i in range(self.ndim):
res += self.sigmas[i] * inp[i] * self.pos[i] res += self.sigmas[i] * inp[i] * self.pos[i]
return Field(self.target, val=res) return Field.from_global_data(self.target, res)
# Adjoint times # Adjoint times
res = np.zeros(self.domain[0].shape) res = np.zeros(self.domain[0].shape)
res[-1] = np.sum(x.val) * self.sigmas[-1] xglob = x.to_global_data()
res[-1] = np.sum(xglob) * self.sigmas[-1]
for i in range(self.ndim): for i in range(self.ndim):
res[i] = np.sum(self.pos[i] * x.val) * self.sigmas[i] res[i] = np.sum(self.pos[i] * xglob) * self.sigmas[i]
return Field(self.domain, val=res) return Field.from_global_data(self.domain, res)
@property @property
def capability(self): def capability(self):
......
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