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

Implementation of jac_T_dot with jutil.fft_adj!

parent 17e386ed
No related branches found
No related tags found
No related merge requests found
...@@ -92,8 +92,11 @@ class Kernel(object): ...@@ -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 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.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 = (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.u_fft = np.fft.rfftn(self.u, self.dim_fft)
self.v_fft = np.fft.rfftn(self.v, 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)) self.LOG.debug('Created '+str(self))
def __repr__(self): def __repr__(self):
...@@ -105,38 +108,3 @@ class Kernel(object): ...@@ -105,38 +108,3 @@ class Kernel(object):
self.LOG.debug('Calling __str__') self.LOG.debug('Calling __str__')
return 'Kernel(a=%s, dim_uv=%s, numcore=%s, geometry=%s)' % \ return 'Kernel(a=%s, dim_uv=%s, numcore=%s, geometry=%s)' % \
(self.a, self.dim_uv, self.numcore, self.geometry) (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
#
...@@ -17,7 +17,7 @@ keys=console,file ...@@ -17,7 +17,7 @@ keys=console,file
[handler_console] [handler_console]
class=logging.StreamHandler class=logging.StreamHandler
level=WARNING level=INFO
formatter=console formatter=console
args=tuple() args=tuple()
......
...@@ -82,7 +82,7 @@ class PhaseMapperRDFC(PhaseMapper): ...@@ -82,7 +82,7 @@ class PhaseMapperRDFC(PhaseMapper):
''' '''
LOG = logging.getLogger(__name__+'.PMConvolve') LOG = logging.getLogger(__name__+'.PhaseMapperRDFC')
def __init__(self, kernel): def __init__(self, kernel):
self.LOG.debug('Calling __init__') self.LOG.debug('Calling __init__')
...@@ -161,10 +161,57 @@ class PhaseMapperRDFC(PhaseMapper): ...@@ -161,10 +161,57 @@ class PhaseMapperRDFC(PhaseMapper):
self.LOG.debug('Calling jac_T_dot') self.LOG.debug('Calling jac_T_dot')
assert len(vector) == self.m, \ assert len(vector) == self.m, \
'vector size not compatible! vector: {}, size: {}'.format(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) 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 return result
...@@ -199,7 +246,7 @@ class PhaseMapperRDRC(PhaseMapper): ...@@ -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): def __init__(self, kernel, threshold=0, numcore=True):
self.LOG.debug('Calling __init__') self.LOG.debug('Calling __init__')
...@@ -349,7 +396,7 @@ class PhaseMapperFDFC(PhaseMapper): ...@@ -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² PHI_0 = -2067.83 # magnetic flux in T*nm²
def __init__(self, a, dim_uv, b_0=1, padding=0): def __init__(self, a, dim_uv, b_0=1, padding=0):
...@@ -458,7 +505,7 @@ class PhaseMapperElectric(PhaseMapper): ...@@ -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 H_BAR = 6.626E-34 # Planck constant in J*s
M_E = 9.109E-31 # electron mass in kg M_E = 9.109E-31 # electron mass in kg
Q_E = 1.602E-19 # electron charge in C Q_E = 1.602E-19 # electron charge in C
......
...@@ -109,7 +109,7 @@ def optimize_linear(data, regularisator=None, max_iter=None): ...@@ -109,7 +109,7 @@ def optimize_linear(data, regularisator=None, max_iter=None):
# Set up necessary objects: # Set up necessary objects:
cost = Costfunction(data, regularisator) cost = Costfunction(data, regularisator)
LOG.info('Cost before optimization: {}'.format(cost(np.zeros(cost.n)))) 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))) LOG.info('Cost after optimization: {}'.format(cost(x_opt)))
# Create and return fitting MagData object: # Create and return fitting MagData object:
mag_opt = MagData(data.a, np.zeros((3,) + data.dim)) mag_opt = MagData(data.a, np.zeros((3,) + data.dim))
......
...@@ -54,7 +54,7 @@ else: ...@@ -54,7 +54,7 @@ else:
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini') 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.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
logging.basicConfig(level=logging.INFO) #logging.basicConfig(level=logging.INFO)
################################################################################################### ###################################################################################################
threshold = 1 threshold = 1
...@@ -62,18 +62,16 @@ a = 1.0 # in nm ...@@ -62,18 +62,16 @@ a = 1.0 # in nm
gain = 5 gain = 5
b_0 = 1 b_0 = 1
inter = 'none' inter = 'none'
dim = (1,) + (64, 64)
dim_small = (64, 64)
smoothed_pictures = True smoothed_pictures = True
lam = 1E-6 lam = 1E-6
log = True log = True
PATH = '../../output/joern/' PATH = '../../output/joern/'
################################################################################################### ###################################################################################################
# Read in files: # Read in files:
phase_map = PhaseMap.load_from_netcdf4(PATH+'phase_map.nc') phase_map = PhaseMap.load_from_netcdf4(PATH+'phase_map_2.nc')
#mask = np.genfromtxt('mask.txt', dtype=bool) with open(PATH + 'mask_2.pickle') as pf:
with open(PATH + 'mask.pickle') as pf:
mask = pickle.load(pf) mask = pickle.load(pf)
dim = mask.shape
# Setup: # Setup:
if not use_mask: if not use_mask:
mask = np.ones_like(mask, dtype=bool) mask = np.ones_like(mask, dtype=bool)
...@@ -109,8 +107,8 @@ plt.savefig(dirname + "/reconstr.png") ...@@ -109,8 +107,8 @@ plt.savefig(dirname + "/reconstr.png")
# Plot the magnetization: # Plot the magnetization:
axis = (mag_data_rec*(1/mag_data_rec.magnitude.max())).quiver_plot(show=False) axis = (mag_data_rec*(1/mag_data_rec.magnitude.max())).quiver_plot(show=False)
axis.set_xlim(20, 45) axis.set_xlim(int(20/64*dim[1], 45/64*dim[2])
axis.set_ylim(20, 45) axis.set_ylim(int(20/64*dim[1], 45/64*dim[2])
plt.savefig(dirname + "/quiver.png") plt.savefig(dirname + "/quiver.png")
# Display the Phase: # Display the Phase:
...@@ -125,14 +123,14 @@ axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1, ...@@ -125,14 +123,14 @@ axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1,
interpolation=inter, show=False) interpolation=inter, show=False)
mag_data_rec.quiver_plot(axis=axis, show=False) mag_data_rec.quiver_plot(axis=axis, show=False)
axis = plt.gca() axis = plt.gca()
axis.set_xlim(20, 45) axis.set_xlim(int(20/64*dim[1], 45/64*dim[2])
axis.set_ylim(20, 45) axis.set_ylim(int(20/64*dim[1], 45/64*dim[2])
plt.savefig(dirname + "/overlay_normal.png") plt.savefig(dirname + "/overlay_normal.png")
axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1, axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1,
interpolation=inter, show=False) interpolation=inter, show=False)
mag_data_rec.quiver_plot(axis=axis, log=log, show=False) mag_data_rec.quiver_plot(axis=axis, log=log, show=False)
axis = plt.gca() axis = plt.gca()
axis.set_xlim(20, 45) axis.set_xlim(int(20/64*dim[1], 45/64*dim[2])
axis.set_ylim(20, 45) axis.set_ylim(int(20/64*dim[1], 45/64*dim[2])
plt.savefig(dirname + "/overlay_log.png") plt.savefig(dirname + "/overlay_log.png")
...@@ -9,6 +9,8 @@ import os ...@@ -9,6 +9,8 @@ import os
import numpy as np import numpy as np
import pickle
from pyDM3reader import DM3lib as dm3 from pyDM3reader import DM3lib as dm3
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
...@@ -34,7 +36,7 @@ a = 1.0 # in nm ...@@ -34,7 +36,7 @@ a = 1.0 # in nm
gain = 5 gain = 5
b_0 = 1 b_0 = 1
inter = 'none' inter = 'none'
dim_small = (64, 64) dim_small = (512, 512)
smoothed_pictures = True smoothed_pictures = True
lam = 1E-4 lam = 1E-4
order = 1 order = 1
...@@ -65,61 +67,68 @@ phase_map_4.display_combined(gain=gain, interpolation=inter) ...@@ -65,61 +67,68 @@ phase_map_4.display_combined(gain=gain, interpolation=inter)
plt.savefig(PATH+'%.0e/phase_map_4part.png'%lam) 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, mask_4 = np.expand_dims(np.where(np.array(dm3_4_ele.resize(dim_small)) >= threshold,
True, False), axis=0) True, False), axis=0)
# Reconstruct the magnetic distribution: phase_map_2.save_to_netcdf4('../../output/joern/phase_map_2.nc')
tic = clock() phase_map_4.save_to_netcdf4('../../output/joern/phase_map_4.nc')
mag_data_rec_2 = rc.optimize_simple_leastsq(phase_map_2, mask_2, b_0, lam=lam, order=order) with open('../../output/joern/mask_2.pickle', 'wb') as pf:
print '2 particle reconstruction time:', clock() - tic pickle.dump(mask_2, pf)
tic = clock() with open('../../output/joern/mask_4.pickle', 'wb') as pf:
mag_data_rec_4 = rc.optimize_simple_leastsq(phase_map_4, mask_4, b_0, lam=lam, order=order) pickle.dump(mask_4, pf)
print '4 particle reconstruction time:', clock() - tic
# Display the reconstructed phase map and holography image: ## Reconstruct the magnetic distribution:
phase_map_rec_2 = pm(mag_data_rec_2) #tic = clock()
phase_map_rec_2.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter) #mag_data_rec_2 = rc.optimize_simple_leastsq(phase_map_2, mask_2, b_0, lam=lam, order=order)
plt.savefig(PATH+'%.0e/phase_map_2part_rec.png'%lam) #print '2 particle reconstruction time:', clock() - tic
phase_map_rec_4 = pm(mag_data_rec_4) #tic = clock()
phase_map_rec_4.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter) #mag_data_rec_4 = rc.optimize_simple_leastsq(phase_map_4, mask_4, b_0, lam=lam, order=order)
plt.savefig(PATH+'%.0e/phase_map_4part_rec.png'%lam) #print '4 particle reconstruction time:', clock() - tic
# Plot the magnetization: ## Display the reconstructed phase map and holography image:
axis = (mag_data_rec_2*(1/mag_data_rec_2.magnitude.max())).quiver_plot() #phase_map_rec_2 = pm(mag_data_rec_2)
axis.set_xlim(20, 45) #phase_map_rec_2.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter)
axis.set_ylim(20, 45) #plt.savefig(PATH+'%.0e/phase_map_2part_rec.png'%lam)
plt.savefig(PATH+'%.0e/mag_data_2part.png'%lam) #phase_map_rec_4 = pm(mag_data_rec_4)
axis = (mag_data_rec_4*(1/mag_data_rec_4.magnitude.max())).quiver_plot() #phase_map_rec_4.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter)
axis.set_xlim(20, 45) #plt.savefig(PATH+'%.0e/phase_map_4part_rec.png'%lam)
axis.set_ylim(20, 45) ## Plot the magnetization:
plt.savefig(PATH+'%.0e/mag_data_4part.png'%lam) #axis = (mag_data_rec_2*(1/mag_data_rec_2.magnitude.max())).quiver_plot()
# Display the Phase: #axis.set_xlim(20, 45)
phase_diff_2 = phase_map_rec_2-phase_map_2 #axis.set_ylim(20, 45)
phase_diff_2.display_phase('Difference') #plt.savefig(PATH+'%.0e/mag_data_2part.png'%lam)
plt.savefig(PATH+'%.0e/phase_map_2part_diff.png'%lam) #axis = (mag_data_rec_4*(1/mag_data_rec_4.magnitude.max())).quiver_plot()
phase_diff_4 = phase_map_rec_4-phase_map_4 #axis.set_xlim(20, 45)
phase_diff_4.display_phase('Difference') #axis.set_ylim(20, 45)
plt.savefig(PATH+'%.0e/phase_map_4part_diff.png'%lam) #plt.savefig(PATH+'%.0e/mag_data_4part.png'%lam)
# Get the average difference from the experimental results: ## Display the Phase:
print 'Average difference (2 cubes):', np.average(phase_diff_2.phase) #phase_diff_2 = phase_map_rec_2-phase_map_2
print 'Average difference (4 cubes):', np.average(phase_diff_4.phase) #phase_diff_2.display_phase('Difference')
# Plot holographic contour maps with overlayed magnetic distributions: #plt.savefig(PATH+'%.0e/phase_map_2part_diff.png'%lam)
axis = phase_map_rec_2.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter) #phase_diff_4 = phase_map_rec_4-phase_map_4
mag_data_rec_2.quiver_plot(axis=axis) #phase_diff_4.display_phase('Difference')
axis = plt.gca() #plt.savefig(PATH+'%.0e/phase_map_4part_diff.png'%lam)
axis.set_xlim(20, 45) ## Get the average difference from the experimental results:
axis.set_ylim(20, 45) #print 'Average difference (2 cubes):', np.average(phase_diff_2.phase)
plt.savefig(PATH+'%.0e/phase_map_2part_holo.png'%lam) #print 'Average difference (4 cubes):', np.average(phase_diff_4.phase)
axis = phase_map_rec_4.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter) ## Plot holographic contour maps with overlayed magnetic distributions:
mag_data_rec_4.quiver_plot(axis=axis) #axis = phase_map_rec_2.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
axis = plt.gca() #mag_data_rec_2.quiver_plot(axis=axis)
axis.set_xlim(20, 45) #axis = plt.gca()
axis.set_ylim(20, 45) #axis.set_xlim(20, 45)
plt.savefig(PATH+'%.0e/phase_map_4part_holo.png'%lam) #axis.set_ylim(20, 45)
axis = phase_map_rec_2.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter) #plt.savefig(PATH+'%.0e/phase_map_2part_holo.png'%lam)
mag_data_rec_2.quiver_plot(axis=axis, log=log) #axis = phase_map_rec_4.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
axis = plt.gca() #mag_data_rec_4.quiver_plot(axis=axis)
axis.set_xlim(20, 45) #axis = plt.gca()
axis.set_ylim(20, 45) #axis.set_xlim(20, 45)
plt.savefig(PATH+'%.0e/phase_map_2part_holo_log.png'%lam) #axis.set_ylim(20, 45)
axis = phase_map_rec_4.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter) #plt.savefig(PATH+'%.0e/phase_map_4part_holo.png'%lam)
mag_data_rec_4.quiver_plot(axis=axis, log=log) #axis = phase_map_rec_2.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
axis = plt.gca() #mag_data_rec_2.quiver_plot(axis=axis, log=log)
axis.set_xlim(20, 45) #axis = plt.gca()
axis.set_ylim(20, 45) #axis.set_xlim(20, 45)
plt.savefig(PATH+'%.0e/phase_map_4part_holo_log.png'%lam) #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)
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