Skip to content
Snippets Groups Projects
Commit fbbbf8d6 authored by Jan Caron's avatar Jan Caron
Browse files

Fixed m and n controversy:

        m: i'm'age
        n: i'n'put
parent 76c17670
No related branches found
No related tags found
No related merge requests found
...@@ -40,6 +40,10 @@ class Costfunction(object): ...@@ -40,6 +40,10 @@ class Costfunction(object):
Se_inv : :class:`~numpy.ndarray` (N=2), optional Se_inv : :class:`~numpy.ndarray` (N=2), optional
Inverted covariance matrix of the measurement errors. The matrix has size `NxN` with N Inverted covariance matrix of the measurement errors. The matrix has size `NxN` with N
being the length of the targetvector y (vectorized phase map information). being the length of the targetvector y (vectorized phase map information).
m: int
Size of the image space.
n: int
Size of the input space.
''' '''
...@@ -53,8 +57,8 @@ class Costfunction(object): ...@@ -53,8 +57,8 @@ class Costfunction(object):
# Extract important information: # Extract important information:
self.y = data_set.phase_vec self.y = data_set.phase_vec
self.Se_inv = data_set.Se_inv self.Se_inv = data_set.Se_inv
self.n = data_set.n * 3 self.n = data_set.n
self.m = len(self.y) self.m = data_set.m
self.LOG.debug('Created '+str(self)) self.LOG.debug('Created '+str(self))
def __repr__(self): def __repr__(self):
...@@ -161,7 +165,7 @@ class CFAdapterScipyCG(LinearOperator): ...@@ -161,7 +165,7 @@ class CFAdapterScipyCG(LinearOperator):
@property @property
def shape(self): def shape(self):
return (self.cost.data_set.m, self.cost.data_set.m) return (self.cost.data_set.n, self.cost.data_set.n)
@property @property
def dtype(self): def dtype(self):
......
...@@ -46,22 +46,23 @@ class DataSet(object): ...@@ -46,22 +46,23 @@ class DataSet(object):
A list of all stored :class:`~.PhaseMap` objects. A list of all stored :class:`~.PhaseMap` objects.
phase_vec: :class:`~numpy.ndarray` (N=1) phase_vec: :class:`~numpy.ndarray` (N=1)
The concatenaded, vectorized phase of all ;class:`~.PhaseMap` objects. The concatenaded, vectorized phase of all ;class:`~.PhaseMap` objects.
n: int
Size of the image space.
m: int m: int
Size of the image space.
n: int
Size of the input space. Size of the input space.
''' '''
LOG = logging.getLogger(__name__+'.DataSet') LOG = logging.getLogger(__name__+'.DataSet')
@property @property
def n(self): def m(self):
return np.sum([len(p.phase_vec) for p in self.phase_maps]) return np.sum([len(p.phase_vec) for p in self.phase_maps])
@property @property
def Se_inv(self): def Se_inv(self):
# TODO: better implementation, maybe get-method? more flexible? input in append? # TODO: better implementation, maybe get-method? more flexible? input in append?
return sp.eye(self.n) return sp.eye(self.m)
@property @property
def phase_vec(self): def phase_vec(self):
...@@ -88,9 +89,9 @@ class DataSet(object): ...@@ -88,9 +89,9 @@ class DataSet(object):
'Dimension has to be a tuple of length 3!' 'Dimension has to be a tuple of length 3!'
if mask is not None: if mask is not None:
assert mask.shape == dim, 'Mask dimensions must match!' assert mask.shape == dim, 'Mask dimensions must match!'
self.m = 3 * np.sum(self.mask) self.n = 3 * np.sum(self.mask)
else: else:
self.m = 3 * np.prod(dim) self.n = 3 * np.prod(dim)
self.a = a self.a = a
self.dim = dim self.dim = dim
self.b_0 = b_0 self.b_0 = b_0
......
...@@ -30,9 +30,9 @@ class ForwardModel(object): ...@@ -30,9 +30,9 @@ class ForwardModel(object):
The grid spacing in nm. The grid spacing in nm.
dim : tuple (N=3) dim : tuple (N=3)
Dimensions of the 3D magnetic distribution. Dimensions of the 3D magnetic distribution.
n: int
Size of the image space. Number of pixels of the 2-dimensional projected grid.
m: int m: int
Size of the image space. Number of pixels of the 2-dimensional projected grid.
n: int
Size of the input space. Number of voxels of the 3-dimensional grid. Size of the input space. Number of voxels of the 3-dimensional grid.
''' '''
...@@ -61,7 +61,7 @@ class ForwardModel(object): ...@@ -61,7 +61,7 @@ class ForwardModel(object):
self.LOG.debug('Calling __call__') self.LOG.debug('Calling __call__')
self.mag_data.set_vector(x, self.data_set.mask) self.mag_data.set_vector(x, self.data_set.mask)
# TODO: Multiprocessing # TODO: Multiprocessing
result = np.zeros(self.n) result = np.zeros(self.m)
hp = self.hook_points hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors): for i, projector in enumerate(self.data_set.projectors):
phase_map = self.phase_mappers[projector.dim_uv](projector(self.mag_data)) phase_map = self.phase_mappers[projector.dim_uv](projector(self.mag_data))
...@@ -90,7 +90,7 @@ class ForwardModel(object): ...@@ -90,7 +90,7 @@ class ForwardModel(object):
''' '''
self.LOG.debug('Calling jac_dot') self.LOG.debug('Calling jac_dot')
self.mag_data.set_vector(vector, self.data_set.mask) self.mag_data.set_vector(vector, self.data_set.mask)
result = np.zeros(self.n) result = np.zeros(self.m)
hp = self.hook_points hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors): for i, projector in enumerate(self.data_set.projectors):
mag_vec = self.mag_data.mag_vec mag_vec = self.mag_data.mag_vec
...@@ -119,7 +119,7 @@ class ForwardModel(object): ...@@ -119,7 +119,7 @@ class ForwardModel(object):
''' '''
self.LOG.debug('Calling jac_T_dot') self.LOG.debug('Calling jac_T_dot')
result = np.zeros(self.m) result = np.zeros(self.n)
hp = self.hook_points hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors): for i, projector in enumerate(self.data_set.projectors):
vec = vector[hp[i]:hp[i+1]] vec = vector[hp[i]:hp[i+1]]
......
...@@ -281,9 +281,9 @@ class MagData(object): ...@@ -281,9 +281,9 @@ class MagData(object):
Order is: first all `x`-, then all `y`-, then all `z`-components. Order is: first all `x`-, then all `y`-, then all `z`-components.
''' '''
return np.reshape([self.magnitude[2][mask], return np.reshape([self.magnitude[0][mask],
self.magnitude[1][mask], self.magnitude[1][mask],
self.magnitude[0][mask]], -1) self.magnitude[2][mask]], -1)
def set_vector(self, vector, mask=None): def set_vector(self, vector, mask=None):
'''Set the magnetic components of the masked pixels to the values specified by `vector`. '''Set the magnetic components of the masked pixels to the values specified by `vector`.
...@@ -304,9 +304,9 @@ class MagData(object): ...@@ -304,9 +304,9 @@ class MagData(object):
assert np.size(vector) % 3 == 0, 'Vector has to contain all 3 components for every pixel!' assert np.size(vector) % 3 == 0, 'Vector has to contain all 3 components for every pixel!'
count = np.size(vector)/3 count = np.size(vector)/3
if mask is not None: if mask is not None:
self.magnitude[2][mask] = vector[:count] # x-component self.magnitude[0][mask] = vector[:count] # x-component
self.magnitude[1][mask] = vector[count:2*count] # y-component self.magnitude[1][mask] = vector[count:2*count] # y-component
self.magnitude[0][mask] = vector[2*count:] # z-component self.magnitude[2][mask] = vector[2*count:] # z-component
else: else:
self.mag_vec = vector self.mag_vec = vector
......
...@@ -75,9 +75,9 @@ class PhaseMapperRDFC(PhaseMapper): ...@@ -75,9 +75,9 @@ class PhaseMapperRDFC(PhaseMapper):
geometry : {'disc', 'slab'}, optional geometry : {'disc', 'slab'}, optional
Elementary geometry which is used for the phase contribution of one pixel. Elementary geometry which is used for the phase contribution of one pixel.
Default is 'disc'. Default is 'disc'.
n: int
Size of the image space.
m: int m: int
Size of the image space.
n: int
Size of the input space. Size of the input space.
''' '''
...@@ -87,8 +87,8 @@ class PhaseMapperRDFC(PhaseMapper): ...@@ -87,8 +87,8 @@ class PhaseMapperRDFC(PhaseMapper):
def __init__(self, kernel): def __init__(self, kernel):
self.LOG.debug('Calling __init__') self.LOG.debug('Calling __init__')
self.kernel = kernel self.kernel = kernel
self.n = np.prod(kernel.dim_uv) self.m = np.prod(kernel.dim_uv)
self.m = 2 * self.n self.n = 2 * self.m
self.LOG.debug('Created '+str(self)) self.LOG.debug('Created '+str(self))
def __repr__(self): def __repr__(self):
...@@ -135,8 +135,8 @@ class PhaseMapperRDFC(PhaseMapper): ...@@ -135,8 +135,8 @@ class PhaseMapperRDFC(PhaseMapper):
''' '''
self.LOG.debug('Calling jac_dot') self.LOG.debug('Calling jac_dot')
assert len(vector) == self.m, \ assert len(vector) == self.n, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m) 'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
u_mag, v_mag = np.reshape(vector, (2,)+self.kernel.dim_uv) u_mag, v_mag = np.reshape(vector, (2,)+self.kernel.dim_uv)
return self._convolve(u_mag, v_mag) return self._convolve(u_mag, v_mag)
...@@ -157,10 +157,10 @@ class PhaseMapperRDFC(PhaseMapper): ...@@ -157,10 +157,10 @@ class PhaseMapperRDFC(PhaseMapper):
''' '''
self.LOG.debug('Calling jac_T_dot') self.LOG.debug('Calling jac_T_dot')
assert len(vector) == self.n, \ assert len(vector) == self.m, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n) 'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
# TODO: directly? ask Jörn again! # TODO: directly? ask Jörn again!
result = np.zeros(self.m) result = np.zeros(self.n)
nc.jac_T_dot_real_convolve(self.kernel.dim_uv[0], self.kernel.dim_uv[1], nc.jac_T_dot_real_convolve(self.kernel.dim_uv[0], self.kernel.dim_uv[1],
self.kernel.u, self.kernel.v, vector, result) self.kernel.u, self.kernel.v, vector, result)
return result return result
...@@ -204,8 +204,8 @@ class PhaseMapperRDRC(PhaseMapper): ...@@ -204,8 +204,8 @@ class PhaseMapperRDRC(PhaseMapper):
self.kernel = kernel self.kernel = kernel
self.threshold = threshold self.threshold = threshold
self.numcore = numcore self.numcore = numcore
self.n = np.prod(kernel.dim_uv) self.m = np.prod(kernel.dim_uv)
self.m = 2 * self.n self.n = 2 * self.m
self.LOG.debug('Created '+str(self)) self.LOG.debug('Created '+str(self))
def __repr__(self): def __repr__(self):
...@@ -266,15 +266,15 @@ class PhaseMapperRDRC(PhaseMapper): ...@@ -266,15 +266,15 @@ class PhaseMapperRDRC(PhaseMapper):
''' '''
self.LOG.debug('Calling jac_dot') self.LOG.debug('Calling jac_dot')
assert len(vector) == self.m, \ assert len(vector) == self.n, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m) 'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
v_dim, u_dim = self.kernel.dim_uv v_dim, u_dim = self.kernel.dim_uv
result = np.zeros(self.n) result = np.zeros(self.m)
if self.numcore: if self.numcore:
nc.jac_dot_real_convolve(v_dim, u_dim, self.kernel.u, self.kernel.v, vector, result) nc.jac_dot_real_convolve(v_dim, u_dim, self.kernel.u, self.kernel.v, vector, result)
else: else:
# Iterate over all contributing pixels (numbered consecutively) # Iterate over all contributing pixels (numbered consecutively)
for s in range(self.n): # column-wise (two columns at a time, u- and v-component) for s in range(self.m): # column-wise (two columns at a time, u- and v-component)
i = s % u_dim # u-coordinate of current contributing pixel i = s % u_dim # u-coordinate of current contributing pixel
j = int(s/u_dim) # v-coordinate of current ccontributing pixel j = int(s/u_dim) # v-coordinate of current ccontributing pixel
u_min = (u_dim-1) - i # u_dim-1: center of the kernel u_min = (u_dim-1) - i # u_dim-1: center of the kernel
...@@ -282,7 +282,7 @@ class PhaseMapperRDRC(PhaseMapper): ...@@ -282,7 +282,7 @@ class PhaseMapperRDRC(PhaseMapper):
v_min = (v_dim-1) - j # v_dim-1: center of the kernel v_min = (v_dim-1) - j # v_dim-1: center of the kernel
v_max = (2*v_dim-1) - j # = v_min + v_dim v_max = (2*v_dim-1) - j # = v_min + v_dim
result += vector[s] * self.kernel.u[v_min:v_max, u_min:u_max].reshape(-1) result += vector[s] * self.kernel.u[v_min:v_max, u_min:u_max].reshape(-1)
result -= vector[s+self.n] * self.kernel.v[v_min:v_max, u_min:u_max].reshape(-1) result -= vector[s+self.m] * self.kernel.v[v_min:v_max, u_min:u_max].reshape(-1)
return result return result
def jac_T_dot(self, vector): def jac_T_dot(self, vector):
...@@ -302,15 +302,15 @@ class PhaseMapperRDRC(PhaseMapper): ...@@ -302,15 +302,15 @@ class PhaseMapperRDRC(PhaseMapper):
''' '''
self.LOG.debug('Calling jac_T_dot') self.LOG.debug('Calling jac_T_dot')
assert len(vector) == self.n, \ assert len(vector) == self.m, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n) 'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
v_dim, u_dim = self.dim_uv v_dim, u_dim = self.dim_uv
result = np.zeros(self.m) result = np.zeros(self.n)
if self.numcore: if self.numcore:
nc.jac_T_dot_real_convolve(v_dim, u_dim, self.kernel.u, self.kernel.v, vector, result) nc.jac_T_dot_real_convolve(v_dim, u_dim, self.kernel.u, self.kernel.v, vector, result)
else: else:
# Iterate over all contributing pixels (numbered consecutively): # Iterate over all contributing pixels (numbered consecutively):
for s in range(self.n): # row-wise (two rows at a time, u- and v-component) for s in range(self.m): # row-wise (two rows at a time, u- and v-component)
i = s % u_dim # u-coordinate of current contributing pixel i = s % u_dim # u-coordinate of current contributing pixel
j = int(s/u_dim) # v-coordinate of current contributing pixel j = int(s/u_dim) # v-coordinate of current contributing pixel
u_min = (u_dim-1) - i # u_dim-1: center of the kernel u_min = (u_dim-1) - i # u_dim-1: center of the kernel
...@@ -318,7 +318,7 @@ class PhaseMapperRDRC(PhaseMapper): ...@@ -318,7 +318,7 @@ class PhaseMapperRDRC(PhaseMapper):
v_min = (v_dim-1) - j # v_dim-1: center of the kernel v_min = (v_dim-1) - j # v_dim-1: center of the kernel
v_max = (2*v_dim-1) - j # = v_min + v_dim v_max = (2*v_dim-1) - j # = v_min + v_dim
result[s] = np.sum(vector*self.u[v_min:v_max, u_min:u_max].reshape(-1)) result[s] = np.sum(vector*self.u[v_min:v_max, u_min:u_max].reshape(-1))
result[s+self.n] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1)) result[s+self.m] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1))
return result return result
...@@ -356,8 +356,8 @@ class PhaseMapperFDFC(PhaseMapper): ...@@ -356,8 +356,8 @@ class PhaseMapperFDFC(PhaseMapper):
self.dim_uv = dim_uv self.dim_uv = dim_uv
self.b_0 = b_0 self.b_0 = b_0
self.padding = padding self.padding = padding
self.n = np.prod(dim_uv) self.m = np.prod(dim_uv)
self.m = 2 * self.n self.n = 2 * self.m
self.LOG.debug('Created '+str(self)) self.LOG.debug('Created '+str(self))
def __repr__(self): def __repr__(self):
...@@ -416,8 +416,8 @@ class PhaseMapperFDFC(PhaseMapper): ...@@ -416,8 +416,8 @@ class PhaseMapperFDFC(PhaseMapper):
''' '''
self.LOG.debug('Calling jac_dot') self.LOG.debug('Calling jac_dot')
assert len(vector) == self.m, \ assert len(vector) == self.n, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m) 'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
mag_proj = MagData(self.a, np.zeros((3, 1)+self.dim_uv)) mag_proj = MagData(self.a, np.zeros((3, 1)+self.dim_uv))
magnitude_proj = self.jac_dot(vector).reshape((2, )+self.dim_uv) magnitude_proj = self.jac_dot(vector).reshape((2, )+self.dim_uv)
mag_proj.magnitude[1:3, 0, ...] = magnitude_proj mag_proj.magnitude[1:3, 0, ...] = magnitude_proj
......
...@@ -45,9 +45,9 @@ class Projector(object): ...@@ -45,9 +45,9 @@ class Projector(object):
coeff : list (N=2) coeff : list (N=2)
List containing the six weighting coefficients describing the influence of the 3 components List containing the six weighting coefficients describing the influence of the 3 components
of a 3-dimensional vector field on the 2 projected components. of a 3-dimensional vector field on the 2 projected components.
n: int
Size of the image space.
m: int m: int
Size of the image space.
n: int
Size of the input space. Size of the input space.
''' '''
...@@ -63,8 +63,8 @@ class Projector(object): ...@@ -63,8 +63,8 @@ class Projector(object):
self.weight = weight self.weight = weight
self.coeff = coeff self.coeff = coeff
self.size_2d, self.size_3d = weight.shape self.size_2d, self.size_3d = weight.shape
self.m = 3 * np.prod(dim) self.n = 3 * np.prod(dim)
self.n = 2 * np.prod(dim_uv) self.m = 2 * np.prod(dim_uv)
self.LOG.debug('Created '+str(self)) self.LOG.debug('Created '+str(self))
def __repr__(self): def __repr__(self):
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment