From 43dc54e7122c0b8641171a29df7947d054c1d2a1 Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Tue, 4 Nov 2014 14:00:34 +0100
Subject: [PATCH] Implementation of jac_T_dot with jutil.fft_adj!
---
pyramid/kernel.py | 38 +------
pyramid/logging.ini | 2 +-
pyramid/phasemapper.py | 61 ++++++++++--
pyramid/reconstruction.py | 2 +-
scripts/collaborations/example_joern.py | 22 ++--
scripts/collaborations/zi_an.py | 127 +++++++++++++-----------
6 files changed, 137 insertions(+), 115 deletions(-)
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index 0c48d21..63f1127 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -92,8 +92,11 @@ class Kernel(object):
dim_combined = 3*np.array(dim_uv) - 2 # (dim_uv-1) + dim_uv + (dim_uv-1) mag + kernel
self.dim_fft = 2 ** np.ceil(np.log2(dim_combined)).astype(int) # next multiple of 2
self.slice_fft = (slice(dim_uv[0]-1, 2*dim_uv[0]-1), slice(dim_uv[1]-1, 2*dim_uv[1]-1))
+ self.slice_fft_compl = (slice(2*dim_uv[0]-1, 2*dim_uv[0]-1), slice(2*dim_uv[1]-1, 2*dim_uv[1]-1))
self.u_fft = np.fft.rfftn(self.u, self.dim_fft)
self.v_fft = np.fft.rfftn(self.v, self.dim_fft)
+ self.u_fft_compl = np.fft.fftn(self.u, self.dim_fft)
+ self.v_fft_compl = np.fft.fftn(self.v, self.dim_fft)
self.LOG.debug('Created '+str(self))
def __repr__(self):
@@ -105,38 +108,3 @@ 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/logging.ini b/pyramid/logging.ini
index d206f57..3a6529f 100644
--- a/pyramid/logging.ini
+++ b/pyramid/logging.ini
@@ -17,7 +17,7 @@ keys=console,file
[handler_console]
class=logging.StreamHandler
-level=WARNING
+level=INFO
formatter=console
args=tuple()
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 991891e..34d5b28 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -82,7 +82,7 @@ class PhaseMapperRDFC(PhaseMapper):
'''
- LOG = logging.getLogger(__name__+'.PMConvolve')
+ LOG = logging.getLogger(__name__+'.PhaseMapperRDFC')
def __init__(self, kernel):
self.LOG.debug('Calling __init__')
@@ -161,10 +161,57 @@ class PhaseMapperRDFC(PhaseMapper):
self.LOG.debug('Calling jac_T_dot')
assert len(vector) == self.m, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
- # TODO: directly? ask Jörn again!
result = np.zeros(self.n)
- nc.jac_T_dot_real_convolve(self.kernel.dim_uv[0], self.kernel.dim_uv[1],
- self.kernel.u, self.kernel.v, vector, result)
+
+ import jutil.fft as jfft
+ from jutil.taketime import TakeTime
+
+ with TakeTime('fft_adjoint'):
+
+ phase = vector.reshape(self.kernel.dim_uv)
+ p0 = self.kernel.dim_fft[0]-self.kernel.dim_uv[0]
+ p1 = self.kernel.dim_fft[1]-self.kernel.dim_uv[1]
+ phase = np.pad(phase, ((0, p0), (0, p1)), 'constant')
+
+ phase_fft = jfft.irfft2_adj(phase)
+ u_mag_fft = phase_fft * self.kernel.u_fft
+ v_mag_fft = phase_fft * -self.kernel.v_fft
+ u_mag = jfft.rfft2_adj(u_mag_fft, self.kernel.dim_fft[1])[self.kernel.slice_fft]
+ v_mag = jfft.rfft2_adj(v_mag_fft, self.kernel.dim_fft[1])[self.kernel.slice_fft]
+ result = -np.concatenate((u_mag.flatten(), v_mag.flatten()))
+
+# # TODO: directly? ask Jörn again!
+# phase = vector.reshape(self.kernel.dim_uv)
+# phase_fft = np.fft.fft(phase, self.kernel.dim_fft)
+# u_mag_fft = phase_fft * self.kernel.u_fft_compl
+# v_mag_fft = phase_fft * -self.kernel.v_fft_compl
+# u_mag = np.fft.ifft(u_mag_fft, self.kernel.dim_fft)[self.kernel.slice_fft]
+# v_mag = np.fft.ifft(u_mag_fft, self.kernel.dim_fft)[self.kernel.slice_fft]
+# result = np.concatenate((u_mag.flatten(), v_mag.flatten()))
+# return result
+
+ #TODO: np.split()
+ #TODO: TakeTime()
+ #TODO: 'constant' in np.pad()
+
+# result[s] = np.sum(vector*self.u[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))
+#
+# # Fourier transform the projected magnetisation:
+# u_mag_fft = np.fft.rfftn(u_mag, self.kernel.dim_fft)
+# v_mag_fft = np.fft.rfftn(v_mag, self.kernel.dim_fft)
+# # Convolve the magnetization with the kernel in Fourier space:
+# phase_fft = u_mag_fft*self.kernel.u_fft - v_mag_fft*self.kernel.v_fft
+# # Return the result:
+# return np.fft.irfftn(phase_fft, self.kernel.dim_fft)[self.kernel.slice_fft]
+
+# with TakeTime('oldschool'):
+# compare = np.zeros(self.n)
+# nc.jac_T_dot_real_convolve(self.kernel.dim_uv[0], self.kernel.dim_uv[1],
+# self.kernel.u, self.kernel.v, vector, compare)
+
+# import pdb; pdb.set_trace()
+# assert np.testing.assert_almost_equal(result, compare)
return result
@@ -199,7 +246,7 @@ class PhaseMapperRDRC(PhaseMapper):
'''
- LOG = logging.getLogger(__name__+'.PMReal')
+ LOG = logging.getLogger(__name__+'.PhaseMapperRDRC')
def __init__(self, kernel, threshold=0, numcore=True):
self.LOG.debug('Calling __init__')
@@ -349,7 +396,7 @@ class PhaseMapperFDFC(PhaseMapper):
'''
- LOG = logging.getLogger(__name__+'.PMFourier')
+ LOG = logging.getLogger(__name__+'.PhaseMapperFDFC')
PHI_0 = -2067.83 # magnetic flux in T*nm²
def __init__(self, a, dim_uv, b_0=1, padding=0):
@@ -458,7 +505,7 @@ class PhaseMapperElectric(PhaseMapper):
'''
- LOG = logging.getLogger(__name__+'.PMElectric')
+ LOG = logging.getLogger(__name__+'.PhaseMapperElectric')
H_BAR = 6.626E-34 # Planck constant in J*s
M_E = 9.109E-31 # electron mass in kg
Q_E = 1.602E-19 # electron charge in C
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index 39948e2..bdbb818 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -109,7 +109,7 @@ def optimize_linear(data, regularisator=None, max_iter=None):
# Set up necessary objects:
cost = Costfunction(data, regularisator)
LOG.info('Cost before optimization: {}'.format(cost(np.zeros(cost.n))))
- x_opt = jcg.conj_grad_minimize(cost, max_iter=max_iter)
+ x_opt = jcg.conj_grad_minimize(cost, max_iter=max_iter).x
LOG.info('Cost after optimization: {}'.format(cost(x_opt)))
# Create and return fitting MagData object:
mag_opt = MagData(data.a, np.zeros((3,) + data.dim))
diff --git a/scripts/collaborations/example_joern.py b/scripts/collaborations/example_joern.py
index 4306d10..a80e2a7 100644
--- a/scripts/collaborations/example_joern.py
+++ b/scripts/collaborations/example_joern.py
@@ -54,7 +54,7 @@ else:
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
-logging.basicConfig(level=logging.INFO)
+#logging.basicConfig(level=logging.INFO)
###################################################################################################
threshold = 1
@@ -62,18 +62,16 @@ a = 1.0 # in nm
gain = 5
b_0 = 1
inter = 'none'
-dim = (1,) + (64, 64)
-dim_small = (64, 64)
smoothed_pictures = True
lam = 1E-6
log = True
PATH = '../../output/joern/'
###################################################################################################
# Read in files:
-phase_map = PhaseMap.load_from_netcdf4(PATH+'phase_map.nc')
-#mask = np.genfromtxt('mask.txt', dtype=bool)
-with open(PATH + 'mask.pickle') as pf:
+phase_map = PhaseMap.load_from_netcdf4(PATH+'phase_map_2.nc')
+with open(PATH + 'mask_2.pickle') as pf:
mask = pickle.load(pf)
+dim = mask.shape
# Setup:
if not use_mask:
mask = np.ones_like(mask, dtype=bool)
@@ -109,8 +107,8 @@ plt.savefig(dirname + "/reconstr.png")
# Plot the magnetization:
axis = (mag_data_rec*(1/mag_data_rec.magnitude.max())).quiver_plot(show=False)
-axis.set_xlim(20, 45)
-axis.set_ylim(20, 45)
+axis.set_xlim(int(20/64*dim[1], 45/64*dim[2])
+axis.set_ylim(int(20/64*dim[1], 45/64*dim[2])
plt.savefig(dirname + "/quiver.png")
# Display the Phase:
@@ -125,14 +123,14 @@ axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1,
interpolation=inter, show=False)
mag_data_rec.quiver_plot(axis=axis, show=False)
axis = plt.gca()
-axis.set_xlim(20, 45)
-axis.set_ylim(20, 45)
+axis.set_xlim(int(20/64*dim[1], 45/64*dim[2])
+axis.set_ylim(int(20/64*dim[1], 45/64*dim[2])
plt.savefig(dirname + "/overlay_normal.png")
axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1,
interpolation=inter, show=False)
mag_data_rec.quiver_plot(axis=axis, log=log, show=False)
axis = plt.gca()
-axis.set_xlim(20, 45)
-axis.set_ylim(20, 45)
+axis.set_xlim(int(20/64*dim[1], 45/64*dim[2])
+axis.set_ylim(int(20/64*dim[1], 45/64*dim[2])
plt.savefig(dirname + "/overlay_log.png")
diff --git a/scripts/collaborations/zi_an.py b/scripts/collaborations/zi_an.py
index 24cb95e..c27b041 100644
--- a/scripts/collaborations/zi_an.py
+++ b/scripts/collaborations/zi_an.py
@@ -9,6 +9,8 @@ import os
import numpy as np
+import pickle
+
from pyDM3reader import DM3lib as dm3
import matplotlib.pyplot as plt
@@ -34,7 +36,7 @@ a = 1.0 # in nm
gain = 5
b_0 = 1
inter = 'none'
-dim_small = (64, 64)
+dim_small = (512, 512)
smoothed_pictures = True
lam = 1E-4
order = 1
@@ -65,61 +67,68 @@ phase_map_4.display_combined(gain=gain, interpolation=inter)
plt.savefig(PATH+'%.0e/phase_map_4part.png'%lam)
mask_4 = np.expand_dims(np.where(np.array(dm3_4_ele.resize(dim_small)) >= threshold,
True, False), axis=0)
-# Reconstruct the magnetic distribution:
-tic = clock()
-mag_data_rec_2 = rc.optimize_simple_leastsq(phase_map_2, mask_2, b_0, lam=lam, order=order)
-print '2 particle reconstruction time:', clock() - tic
-tic = clock()
-mag_data_rec_4 = rc.optimize_simple_leastsq(phase_map_4, mask_4, b_0, lam=lam, order=order)
-print '4 particle reconstruction time:', clock() - tic
-# Display the reconstructed phase map and holography image:
-phase_map_rec_2 = pm(mag_data_rec_2)
-phase_map_rec_2.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter)
-plt.savefig(PATH+'%.0e/phase_map_2part_rec.png'%lam)
-phase_map_rec_4 = pm(mag_data_rec_4)
-phase_map_rec_4.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter)
-plt.savefig(PATH+'%.0e/phase_map_4part_rec.png'%lam)
-# Plot the magnetization:
-axis = (mag_data_rec_2*(1/mag_data_rec_2.magnitude.max())).quiver_plot()
-axis.set_xlim(20, 45)
-axis.set_ylim(20, 45)
-plt.savefig(PATH+'%.0e/mag_data_2part.png'%lam)
-axis = (mag_data_rec_4*(1/mag_data_rec_4.magnitude.max())).quiver_plot()
-axis.set_xlim(20, 45)
-axis.set_ylim(20, 45)
-plt.savefig(PATH+'%.0e/mag_data_4part.png'%lam)
-# Display the Phase:
-phase_diff_2 = phase_map_rec_2-phase_map_2
-phase_diff_2.display_phase('Difference')
-plt.savefig(PATH+'%.0e/phase_map_2part_diff.png'%lam)
-phase_diff_4 = phase_map_rec_4-phase_map_4
-phase_diff_4.display_phase('Difference')
-plt.savefig(PATH+'%.0e/phase_map_4part_diff.png'%lam)
-# Get the average difference from the experimental results:
-print 'Average difference (2 cubes):', np.average(phase_diff_2.phase)
-print 'Average difference (4 cubes):', np.average(phase_diff_4.phase)
-# Plot holographic contour maps with overlayed magnetic distributions:
-axis = phase_map_rec_2.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
-mag_data_rec_2.quiver_plot(axis=axis)
-axis = plt.gca()
-axis.set_xlim(20, 45)
-axis.set_ylim(20, 45)
-plt.savefig(PATH+'%.0e/phase_map_2part_holo.png'%lam)
-axis = phase_map_rec_4.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
-mag_data_rec_4.quiver_plot(axis=axis)
-axis = plt.gca()
-axis.set_xlim(20, 45)
-axis.set_ylim(20, 45)
-plt.savefig(PATH+'%.0e/phase_map_4part_holo.png'%lam)
-axis = phase_map_rec_2.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
-mag_data_rec_2.quiver_plot(axis=axis, log=log)
-axis = plt.gca()
-axis.set_xlim(20, 45)
-axis.set_ylim(20, 45)
-plt.savefig(PATH+'%.0e/phase_map_2part_holo_log.png'%lam)
-axis = phase_map_rec_4.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
-mag_data_rec_4.quiver_plot(axis=axis, log=log)
-axis = plt.gca()
-axis.set_xlim(20, 45)
-axis.set_ylim(20, 45)
-plt.savefig(PATH+'%.0e/phase_map_4part_holo_log.png'%lam)
+phase_map_2.save_to_netcdf4('../../output/joern/phase_map_2.nc')
+phase_map_4.save_to_netcdf4('../../output/joern/phase_map_4.nc')
+with open('../../output/joern/mask_2.pickle', 'wb') as pf:
+ pickle.dump(mask_2, pf)
+with open('../../output/joern/mask_4.pickle', 'wb') as pf:
+ pickle.dump(mask_4, pf)
+
+## Reconstruct the magnetic distribution:
+#tic = clock()
+#mag_data_rec_2 = rc.optimize_simple_leastsq(phase_map_2, mask_2, b_0, lam=lam, order=order)
+#print '2 particle reconstruction time:', clock() - tic
+#tic = clock()
+#mag_data_rec_4 = rc.optimize_simple_leastsq(phase_map_4, mask_4, b_0, lam=lam, order=order)
+#print '4 particle reconstruction time:', clock() - tic
+## Display the reconstructed phase map and holography image:
+#phase_map_rec_2 = pm(mag_data_rec_2)
+#phase_map_rec_2.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter)
+#plt.savefig(PATH+'%.0e/phase_map_2part_rec.png'%lam)
+#phase_map_rec_4 = pm(mag_data_rec_4)
+#phase_map_rec_4.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter)
+#plt.savefig(PATH+'%.0e/phase_map_4part_rec.png'%lam)
+## Plot the magnetization:
+#axis = (mag_data_rec_2*(1/mag_data_rec_2.magnitude.max())).quiver_plot()
+#axis.set_xlim(20, 45)
+#axis.set_ylim(20, 45)
+#plt.savefig(PATH+'%.0e/mag_data_2part.png'%lam)
+#axis = (mag_data_rec_4*(1/mag_data_rec_4.magnitude.max())).quiver_plot()
+#axis.set_xlim(20, 45)
+#axis.set_ylim(20, 45)
+#plt.savefig(PATH+'%.0e/mag_data_4part.png'%lam)
+## Display the Phase:
+#phase_diff_2 = phase_map_rec_2-phase_map_2
+#phase_diff_2.display_phase('Difference')
+#plt.savefig(PATH+'%.0e/phase_map_2part_diff.png'%lam)
+#phase_diff_4 = phase_map_rec_4-phase_map_4
+#phase_diff_4.display_phase('Difference')
+#plt.savefig(PATH+'%.0e/phase_map_4part_diff.png'%lam)
+## Get the average difference from the experimental results:
+#print 'Average difference (2 cubes):', np.average(phase_diff_2.phase)
+#print 'Average difference (4 cubes):', np.average(phase_diff_4.phase)
+## Plot holographic contour maps with overlayed magnetic distributions:
+#axis = phase_map_rec_2.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
+#mag_data_rec_2.quiver_plot(axis=axis)
+#axis = plt.gca()
+#axis.set_xlim(20, 45)
+#axis.set_ylim(20, 45)
+#plt.savefig(PATH+'%.0e/phase_map_2part_holo.png'%lam)
+#axis = phase_map_rec_4.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
+#mag_data_rec_4.quiver_plot(axis=axis)
+#axis = plt.gca()
+#axis.set_xlim(20, 45)
+#axis.set_ylim(20, 45)
+#plt.savefig(PATH+'%.0e/phase_map_4part_holo.png'%lam)
+#axis = phase_map_rec_2.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
+#mag_data_rec_2.quiver_plot(axis=axis, log=log)
+#axis = plt.gca()
+#axis.set_xlim(20, 45)
+#axis.set_ylim(20, 45)
+#plt.savefig(PATH+'%.0e/phase_map_2part_holo_log.png'%lam)
+#axis = phase_map_rec_4.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
+#mag_data_rec_4.quiver_plot(axis=axis, log=log)
+#axis = plt.gca()
+#axis.set_xlim(20, 45)
+#axis.set_ylim(20, 45)
+#plt.savefig(PATH+'%.0e/phase_map_4part_holo_log.png'%lam)
--
GitLab