From e240a3e2d724f694752d23da0fee3d9208ae00bc Mon Sep 17 00:00:00 2001
From: caron <j.caron@fz-juelich.de>
Date: Tue, 12 Dec 2017 17:00:02 +0100
Subject: [PATCH] forwardmodel: DistributedForwardModel now correctly
distributes images! also: now works with logging instead of printing!
projector: Fixed bug with shape of impacts (was erroneously transposed)!
regularisator: Regularisator base class now in __all__!
---
pyramid/dataset.py | 11 ++++----
pyramid/forwardmodel.py | 54 ++++++++++++++++++++++------------------
pyramid/kernel.py | 2 +-
pyramid/projector.py | 22 ++++++++--------
pyramid/regularisator.py | 4 +--
5 files changed, 50 insertions(+), 43 deletions(-)
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index cae8f64..3e6a7a3 100644
--- a/pyramid/dataset.py
+++ b/pyramid/dataset.py
@@ -142,7 +142,7 @@ class DataSet(object):
self.dim = dim
self.b_0 = b_0 # TODO: Not very general!!! get rid off!
self.mask = mask
- self.Se_inv = Se_inv
+ self.Se_inv = Se_inv # TODO: Should be constructed in the costfunction, not here!!
self._phasemaps = []
self._projectors = []
self._phasemappers = []
@@ -219,7 +219,7 @@ class DataSet(object):
for i in range(len(phasemap)):
self._append_single(phasemap[i], projector[i], phasemapper[i])
# Reset the Se_inv matrix from phasemaps confidence matrices:
- self.set_Se_inv_diag_with_conf()
+ self.set_Se_inv_diag_with_conf() # TODO: Very clunky and redundant when adding a lot!
def create_phasemaps(self, magdata, difference=False, ramp=None):
"""Create a list of phasemaps with the projectors in the dataset for a given
@@ -560,12 +560,13 @@ class DataSetCharge(object):
def __repr__(self):
self._log.debug('Calling __repr__')
- return '%s(a=%r, dim=%r, electrode_vec=%r, mask=%r, Se_inv=%r)' % (self.__class__, self.a, self.dim,
- self.electrode_vec, self.mask, self.Se_inv)
+ return '%s(a=%r, dim=%r, electrode_vec=%r, mask=%r, Se_inv=%r)' % \
+ (self.__class__, self.a, self.dim, self.electrode_vec, self.mask, self.Se_inv)
def __str__(self):
self._log.debug('Calling __str__')
- return 'DataSetCharge(a=%s, dim=%s, electrode_vec=%s)' % (self.a, self.dim, self.electrode_vec)
+ return 'DataSetCharge(a=%s, dim=%s, electrode_vec=%s)' % \
+ (self.a, self.dim, self.electrode_vec)
def _append_single(self, phasemap, projector, phasemapper=None):
self._log.debug('Calling _append')
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index 90f08b3..806465e 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -365,28 +365,41 @@ class DistributedForwardModel(ForwardModel):
All ramp parameters have to be at the end of the input vector and are split automatically.
Default is None (no ramps are added).
nprocs: int
- Number of processes which should be created. Default is 1 (not recommended).
+ Number of processes which should be created. Default is 1 (not recommended). # TODO: <<<!!!
"""
- def __init__(self, data_set, ramp_order=None, nprocs=1):
+ _log = mp.log_to_stderr()
+
+ def __init__(self, data_set, ramp_order=None, nprocs='auto'):
# Evoke super constructor to set up the normal ForwardModel:
super().__init__(data_set, ramp_order)
# Initialize multirocessing specific stuff:
+ if nprocs == 'auto':
+ nprocs = mp.cpu_count() - 2 # Use two cores less to reserve cpu for the system.
self.nprocs = nprocs
- img_per_proc = np.ceil(self.data_set.count / self.nprocs).astype(np.int)
- hp = self.data_set.hook_points
- self.proc_hook_points = [0]
self.pipes = []
self.processes = []
- print('NPROCS:', self.nprocs) # TODO: Logging instead of printing!
+ self.proc_hook_points = [0] # Hook points of the processes in the output vector
+ hp = self.data_set.hook_points # Hook points of the images in the output vector
+ # Calculate the best distribution of images to the processes:
+ proc_img_range = []
+ n = self.data_set.count // nprocs # min items per process
+ r = self.data_set.count % nprocs # remainder items
+ start, stop = 0, n + (r != 0)
for proc_id in range(self.nprocs):
+ proc_img_range.append((start, stop))
+ if r > 0: # use up remainder:
+ r -= 1
+ start = stop
+ stop = stop + n + (r != 0)
+ # Set up the workers:
+ self._log.info('Creating {} processes'.format(self.nprocs))
+ for proc_id, (start, stop) in enumerate(proc_img_range):
# Create SubDataSets:
sub_data = DataSet(self.data_set.a, self.data_set.dim, self.data_set.b_0,
- self.data_set.mask, self.data_set.Se_inv)
+ self.data_set.mask, Se_inv=None) # Se_inv is set later!
# Distribute data to SubDataSets:
- start = proc_id * img_per_proc
- stop = np.min(((proc_id + 1) * img_per_proc, self.data_set.count))
self.proc_hook_points.append(hp[stop])
phasemaps = self.data_set.phasemaps[start:stop]
projectors = self.data_set.projectors[start:stop]
@@ -396,7 +409,7 @@ class DistributedForwardModel(ForwardModel):
# Create communication pipe:
self.pipes.append(mp.Pipe())
# Create process:
- p = mp.Process(name='Worker {}'.format(proc_id), target=_worker,
+ p = mp.Process(name='Worker {}'.format(proc_id), target=self._worker,
args=(sub_fwd_model, self.pipes[proc_id][1]))
self.processes.append(p)
# Start process:
@@ -423,6 +436,13 @@ class DistributedForwardModel(ForwardModel):
# Return result:
return result
+ def _worker(self, fwd_model, pipe):
+ for method, arguments in iter(pipe.recv, 'STOP'):
+ sys.stdout.flush()
+ result = getattr(fwd_model, method)(*arguments)
+ pipe.send(result)
+ sys.stdout.flush()
+
def jac_dot(self, x, vector):
"""Calculate the product of the Jacobi matrix with a given `vector`.
@@ -511,17 +531,3 @@ class DistributedForwardModel(ForwardModel):
# Exit the completed processes:
for p in self.processes:
p.join()
-
-
-def _worker(fwd_model, pipe):
- # Has to be directly accessible in the module as a function, NOT a method of a class instance!
- # TODO: Logging instead of printing!
- print('... {} starting!'.format(mp.current_process().name))
- sys.stdout.flush()
- for method, arguments in iter(pipe.recv, 'STOP'):
- # '... {} processes method {}'.format(mp.current_process().name, method)
- sys.stdout.flush()
- result = getattr(fwd_model, method)(*arguments)
- pipe.send(result)
- print('... ', mp.current_process().name, 'exiting!')
- sys.stdout.flush()
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index 1e1d149..c78aa4d 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -240,7 +240,7 @@ class KernelCharge(object):
self.slice_c = (slice(0, dim_uv[0]), # Charge is padded on the far end!
slice(0, dim_uv[1])) # (Phase cutout is shifted as listed above)
# Calculate kernel (single pixel phase):
- coeff = a * C_e * Q_E / (4 * np.pi * EPS_0) # Minus is gone because of negative z-direction
+ coeff = a * C_e * Q_E / (4 * np.pi * EPS_0) # Minus gone because of negative z-direction
v_dim, u_dim = dim_uv
u = np.linspace(-(u_dim - 1), u_dim - 1, num=2 * u_dim - 1)
v = np.linspace(-(v_dim - 1), v_dim - 1, num=2 * v_dim - 1)
diff --git a/pyramid/projector.py b/pyramid/projector.py
index 99ff190..9830be2 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -427,17 +427,17 @@ class RotTiltProjector(Projector):
dim_v, dim_u = dim_uv
# Creating coordinate list of all voxels:
voxels = list(itertools.product(range(dim_z), range(dim_y), range(dim_x)))
- # Calculate vectors to voxels relative to rotation center (each row contains (z, y, x)):
- voxel_vecs = np.asarray(voxels) + 0.5 - np.asarray(center)
- # Change to coordinate order of quaternions (x, y, z) instead of (z, y, x):
- voxel_vecs = np.fliplr(voxel_vecs)
+ # Calculate vectors to voxels relative to rotation center (each column contains (z, y, x)):
+ voxel_vecs = (np.asarray(voxels) + 0.5 - np.asarray(center)).T # .T: row to column!
+ # Change to coordinate order of quaternions (x, y, z) instead of (z, y, x) per column:
+ voxel_vecs = np.flipud(voxel_vecs)
# Calculate impact positions (x, y) on the projected pixel coordinate grid (z is dropped):
- impacts = quat.matrix[:2, :].dot(voxel_vecs.T).T # Only x and y row of matrix is used!
- # Reverse transpose and change back coordinate order from (x, y) to (y, x)/(v, u):
- impacts = np.fliplr(impacts.T) # Now contains rows with (v, u) entries as desired!
- # First index: voxel, second index: 0 -> v, 1 -> u!
- impacts[:, 1] += dim_u / 2. # Shift back to normal indices
- impacts[:, 0] += dim_v / 2. # Shift back to normal indices
+ impacts = quat.matrix[:2, :].dot(voxel_vecs) # Only x and y row of matrix is used!
+ # Change back coordinate order from (x, y) to (y, x)/(v, u) per column:
+ impacts = np.flipud(impacts) # Now contains columns with (v, u) entries (pythonic index)!
+ # First index (column): 0 -> v, 1 -> u, second index (row): voxel!
+ impacts[0, :] += dim_v / 2. # Shift back to normal indices
+ impacts[1, :] += dim_u / 2. # Shift back to normal indices
# Calculate equivalence radius:
R = (3 / (4 * np.pi)) ** (1 / 3.)
# Prepare weight matrix calculation:
@@ -451,7 +451,7 @@ class RotTiltProjector(Projector):
for i, voxel in enumerate(tqdm(voxels, disable=disable, leave=False,
desc='Set up projector')):
column_index = voxel[0] * dim_y * dim_x + voxel[1] * dim_x + voxel[2]
- remainder, impact = np.modf(impacts[i, :]) # split index of impact and remainder!
+ remainder, impact = np.modf(impacts[:, i]) # split index of impact and remainder!
sub_pixel = (remainder * subcount).astype(dtype=np.int) # sub_pixel inside impact px.
# Go over all influenced pixels (impact and neighbours, indices are [0, 1, 2]!):
for px_ind in list(itertools.product(range(3), range(3))):
diff --git a/pyramid/regularisator.py b/pyramid/regularisator.py
index 8088ccd..4c4caf2 100644
--- a/pyramid/regularisator.py
+++ b/pyramid/regularisator.py
@@ -14,8 +14,8 @@ from scipy import sparse
import jutil.diff as jdiff
import jutil.norms as jnorm
-__all__ = ['NoneRegularisator', 'ZeroOrderRegularisator', 'FirstOrderRegularisator',
- 'ComboRegularisator']
+__all__ = ['Regularisator', 'NoneRegularisator', 'ZeroOrderRegularisator',
+ 'FirstOrderRegularisator', 'ComboRegularisator']
class Regularisator(object):
--
GitLab