From ecfe52f674c98004bc5c752355bddbe320a3d6dc Mon Sep 17 00:00:00 2001
From: "Joern Ungermann (IEK-7 FZ-Juelich)" <j.ungermann@fz-juelich.de>
Date: Thu, 30 Oct 2014 00:37:09 +0100
Subject: [PATCH] fixed some bugs, got reconstruction to work.
---
pyramid/costfunction.py | 7 +++++--
pyramid/dataset.py | 2 +-
pyramid/forwardmodel.py | 4 +++-
pyramid/kernel.py | 35 +++++++++++++++++++++++++++++++++++
pyramid/magdata.py | 15 ++++++++++-----
pyramid/phasemap.py | 1 +
pyramid/phasemapper.py | 12 +++++++-----
pyramid/reconstruction.py | 33 ++++++++++++++++++++++++++++-----
pyramid/regularisator.py | 36 +++++++++++-------------------------
9 files changed, 101 insertions(+), 44 deletions(-)
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index ad7c226..50df34d 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -53,8 +53,9 @@ class Costfunction(object):
# Extract important information:
self.y = data_set.phase_vec
self.Se_inv = data_set.Se_inv
- self.n = data_set.n * 3
+ self.n = data_set.m
self.m = len(self.y)
+ self.chisq, self.chisq_a, self.chisq_m = None, None, None
self.LOG.debug('Created '+str(self))
def __repr__(self):
@@ -73,7 +74,9 @@ class Costfunction(object):
def __call__(self, x):
self.LOG.debug('Calling __call__')
delta_y = self.fwd_model(x) - self.y
- self.chisq = delta_y.dot(self.Se_inv.dot(delta_y)) + self.regularisator(x)
+ self.chisq_m = delta_y.dot(self.Se_inv.dot(delta_y))
+ self.chisq_a = self.regularisator(x)
+ self.chisq = self.chisq_m + self.chisq_a
return self.chisq
def jac(self, x):
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index fd5fb49..914d306 100644
--- a/pyramid/dataset.py
+++ b/pyramid/dataset.py
@@ -88,7 +88,7 @@ class DataSet(object):
'Dimension has to be a tuple of length 3!'
if mask is not None:
assert mask.shape == dim, 'Mask dimensions must match!'
- self.m = 3 * np.sum(self.mask)
+ self.m = 3 * np.sum(mask)
else:
self.m = 3 * np.prod(dim)
self.a = a
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index c48a7f7..758a692 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -59,6 +59,7 @@ class ForwardModel(object):
def __call__(self, x):
self.LOG.debug('Calling __call__')
+ self.mag_data.magnitude[:] = 0
self.mag_data.set_vector(x, self.data_set.mask)
# TODO: Multiprocessing
result = np.zeros(self.n)
@@ -89,6 +90,7 @@ class ForwardModel(object):
'''
self.LOG.debug('Calling jac_dot')
+ self.mag_data.magnitude[:] = 0
self.mag_data.set_vector(vector, self.data_set.mask)
result = np.zeros(self.n)
hp = self.hook_points
@@ -119,7 +121,7 @@ class ForwardModel(object):
'''
self.LOG.debug('Calling jac_T_dot')
- result = np.zeros(self.m)
+ result = np.zeros(3*np.prod(self.data_set.dim))
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
vec = vector[hp[i]:hp[i+1]]
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index be8502d..0c48d21 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -105,3 +105,38 @@ class Kernel(object):
self.LOG.debug('Calling __str__')
return 'Kernel(a=%s, dim_uv=%s, numcore=%s, geometry=%s)' % \
(self.a, self.dim_uv, self.numcore, self.geometry)
+
+ def _multiply_jacobi(self, vector):
+ self.LOG.debug('Calling _multiply_jacobi')
+ v_dim, u_dim = self.dim_uv
+ assert len(vector) == 2 * self.size, \
+ 'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.size)
+ result = np.zeros(self.size)
+ # Iterate over all contributing pixels (numbered consecutively)
+ for s in range(self.size): # column-wise (two columns at a time, u- and v-component)
+ i = s % u_dim # u-coordinate of current contributing 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_max = (2*u_dim-1) - i # = u_min + u_dim
+ v_min = (v_dim-1) - j # v_dim-1: center of the kernel
+ v_max = (2*v_dim-1) - j # = v_min + v_dim
+ result += vector[s] * self.u[v_min:v_max, u_min:u_max].reshape(-1) # u
+ result -= vector[s+self.size] * self.v[v_min:v_max, u_min:u_max].reshape(-1) # v
+ return result
+
+ def _multiply_jacobi_T(self, vector):
+ self.LOG.debug('Calling _multiply_jacobi_T')
+ v_dim, u_dim = self.dim_uv
+ result = np.zeros(2*self.size)
+ # Iterate over all contributing pixels (numbered consecutively):
+ for s in range(self.size): # row-wise (two rows at a time, u- and v-component)
+ i = s % u_dim # u-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_max = (2*u_dim-1) - i # = u_min + u_dim
+ v_min = (v_dim-1) - j # v_dim-1: center of the kernel
+ 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)) # u
+ result[s+self.size] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1)) # v
+ return result
+#
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index 5989153..fc6f4d4 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -83,7 +83,8 @@ class MagData(object):
@mag_vec.setter
def mag_vec(self, mag_vec):
assert isinstance(mag_vec, np.ndarray), 'Vector has to be a numpy array!'
- assert np.size(mag_vec) == 3*np.prod(self.dim), 'Vector has to match magnitude dimensions!'
+ assert np.size(mag_vec) == 3*np.prod(self.dim), \
+ 'Vector has to match magnitude dimensions! {} {}'.format(mag_vec.shape, 3*np.prod(self.dim))
self.magnitude = mag_vec.reshape((3,)+self.dim)
def __init__(self, a, magnitude):
@@ -281,9 +282,13 @@ class MagData(object):
Order is: first all `x`-, then all `y`-, then all `z`-components.
'''
- return np.reshape([self.magnitude[2][mask],
+ if mask is not None:
+ return np.reshape([self.magnitude[0][mask],
self.magnitude[1][mask],
- self.magnitude[0][mask]], -1)
+ self.magnitude[2][mask]], -1)
+ else:
+ return self.mag_vec
+
def set_vector(self, vector, mask=None):
'''Set the magnetic components of the masked pixels to the values specified by `vector`.
@@ -304,9 +309,9 @@ class MagData(object):
assert np.size(vector) % 3 == 0, 'Vector has to contain all 3 components for every pixel!'
count = np.size(vector)/3
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[0][mask] = vector[2*count:] # z-component
+ self.magnitude[2][mask] = vector[2*count:] # z-component
else:
self.mag_vec = vector
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index af1e708..f5a3a4a 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -293,6 +293,7 @@ class PhaseMap(object):
phase_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
a = phase_file.a
phase = phase_file.variables['phase'][:]
+ #phase = phase[28:36, 28:36]
phase_file.close()
return PhaseMap(a, phase)
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 1d4205a..343a876 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -138,7 +138,9 @@ class PhaseMapperRDFC(PhaseMapper):
assert len(vector) == self.m, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
u_mag, v_mag = np.reshape(vector, (2,)+self.kernel.dim_uv)
- return self._convolve(u_mag, v_mag)
+ result = self._convolve(u_mag, v_mag).reshape(-1)
+ return result
+ #return self.kernel._multiply_jacobi(vector)
def jac_T_dot(self, vector):
'''Calculate the product of the transposed Jacobi matrix with a given `vector`.
@@ -161,10 +163,10 @@ class PhaseMapperRDFC(PhaseMapper):
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
# TODO: directly? ask Jörn again!
result = np.zeros(self.m)
+ v_dim, u_dim = self.kernel.dim_uv
nc.jac_T_dot_real_convolve(self.kernel.dim_uv[0], self.kernel.dim_uv[1],
self.kernel.u, self.kernel.v, vector, result)
- return result
-
+ return self.kernel._multiply_jacobi_T(vector)
class PhaseMapperRDRC(PhaseMapper):
@@ -317,8 +319,8 @@ class PhaseMapperRDRC(PhaseMapper):
u_max = (2*u_dim-1) - i # = u_min + u_dim
v_min = (v_dim-1) - j # v_dim-1: center of the kernel
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+self.n] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1))
+ result[s] = np.sum(vector*self.kernel.u[v_min:v_max, u_min:u_max].reshape(-1))
+ result[s+self.n] = np.sum(vector*-self.kernel.v[v_min:v_max, u_min:u_max].reshape(-1))
return result
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index 1eff2ea..d8d6640 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -119,11 +119,11 @@ def optimize_linear(data, regularisator=None, maxiter=1000, verbosity=0):
fwd_model = ForwardModel(data)
cost = Costfunction(data, regularisator)
print cost(np.zeros(cost.n))
- x_opt = jutil.cg.conj_grad_minimize(cost)
+ x_opt = jutil.cg.conj_grad_minimize(cost, max_iter=20)
print cost(x_opt)
# Create and return fitting MagData object:
mag_opt = MagData(data.a, np.zeros((3,)+data.dim))
- mag_opt.mag_vec = x_opt
+ mag_opt.set_vector(x_opt, data.mask)
return mag_opt
@@ -152,15 +152,38 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
LOG.debug('Calling optimize_cg')
if first_guess is None:
first_guess = MagData(data.a, np.zeros((3,)+data.dim))
- x_0 = first_guess.mag_vec
+ x_0 = first_guess.get_vector(data.mask)
fwd_model = ForwardModel(data)
cost = Costfunction(data, regularisator)
+
+# proj = fwd_model.data_set.projectors[0]
+# proj_jac1 = np.array([proj.jac_dot(np.eye(proj.m, 1, -i).squeeze()) for i in range(proj.m)])
+# proj_jac2 = np.array([proj.jac_T_dot(np.eye(proj.n, 1, -i).squeeze()) for i in range(proj.n)])
+# print jac1, jac2.T, abs(jac1-jac2.T).sum()
+# print jac1.shape, jac2.shape
+
+# pm = fwd_model.phase_mappers[proj.dim_uv]
+# pm_jac1 = np.array([pm.jac_dot(np.eye(pm.m)[:, i]) for i in range(pm.m)])
+# pm_jac2 = np.array([pm.jac_T_dot(np.eye(pm.n)[:, i]) for i in range(pm.n)])
+# print jac1, jac2.T, abs(jac1-jac2.T).sum()
+# print jac1.shape, jac2.shape
+
+
+ # jac1 = np.array([fwd_model.jac_dot(x_0, np.eye(fwd_model.m)[:, i]) for i in range(fwd_model.m)])
+ # jac2 = np.array([fwd_model.jac_T_dot(x_0, np.eye(fwd_model.n)[:, i]) for i in range(fwd_model.n)])
+ # print proj_jac1.dot(pm_jac1)
+ # print (pm_jac2.dot(proj_jac2)).T
+ # print jac1
+# print jac2.T
+# print abs(jac1-jac2.T).sum()
+# print jac1.shape, jac2.shape
+
assert len(x_0) == cost.n, (len(x_0), cost.m, cost.n)
- result = jutil.minimizer.minimize(cost, x_0, options={"conv_rel":1e-20}, tol={"max_iteration":1})
+ result = jutil.minimizer.minimize(cost, x_0, options={"conv_rel":1e-2}, tol={"max_iteration":4})
x_opt = result.x
print cost(x_opt)
mag_opt = MagData(data.a, np.zeros((3,)+data.dim))
- mag_opt.mag_vec = x_opt
+ mag_opt.set_vector(x_opt, data.mask)
return mag_opt
diff --git a/pyramid/regularisator.py b/pyramid/regularisator.py
index bd06648..5f4c2cf 100644
--- a/pyramid/regularisator.py
+++ b/pyramid/regularisator.py
@@ -110,28 +110,14 @@ class ZeroOrderRegularisator(Regularisator):
self.LOG.debug('Created '+str(self))
-#class FirstOrderRegularisator(Regularisator):
-# # TODO: Docstring!
-#
-# def __init__(self, fwd_model, lam, x_a=None):
-# size_3d = fwd_model.size_3d
-# dim = fwd_model.dim
-# converter = IndexConverter(dim)
-# row = []
-# col = []
-# data = []
-#
-# for i in range(size_3d):
-# neighbours = converter.find_neighbour_ind(i)
-#
-# Sa_inv = csr_matrix(coo_matrix(data, (rows, columns)), shape=(3*size_3d, 3*size_3d))
-#
-# term2 = []
-# for i in range(3):
-# component = mag_data[i, ...]
-# for j in range(3):
-# if component.shape[j] > 1:
-# term2.append(np.diff(component, axis=j).reshape(-1))
-#
-# super(FirstOrderRegularisator, self).__init__(Sa_inv, x_a)
-# self.LOG.debug('Created '+str(self))
+class FirstOrderRegularisator(Regularisator):
+ # TODO: Docstring!
+
+ def __init__(self, mask, lam, x_a=None):
+ import jutil
+ D0 = jutil.diff.get_diff_operator(mask, 0, 3)
+ D1 = jutil.diff.get_diff_operator(mask, 1, 3)
+ D = jutil.operator.VStack([D0, D1])
+ norm = jutil.norms.WeightedL2Square(D)
+ super(FirstOrderRegularisator, self).__init__(norm, lam)
+ self.LOG.debug('Created '+str(self))
--
GitLab