Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • empyre/empyre
  • weber/empyre
  • wessels/empyre
  • bryan/empyre
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1504 additions and 2158 deletions
scripts/publications/paper 1/ch5-0-evaluation_and_comparison.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Created on Fri Jul 26 14:37:20 2013 @author: Jan"""
import os
import numpy as np
from numpy import pi
import matplotlib.pyplot as plt
from matplotlib.ticker import FixedFormatter, IndexLocator
import pyramid
import pyramid.magcreator as mc
from pyramid.magdata import MagData
import shelve
import logging
import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
force_calculation = False
print '\nACCESS SHELVE'
# Create / Open databank:
directory = '../../../output/paper 1'
if not os.path.exists(directory):
os.makedirs(directory)
data_shelve = shelve.open(directory + '/paper_1_shelve')
###############################################################################################
print 'CH5-0 KERNELS'
x = np.linspace(-5, 5, 1000)
y_r = x/np.abs(x)**3
y_k = x/np.abs(x)**2
fig = plt.figure()
axis = fig.add_subplot(1, 1, 1, aspect='equal')
axis.plot(x, y_r, 'r', label=r'$r/|r|^3$')
axis.plot(x, y_k, 'b', label=r'$k/|k|^2$')
axis.set_xlim(-5, 5)
axis.set_ylim(-5, 5)
axis.axvline(0, linewidth=2, color='k')
axis.axhline(0, linewidth=2, color='k')
axis.legend()
###############################################################################################
print 'CH5-0 MAGNETIC DISTRIBUTIONS'
key = 'ch5-0-magnetic_distributions'
if key in data_shelve and not force_calculation:
print '--LOAD MAGNETIC DISTRIBUTIONS'
(mag_data_disc, mag_data_vort) = data_shelve[key]
else:
print '--CREATE MAGNETIC DISTRIBUTIONS'
# Input parameters:
a = 1.0 # in nm
phi = pi/2
dim = (16, 128, 128) # in px (z, y, x)
# Create magnetic shape:
center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts at 0!
radius = dim[1]/4 # in px
height = dim[0]/2 # in px
mag_shape = mc.Shapes.disc(dim, center, radius, height)
print '--CREATE MAGN. DISTR. OF HOMOG. MAG. DISC'
mag_data_disc = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
mag_data_disc.scale_down(2)
print '--CREATE MAGN. DISTR. OF VORTEX STATE DISC'
mag_data_vort = MagData(a, mc.create_mag_dist_vortex(mag_shape, center))
mag_data_vort.scale_down(2)
# Mayavi-Plots:
mag_data_disc.quiver_plot3d()
mag_data_vort.quiver_plot3d()
print '--SHELVE MAGNETIC DISTRIBUTIONS'
data_shelve[key] = (mag_data_disc, mag_data_vort)
print '--PLOT/SAVE MAGNETIC DISTRIBUTIONS'
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Magnetic Distributions', fontsize=20)
# Plot MagData (Disc):
mag_data_disc.quiver_plot('Homog. magn. disc', axis=axes[0])
axes[0].set_aspect('equal')
axes[0].xaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
axes[0].yaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
axes[0].xaxis.set_major_formatter(FixedFormatter([16*i for i in range(10)]))
axes[0].yaxis.set_major_formatter(FixedFormatter([16*i for i in range(10)]))
axes[0].set_xlabel('x [nm]', fontsize=15)
axes[0].set_ylabel('x [nm]', fontsize=15)
# Plot MagData (Disc):
mag_data_vort.quiver_plot('Vortex state disc', axis=axes[1])
axes[1].set_aspect('equal')
axes[1].xaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
axes[1].yaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
axes[1].xaxis.set_major_formatter(FixedFormatter([16*i for i in range(10)]))
axes[1].yaxis.set_major_formatter(FixedFormatter([16*i for i in range(10)]))
axes[1].set_xlabel('x [nm]', fontsize=15)
axes[1].set_ylabel('x [nm]', fontsize=15)
# Save Plots:
plt.figtext(0.15, 0.15, 'a)', fontsize=30)
plt.figtext(0.57, 0.15, 'b)', fontsize=30)
plt.savefig(directory + '/ch5-0-magnetic_distributions.png', bbox_inches='tight')
###############################################################################################
print 'CLOSING SHELVE\n'
# Close shelve:
data_shelve.close()
scripts/publications/paper 1/ch5-1-evaluation_real_space.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Created on Fri Jul 26 14:37:30 2013 @author: Jan"""
import os
import numpy as np
from numpy import pi
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm
from matplotlib.ticker import MaxNLocator
from matplotlib.cm import RdBu
from matplotlib.patches import Rectangle
import pyramid
import pyramid.magcreator as mc
import pyramid.analytic as an
from pyramid.phasemapper import pm
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
import shelve
import logging
import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
PHI_0 = -2067.83 # magnetic flux in T*nm²
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
force_calculation = False
print '\nACCESS SHELVE'
# Create / Open databank:
directory = '../../../output/paper 1'
if not os.path.exists(directory):
os.makedirs(directory)
data_shelve = shelve.open(directory + '/paper_1_shelve')
###############################################################################################
print 'CH5-1 ANALYTIC SOLUTIONS'
# Input parameters:
a = 0.125 # in nm
phi = pi/2
dim = (128, 1024, 1024) # in px (z, y, x)
gain = 100
# Create magnetic shape:
center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts with 0!
radius = dim[1]/4 # in px
height = dim[0]/2 # in px
print '--CALCULATE ANALYTIC SOLUTIONS'
# Get analytic solution:
phase_map_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
phase_map_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
print '--PLOT/SAVE ANALYTIC SOLUTIONS'
phase_map_ana_disc.display_combined('Analytic solution: hom. magn. disc', 'bilinear', gain=gain)
axis = plt.gcf().add_subplot(1, 2, 2, aspect='equal')
axis.axhline(y=512, linewidth=3, linestyle='--', color='r')
plt.figtext(0.15, 0.2, 'a)', fontsize=30, color='w')
#plt.figtext(0.52, 0.2, 'b)', fontsize=30)
plt.savefig(directory + '/ch5-1-analytic_solution_disc.png', bbox_inches='tight')
phase_map_ana_vort.display_combined('Analytic solution: vortex state', 'bilinear', gain=gain)
axis = plt.gcf().add_subplot(1, 2, 2, aspect='equal')
axis.axhline(y=512, linewidth=3, linestyle='--', color='r')
plt.figtext(0.15, 0.2, 'c)', fontsize=30, color='w')
plt.figtext(0.52, 0.2, 'd)', fontsize=30)
plt.savefig(directory + '/ch5-1-analytic_solution_vort.png', bbox_inches='tight')
# Get colorwheel:
PhaseMap.make_color_wheel()
plt.figtext(0.15, 0.14, 'e)', fontsize=30, color='w')
plt.savefig(directory + '/ch5-1-colorwheel.png', bbox_inches='tight')
###############################################################################################
print 'CH5-1 PHASE SLICES REAL SPACE'
# Input parameters:
a = 0.5 # in nm
phi = pi/2
gain = 100
dim = (32, 256, 256) # in px (z, y, x)
# Create magnetic shape:
center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts with 0!
radius = dim[1]/4 # in px
height = dim[0]/2 # in px
key = 'ch5-1-phase_slices_real'
if key in data_shelve and not force_calculation:
print '--LOAD MAGNETIC DISTRIBUTION'
(x_d, y_d, dy_d, x_v, y_v, dy_v) = data_shelve[key]
else:
print '--CREATE MAGNETIC DISTRIBUTION'
mag_shape = mc.Shapes.disc(dim, center, radius, height)
print '--CREATE PHASE SLICES HOMOG. MAGN. DISC'
# Arrays for plotting:
x_d = []
y_d = []
dy_d = []
# Analytic solution:
L = dim[1] * a # in px/nm
Lz = 0.5 * dim[0] * a # in px/nm
R = 0.25 * L # in px/nm
x0 = L / 2 # in px/nm
def F_disc(x):
coeff = - pi * Lz / (2*PHI_0) * 1E3 # in mrad -> *1000
result = coeff * (- (x - x0) * np.sin(phi))
result *= np.where(np.abs(x - x0) <= R, 1, (R / (x - x0)) ** 2)
return result
x_d.append(np.linspace(0, L, 5000))
y_d.append(F_disc(x_d[0]))
dy_d.append(np.zeros_like(x_d[0]))
# Create MagData (Disc):
mag_data_disc = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
for i in range(5):
print '----a =', mag_data_disc.a, 'nm', 'dim =', mag_data_disc.dim
phase_map = pm(mag_data_disc)
phase_map.display_combined('Disc, a = {} nm'.format(mag_data_disc.a), gain=gain)
x_d.append(np.linspace(mag_data_disc.a * 0.5,
mag_data_disc.a * (mag_data_disc.dim[1]-0.5),
mag_data_disc.dim[1]))
slice_pos = int(mag_data_disc.dim[1]/2)
y_d.append(phase_map.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
dy_d.append(phase_map.phase[slice_pos, :]*1E3 - F_disc(x_d[-1])) # *1E3: rad to mrad
if i < 4:
mag_data_disc.scale_down()
print '--CREATE PHASE SLICES VORTEX STATE DISC'
x_v = []
y_v = []
dy_v = []
# Analytic solution:
L = dim[1] * a # in px/nm
Lz = 0.5 * dim[0] * a # in px/nm
R = 0.25 * L # in px/nm
x0 = L / 2 # in px/nm
def F_vort(x):
coeff = pi*Lz/PHI_0 * 1E3 # in mrad -> *1000
result = coeff * np.where(np.abs(x - x0) <= R, (np.abs(x-x0)-R), 0)
return result
x_v.append(np.linspace(0, L, 5001))
y_v.append(F_vort(x_v[0]))
dy_v.append(np.zeros_like(x_v[0]))
# Create MagData (Vortex):
mag_data_vort = MagData(a, mc.create_mag_dist_vortex(mag_shape))
for i in range(5):
print '----a =', mag_data_vort.a, 'nm', 'dim =', mag_data_vort.dim
phase_map = pm(mag_data_vort)
phase_map.display_combined('Disc, a = {} nm'.format(mag_data_vort.a), gain=gain)
x_v.append(np.linspace(mag_data_vort.a * 0.5,
mag_data_vort.a * (mag_data_vort.dim[1]-0.5),
mag_data_vort.dim[1]))
slice_pos = int(mag_data_vort.dim[1]/2)
y_v.append(phase_map.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
dy_v.append(phase_map.phase[slice_pos, :]*1E3 - F_vort(x_v[-1])) # *1E3: rad to mrad
if i < 4:
mag_data_vort.scale_down()
# Shelve x, y and dy:
print '--SAVE PHASE SLICES'
data_shelve[key] = (x_d, y_d, dy_d, x_v, y_v, dy_v)
# Create figure:
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Central phase slices', fontsize=20)
print '--PLOT/SAVE PHASE SLICES HOMOG. MAGN. DISC'
# Plot phase slices:
axes[0].plot(x_d[0], y_d[0], '-k', linewidth=1.5, label='analytic')
axes[0].plot(x_d[1], y_d[1], '-r', linewidth=1.5, label='0.5 nm')
axes[0].plot(x_d[2], y_d[2], '-m', linewidth=1.5, label='1 nm')
axes[0].plot(x_d[3], y_d[3], '-y', linewidth=1.5, label='2 nm')
axes[0].plot(x_d[4], y_d[4], '-g', linewidth=1.5, label='4 nm')
axes[0].plot(x_d[5], y_d[5], '-c', linewidth=1.5, label='8 nm')
axes[0].tick_params(axis='both', which='major', labelsize=14)
axes[0].set_title('Homog. magn. disc', fontsize=18)
axes[0].set_xlabel('x [nm]', fontsize=15)
axes[0].set_ylabel('phase [mrad]', fontsize=15)
axes[0].set_xlim(0, 128)
axes[0].set_ylim(-220, 220)
# Plot Zoombox and Arrow:
zoom = (23.5, 160, 15, 40)
rect = Rectangle((zoom[0], zoom[1]), zoom[2], zoom[3], fc='w', ec='k')
axes[0].add_patch(rect)
axes[0].arrow(zoom[0]+zoom[2], zoom[1]+zoom[3]/2, 36, 0, length_includes_head=True,
head_width=10, head_length=4, fc='k', ec='k')
# Plot zoom inset:
ins_axis_d = plt.axes([0.33, 0.57, 0.14, 0.3])
ins_axis_d.plot(x_d[0], y_d[0], '-k', linewidth=1.5, label='analytic')
ins_axis_d.plot(x_d[1], y_d[1], '-r', linewidth=1.5, label='0.5 nm')
ins_axis_d.plot(x_d[2], y_d[2], '-m', linewidth=1.5, label='1 nm')
ins_axis_d.plot(x_d[3], y_d[3], '-y', linewidth=1.5, label='2 nm')
ins_axis_d.plot(x_d[4], y_d[4], '-g', linewidth=1.5, label='4 nm')
ins_axis_d.plot(x_d[5], y_d[5], '-c', linewidth=1.5, label='8 nm')
ins_axis_d.tick_params(axis='both', which='major', labelsize=14)
ins_axis_d.set_xlim(zoom[0], zoom[0]+zoom[2])
ins_axis_d.set_ylim(zoom[1], zoom[1]+zoom[3])
ins_axis_d.xaxis.set_major_locator(MaxNLocator(nbins=4, integer=True))
ins_axis_d.yaxis.set_major_locator(MaxNLocator(nbins=3))
print '--PLOT/SAVE PHASE SLICES VORTEX STATE DISC'
# Plot phase slices:
axes[1].plot(x_v[0], y_v[0], '-k', linewidth=1.5, label='analytic')
axes[1].plot(x_v[1], y_v[1], '-r', linewidth=1.5, label='0.5 nm')
axes[1].plot(x_v[2], y_v[2], '-m', linewidth=1.5, label='1 nm')
axes[1].plot(x_v[3], y_v[3], '-y', linewidth=1.5, label='2 nm')
axes[1].plot(x_v[4], y_v[4], '-g', linewidth=1.5, label='4 nm')
axes[1].plot(x_v[5], y_v[5], '-c', linewidth=1.5, label='8 nm')
axes[1].tick_params(axis='both', which='major', labelsize=14)
axes[1].set_title('Vortex state disc', fontsize=18)
axes[1].set_xlabel('x [nm]', fontsize=15)
axes[1].set_ylabel('phase [mrad]', fontsize=15)
axes[1].set_xlim(0, 128)
axes[1].yaxis.set_major_locator(MaxNLocator(nbins=6))
axes[1].legend()
# Plot Zoombox and Arrow:
zoom = (59, 340, 10, 55)
rect = Rectangle((zoom[0], zoom[1]), zoom[2], zoom[3], fc='w', ec='k')
axes[1].add_patch(rect)
axes[1].arrow(zoom[0]+zoom[2]/2, zoom[1], 0, -193, length_includes_head=True,
head_width=2, head_length=20, fc='k', ec='k')
# Plot zoom inset:
ins_axis_v = plt.axes([0.695, 0.15, 0.075, 0.3])
ins_axis_v.plot(x_v[0], y_v[0], '-k', linewidth=1.5, label='analytic')
ins_axis_v.plot(x_v[1], y_v[1], '-r', linewidth=1.5, label='0.5 nm')
ins_axis_v.plot(x_v[2], y_v[2], '-m', linewidth=1.5, label='1 nm')
ins_axis_v.plot(x_v[3], y_v[3], '-y', linewidth=1.5, label='2 nm')
ins_axis_v.plot(x_v[4], y_v[4], '-g', linewidth=1.5, label='4 nm')
ins_axis_v.plot(x_v[5], y_v[5], '-c', linewidth=1.5, label='8 nm')
ins_axis_v.tick_params(axis='both', which='major', labelsize=14)
ins_axis_v.set_xlim(zoom[0], zoom[0]+zoom[2])
ins_axis_v.set_ylim(zoom[1], zoom[1]+zoom[3])
ins_axis_v.xaxis.set_major_locator(MaxNLocator(nbins=4, integer=True))
ins_axis_v.yaxis.set_major_locator(MaxNLocator(nbins=4))
plt.show()
#plt.figtext(0.15, 0.13, 'a)', fontsize=30)
#plt.figtext(0.57, 0.13, 'b)', fontsize=30)
plt.savefig(directory + '/ch5-1-phase_slice_comparison.png', bbox_inches='tight')
# Create figure:
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Central phase slice errors', fontsize=20)
print '--PLOT/SAVE PHASE SLICE ERRORS HOMOG. MAGN. DISC'
# Plot phase slices:
axes[0].plot(x_d[0], dy_d[0], '-k', linewidth=1.5, label='analytic')
axes[0].plot(x_d[1], dy_d[1], '-r', linewidth=1.5, label='0.5 nm')
axes[0].plot(x_d[2], dy_d[2], '-m', linewidth=1.5, label='1 nm')
axes[0].plot(x_d[3], dy_d[3], '-y', linewidth=1.5, label='2 nm')
axes[0].plot(x_d[4], dy_d[4], '-g', linewidth=1.5, label='4 nm')
axes[0].plot(x_d[5], dy_d[5], '-c', linewidth=1.5, label='8 nm')
axes[0].tick_params(axis='both', which='major', labelsize=14)
axes[0].set_title('Homog. magn. disc', fontsize=18)
axes[0].set_xlabel('x [nm]', fontsize=15)
axes[0].set_ylabel('phase [mrad]', fontsize=15)
axes[0].set_xlim(0, 128)
print '--PLOT/SAVE PHASE SLICE ERRORS VORTEX STATE DISC'
# Plot phase slices:
axes[1].plot(x_v[0], dy_v[0], '-k', linewidth=1.5, label='analytic')
axes[1].plot(x_v[1], dy_v[1], '-r', linewidth=1.5, label='0.5 nm')
axes[1].plot(x_v[2], dy_v[2], '-m', linewidth=1.5, label='1 nm')
axes[1].plot(x_v[3], dy_v[3], '-y', linewidth=1.5, label='2 nm')
axes[1].plot(x_v[4], dy_v[4], '-g', linewidth=1.5, label='4 nm')
axes[1].plot(x_v[5], dy_v[5], '-c', linewidth=1.5, label='8 nm')
axes[1].tick_params(axis='both', which='major', labelsize=14)
axes[1].set_title('Vortex state disc', fontsize=18)
axes[1].set_xlabel('x [nm]', fontsize=15)
axes[1].set_ylabel('phase [mrad]', fontsize=15)
axes[1].set_xlim(0, 128)
axes[1].legend(loc=4)
plt.show()
#plt.figtext(0.15, 0.13, 'a)', fontsize=30)
#plt.figtext(0.57, 0.13, 'b)', fontsize=30)
plt.savefig(directory + '/ch5-1-phase_slice_errors.png', bbox_inches='tight')
###############################################################################################
print 'CH5-1 PHASE DIFFERENCES REAL SPACE'
# Input parameters:
a = 1.0 # in nm
phi = pi/2
dim = (16, 128, 128) # in px (z, y, x)
center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts with 0!
radius = dim[1]/4 # in px
height = dim[0]/2 # in px
key = 'ch5-1-phase_diff_mag_dist'
if key in data_shelve and not force_calculation:
print '--LOAD MAGNETIC DISTRIBUTIONS'
(mag_data_disc, mag_data_vort) = data_shelve[key]
else:
print '--CREATE MAGNETIC DISTRIBUTIONS'
# Create magnetic shape (4 times the size):
a_big = a / 2
dim_big = (dim[0]*2, dim[1]*2, dim[2]*2)
center_big = (dim_big[0]/2-0.5, dim_big[1]/2.-0.5, dim_big[2]/2.-0.5)
radius_big = dim_big[1]/4 # in px
height_big = dim_big[0]/2 # in px
mag_shape = mc.Shapes.disc(dim_big, center_big, radius_big, height_big)
# Create MagData (4 times the size):
mag_data_disc = MagData(a_big, mc.create_mag_dist_homog(mag_shape, phi))
mag_data_vort = MagData(a_big, mc.create_mag_dist_vortex(mag_shape, center_big))
# Scale mag_data, grid spacing and dimensions:
mag_data_disc.scale_down()
mag_data_vort.scale_down()
print '--SAVE MAGNETIC DISTRIBUTIONS'
# Shelve magnetic distributions:
data_shelve[key] = (mag_data_disc, mag_data_vort)
print '--CALCULATE PHASE DIFFERENCES'
# Get analytic solutions:
phase_map_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
phase_map_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
# Create norm for the plots:
bounds = np.array([-3, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 3])
norm = BoundaryNorm(bounds, RdBu.N)
# Calculations (Disc):
phase_map_num_disc = pm(mag_data_disc)
phase_map_num_disc.unit = 'mrad'
phase_diff_disc = phase_map_num_disc - phase_map_ana_disc
RMS_disc = np.sqrt(np.mean(phase_diff_disc.phase**2))
# Calculations (Vortex):
phase_map_num_vort = pm(mag_data_vort)
phase_map_num_vort.unit = 'mrad'
phase_diff_vort = phase_map_num_vort - phase_map_ana_vort
RMS_vort = np.sqrt(np.mean(phase_diff_vort.phase**2))
print '--PLOT/SAVE PHASE DIFFERENCES'
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Difference of the real space approach from the analytical solution', fontsize=20)
# Plot MagData (Disc):
phase_diff_disc.display_phase(u'Homog. magn. disc, RMS = {:3.2f} µrad'.format(RMS_disc*1000),
limit=np.max(bounds), norm=norm, axis=axes[0])
axes[0].set_aspect('equal')
# Plot MagData (Disc):
phase_diff_vort.display_phase(u'Vortex state disc, RMS = {:3.2f} µrad'.format(RMS_vort*1000),
limit=np.max(bounds), norm=norm, axis=axes[1])
axes[1].set_aspect('equal')
# Save Plots:
plt.figtext(0.15, 0.2, 'a)', fontsize=30)
plt.figtext(0.52, 0.2, 'b)', fontsize=30)
plt.savefig(directory + '/ch5-1-phase_differences.png', bbox_inches='tight')
###############################################################################################
print 'CLOSING SHELVE\n'
# Close shelve:
data_shelve.close()
scripts/publications/paper 1/ch5-2-evaluation_fourier_space.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Created on Fri Jul 26 14:37:30 2013 @author: Jan"""
import os
import numpy as np
from numpy import pi
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm
from matplotlib.ticker import MaxNLocator
from matplotlib.cm import RdBu
import pyramid
import pyramid.magcreator as mc
import pyramid.analytic as an
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector
from pyramid.phasemapper import PhaseMapperFDFC
import time
import shelve
import logging
import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
PHI_0 = -2067.83 # magnetic flux in T*nm²
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
force_calculation = False
print '\nACCESS SHELVE'
# Create / Open databank:
directory = '../../../output/paper 1'
if not os.path.exists(directory):
os.makedirs(directory)
data_shelve = shelve.open(directory + '/paper_1_shelve')
###############################################################################################
print 'CH5-1 PHASE SLICES FOURIER SPACE'
# Input parameters:
a = 0.5 # in nm
phi = pi/2
gain = 100
dim = (32, 256, 256) # in px (z, y, x)
# Create magnetic shape:
center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts with 0!
radius = dim[1]/4 # in px
height = dim[0]/2 # in px
key = 'ch5-2-phase_slices_fourier'
if key in data_shelve and not force_calculation:
print '--LOAD MAGNETIC DISTRIBUTION'
(x_d, y_d0, y_d10, dy_d0, dy_d10, x_v, y_v0, y_v10, dy_v0, dy_v10) = data_shelve[key]
else:
print '--CREATE MAGNETIC DISTRIBUTION'
mag_shape = mc.Shapes.disc(dim, center, radius, height)
print '--CREATE PHASE SLICES HOMOG. MAGN. DISC'
# Arrays for plotting:
x_d = []
y_d0 = []
y_d10 = []
dy_d0 = []
dy_d10 = []
# Analytic solution:
L = dim[1] * a # in px/nm
Lz = 0.5 * dim[0] * a # in px/nm
R = 0.25 * L # in px/nm
x0 = L / 2 # in px/nm
def F_disc(x):
coeff = - pi * Lz / (2*PHI_0) * 1E3 # in mrad -> *1000
result = coeff * (- (x - x0) * np.sin(phi))
result *= np.where(np.abs(x - x0) <= R, 1, (R / (x - x0)) ** 2)
return result
x_d.append(np.linspace(0, L, 5000))
y_d0.append(F_disc(x_d[0]))
y_d10.append(F_disc(x_d[0]))
dy_d0.append(np.zeros_like(x_d[0]))
dy_d10.append(np.zeros_like(x_d[0]))
# Create MagData (Disc):
mag_data_disc = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
for i in range(5):
print '----a =', mag_data_disc.a, 'nm', 'dim =', mag_data_disc.dim
projector = SimpleProjector(mag_data_disc.dim)
phase_map0 = PhaseMapperFDFC(mag_data_disc.a, projector.dim_uv,
padding=0)(projector(mag_data_disc))
phase_map10 = PhaseMapperFDFC(mag_data_disc.a, projector.dim_uv,
padding=10)(projector(mag_data_disc))
phase_map0.unit = 'mrad'
phase_map10.unit = 'mrad'
phase_map0.display_combined('Disc, a = {} nm'.format(mag_data_disc.a), gain=gain)
phase_map10.display_combined('Disc, a = {} nm'.format(mag_data_disc.a), gain=gain)
x_d.append(np.linspace(mag_data_disc.a * 0.5,
mag_data_disc.a * (mag_data_disc.dim[1]-0.5),
mag_data_disc.dim[1]))
slice_pos = int(mag_data_disc.dim[1]/2)
y_d0.append(phase_map0.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
y_d10.append(phase_map10.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
dy_d0.append(phase_map0.phase[slice_pos, :]*1E3 - F_disc(x_d[-1])) # *1E3: in mrad
dy_d10.append(phase_map10.phase[slice_pos, :]*1E3 - F_disc(x_d[-1])) # *1E3: in mrad
if i < 4:
mag_data_disc.scale_down()
print '--CREATE PHASE SLICES VORTEX STATE DISC'
x_v = []
y_v0 = []
y_v10 = []
dy_v0 = []
dy_v10 = []
# Analytic solution:
L = dim[1] * a # in px/nm
Lz = 0.5 * dim[0] * a # in px/nm
R = 0.25 * L # in px/nm
x0 = L / 2 # in px/nm
def F_vort(x):
coeff = pi*Lz/PHI_0 * 1E3 # in mrad -> *1000
result = coeff * np.where(np.abs(x - x0) <= R, (np.abs(x-x0)-R), 0)
return result
x_v.append(np.linspace(0, L, 5000))
y_v0.append(F_vort(x_v[0]))
y_v10.append(F_vort(x_v[0]))
dy_v0.append(np.zeros_like(x_v[0]))
dy_v10.append(np.zeros_like(x_v[0]))
# Create MagData (Vortex):
mag_data_vort = MagData(a, mc.create_mag_dist_vortex(mag_shape))
for i in range(5):
print '----a =', mag_data_vort.a, 'nm', 'dim =', mag_data_vort.dim
projector = SimpleProjector(mag_data_vort.dim)
phase_map0 = PhaseMapperFDFC(mag_data_vort.a, projector.dim_uv,
padding=0)(projector(mag_data_vort))
phase_map10 = PhaseMapperFDFC(mag_data_vort.a, projector.dim_uv,
padding=10)(projector(mag_data_vort))
phase_map0.unit = 'mrad'
phase_map10.unit = 'mrad'
print np.mean(phase_map0.phase)
phase_map0 -= phase_map0.phase[0, 0]
phase_map10 -= phase_map10.phase[0, 0]
phase_map0.display_combined('Vortex, a = {} nm'.format(mag_data_vort.a), gain=gain)
phase_map10.display_combined('Vortex, a = {} nm'.format(mag_data_vort.a), gain=gain)
x_v.append(np.linspace(mag_data_vort.a * 0.5,
mag_data_vort.a * (mag_data_vort.dim[1]-0.5),
mag_data_vort.dim[1]))
slice_pos = int(mag_data_vort.dim[1]/2)
y_v0.append(phase_map0.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
y_v10.append(phase_map10.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
dy_v0.append(phase_map0.phase[slice_pos, :]*1E3 - F_vort(x_v[-1])) # *1E3: in mrad
dy_v10.append(phase_map10.phase[slice_pos, :]*1E3 - F_vort(x_v[-1])) # *1E3: in mrad
if i < 4:
mag_data_vort.scale_down()
# Shelve x, y and dy:
print '--SAVE PHASE SLICES'
data_shelve[key] = (x_d, y_d0, y_d10, dy_d0, dy_d10, x_v, y_v0, y_v10, dy_v0, dy_v10)
# Create figure:
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Central phase slices (padding = 0)', fontsize=20)
print '--PLOT/SAVE PHASE SLICES HOMOG. MAGN. DISC PADDING = 0'
# Plot phase slices:
axes[0].plot(x_d[0], y_d0[0], '-k', linewidth=1.5, label='analytic')
axes[0].plot(x_d[1], y_d0[1], '-r', linewidth=1.5, label='0.5 nm')
axes[0].plot(x_d[2], y_d0[2], '-m', linewidth=1.5, label='1 nm')
axes[0].plot(x_d[3], y_d0[3], '-y', linewidth=1.5, label='2 nm')
axes[0].plot(x_d[4], y_d0[4], '-g', linewidth=1.5, label='4 nm')
axes[0].plot(x_d[5], y_d0[5], '-c', linewidth=1.5, label='8 nm')
axes[0].tick_params(axis='both', which='major', labelsize=14)
axes[0].set_title('Homog. magn. disc', fontsize=18)
axes[0].set_xlabel('x [nm]', fontsize=15)
axes[0].set_ylabel('phase [mrad]', fontsize=15)
axes[0].set_xlim(0, 128)
axes[0].set_ylim(-220, 220)
# # Plot Zoombox and Arrow:
# zoom = (23.5, 160, 15, 40)
# rect = Rectangle((zoom[0], zoom[1]), zoom[2], zoom[3], fc='w', ec='k')
# axes[0].add_patch(rect)
# axes[0].arrow(zoom[0]+zoom[2], zoom[1]+zoom[3]/2, 36, 0, length_includes_head=True,
# head_width=10, head_length=4, fc='k', ec='k')
# # Plot zoom inset:
# ins_axis_d = plt.axes([0.33, 0.57, 0.14, 0.3])
# ins_axis_d.plot(x_d[0], y_d0[0], '-k', linewidth=1.5, label='analytic')
# ins_axis_d.plot(x_d[1], y_d0[1], '-r', linewidth=1.5, label='0.5 nm')
# ins_axis_d.plot(x_d[2], y_d0[2], '-m', linewidth=1.5, label='1 nm')
# ins_axis_d.plot(x_d[3], y_d0[3], '-y', linewidth=1.5, label='2 nm')
# ins_axis_d.plot(x_d[4], y_d0[4], '-g', linewidth=1.5, label='4 nm')
# ins_axis_d.plot(x_d[5], y_d0[5], '-c', linewidth=1.5, label='8 nm')
# ins_axis_d.tick_params(axis='both', which='major', labelsize=14)
# ins_axis_d.set_xlim(zoom[0], zoom[0]+zoom[2])
# ins_axis_d.set_ylim(zoom[1], zoom[1]+zoom[3])
# ins_axis_d.xaxis.set_major_locator(MaxNLocator(nbins=4, integer= True))
# ins_axis_d.yaxis.set_major_locator(MaxNLocator(nbins=3))
print '--PLOT/SAVE PHASE SLICES VORTEX STATE DISC PADDING = 0'
# Plot phase slices:
axes[1].plot(x_v[0], y_v0[0], '-k', linewidth=1.5, label='analytic')
axes[1].plot(x_v[1], y_v0[1], '-r', linewidth=1.5, label='0.5 nm')
axes[1].plot(x_v[2], y_v0[2], '-m', linewidth=1.5, label='1 nm')
axes[1].plot(x_v[3], y_v0[3], '-y', linewidth=1.5, label='2 nm')
axes[1].plot(x_v[4], y_v0[4], '-g', linewidth=1.5, label='4 nm')
axes[1].plot(x_v[5], y_v0[5], '-c', linewidth=1.5, label='8 nm')
axes[1].tick_params(axis='both', which='major', labelsize=14)
axes[1].set_title('Vortex state disc', fontsize=18)
axes[1].set_xlabel('x [nm]', fontsize=15)
axes[1].set_ylabel('phase [mrad]', fontsize=15)
axes[1].set_xlim(0, 128)
axes[1].yaxis.set_major_locator(MaxNLocator(nbins=6))
axes[1].legend()
# # Plot Zoombox and Arrow:
# zoom = (59, 340, 10, 55)
# rect = Rectangle((zoom[0], zoom[1]), zoom[2], zoom[3], fc='w', ec='k')
# axes[1].add_patch(rect)
# axes[1].arrow(zoom[0]+zoom[2]/2, zoom[1], 0, -193, length_includes_head=True,
# head_width=2, head_length=20, fc='k', ec='k')
# # Plot zoom inset:
# ins_axis_v = plt.axes([0.695, 0.15, 0.075, 0.3])
# ins_axis_v.plot(x_v[0], y_v0[0], '-k', linewidth=1.5, label='analytic')
# ins_axis_v.plot(x_v[1], y_v0[1], '-r', linewidth=1.5, label='0.5 nm')
# ins_axis_v.plot(x_v[2], y_v0[2], '-m', linewidth=1.5, label='1 nm')
# ins_axis_v.plot(x_v[3], y_v0[3], '-y', linewidth=1.5, label='2 nm')
# ins_axis_v.plot(x_v[4], y_v0[4], '-g', linewidth=1.5, label='4 nm')
# ins_axis_v.plot(x_v[5], y_v0[5], '-c', linewidth=1.5, label='8 nm')
# ins_axis_v.tick_params(axis='both', which='major', labelsize=14)
# ins_axis_v.set_xlim(zoom[0], zoom[0]+zoom[2])
# ins_axis_v.set_ylim(zoom[1], zoom[1]+zoom[3])
# ins_axis_v.xaxis.set_major_locator(MaxNLocator(nbins=4, integer= True))
# ins_axis_v.yaxis.set_major_locator(MaxNLocator(nbins=4))
plt.show()
plt.figtext(0.15, 0.13, 'a)', fontsize=30)
plt.figtext(0.57, 0.13, 'b)', fontsize=30)
plt.savefig(directory + '/ch5-1-phase_slice_padding_0.png', bbox_inches='tight')
# Create figure:
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Central phase slices (padding = 10)', fontsize=20)
print '--PLOT/SAVE PHASE SLICES HOMOG. MAGN. DISC PADDING = 10'
# Plot phase slices:
axes[0].plot(x_d[0], y_d10[0], '-k', linewidth=1.5, label='analytic')
axes[0].plot(x_d[1], y_d10[1], '-r', linewidth=1.5, label='0.5 nm')
axes[0].plot(x_d[2], y_d10[2], '-m', linewidth=1.5, label='1 nm')
axes[0].plot(x_d[3], y_d10[3], '-y', linewidth=1.5, label='2 nm')
axes[0].plot(x_d[4], y_d10[4], '-g', linewidth=1.5, label='4 nm')
axes[0].plot(x_d[5], y_d10[5], '-c', linewidth=1.5, label='8 nm')
axes[0].tick_params(axis='both', which='major', labelsize=14)
axes[0].set_title('Homog. magn. disc', fontsize=18)
axes[0].set_xlabel('x [nm]', fontsize=15)
axes[0].set_ylabel('phase [mrad]', fontsize=15)
axes[0].set_xlim(0, 128)
axes[0].set_ylim(-220, 220)
# # Plot Zoombox and Arrow:
# zoom = (23.5, 160, 15, 40)
# rect = Rectangle((zoom[0], zoom[1]), zoom[2], zoom[3], fc='w', ec='k')
# axes[0].add_patch(rect)
# axes[0].arrow(zoom[0]+zoom[2], zoom[1]+zoom[3]/2, 36, 0, length_includes_head=True,
# head_width=10, head_length=4, fc='k', ec='k')
# # Plot zoom inset:
# ins_axis_d = plt.axes([0.33, 0.57, 0.14, 0.3])
# ins_axis_d.plot(x_d[0], y_d10[0], '-k', linewidth=1.5, label='analytic')
# ins_axis_d.plot(x_d[1], y_d10[1], '-r', linewidth=1.5, label='0.5 nm')
# ins_axis_d.plot(x_d[2], y_d10[2], '-m', linewidth=1.5, label='1 nm')
# ins_axis_d.plot(x_d[3], y_d10[3], '-y', linewidth=1.5, label='2 nm')
# ins_axis_d.plot(x_d[4], y_d10[4], '-g', linewidth=1.5, label='4 nm')
# ins_axis_d.plot(x_d[5], y_d10[5], '-c', linewidth=1.5, label='8 nm')
# ins_axis_d.tick_params(axis='both', which='major', labelsize=14)
# ins_axis_d.set_xlim(zoom[0], zoom[0]+zoom[2])
# ins_axis_d.set_ylim(zoom[1], zoom[1]+zoom[3])
# ins_axis_d.xaxis.set_major_locator(MaxNLocator(nbins=4, integer= True))
# ins_axis_d.yaxis.set_major_locator(MaxNLocator(nbins=3))
print '--PLOT/SAVE PHASE SLICES VORTEX STATE DISC PADDING = 10'
# Plot phase slices:
axes[1].plot(x_v[0], y_v10[0], '-k', linewidth=1.5, label='analytic')
axes[1].plot(x_v[1], y_v10[1], '-r', linewidth=1.5, label='0.5 nm')
axes[1].plot(x_v[2], y_v10[2], '-m', linewidth=1.5, label='1 nm')
axes[1].plot(x_v[3], y_v10[3], '-y', linewidth=1.5, label='2 nm')
axes[1].plot(x_v[4], y_v10[4], '-g', linewidth=1.5, label='4 nm')
axes[1].plot(x_v[5], y_v10[5], '-c', linewidth=1.5, label='8 nm')
axes[1].tick_params(axis='both', which='major', labelsize=14)
axes[1].set_title('Vortex state disc', fontsize=18)
axes[1].set_xlabel('x [nm]', fontsize=15)
axes[1].set_ylabel('phase [mrad]', fontsize=15)
axes[1].set_xlim(0, 128)
axes[1].yaxis.set_major_locator(MaxNLocator(nbins=6))
axes[1].legend()
# # Plot Zoombox and Arrow:
# zoom = (59, 340, 10, 55)
# rect = Rectangle((zoom[0], zoom[1]), zoom[2], zoom[3], fc='w', ec='k')
# axes[1].add_patch(rect)
# axes[1].arrow(zoom[0]+zoom[2]/2, zoom[1], 0, -193, length_includes_head=True,
# head_width=2, head_length=20, fc='k', ec='k')
# # Plot zoom inset:
# ins_axis_v = plt.axes([0.695, 0.15, 0.075, 0.3])
# ins_axis_v.plot(x_v[0], y_v10[0], '-k', linewidth=1.5, label='analytic')
# ins_axis_v.plot(x_v[1], y_v10[1], '-r', linewidth=1.5, label='0.5 nm')
# ins_axis_v.plot(x_v[2], y_v10[2], '-m', linewidth=1.5, label='1 nm')
# ins_axis_v.plot(x_v[3], y_v10[3], '-y', linewidth=1.5, label='2 nm')
# ins_axis_v.plot(x_v[4], y_v10[4], '-g', linewidth=1.5, label='4 nm')
# ins_axis_v.plot(x_v[5], y_v10[5], '-c', linewidth=1.5, label='8 nm')
# ins_axis_v.tick_params(axis='both', which='major', labelsize=14)
# ins_axis_v.set_xlim(zoom[0], zoom[0]+zoom[2])
# ins_axis_v.set_ylim(zoom[1], zoom[1]+zoom[3])
# ins_axis_v.xaxis.set_major_locator(MaxNLocator(nbins=4, integer= True))
# ins_axis_v.yaxis.set_major_locator(MaxNLocator(nbins=4))
plt.show()
plt.figtext(0.15, 0.13, 'c)', fontsize=30)
plt.figtext(0.57, 0.13, 'd)', fontsize=30)
plt.savefig(directory + '/ch5-1-phase_slice_padding_10.png', bbox_inches='tight')
# Create figure:
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Central phase slice errors (padding = 0)', fontsize=20)
print '--PLOT/SAVE PHASE SLICE ERRORS HOMOG. MAGN. DISC PADDING = 0'
# Plot phase slices:
axes[0].plot(x_d[0], dy_d0[0], '-k', linewidth=1.5, label='analytic')
axes[0].plot(x_d[1], dy_d0[1], '-r', linewidth=1.5, label='0.5 nm')
axes[0].plot(x_d[2], dy_d0[2], '-m', linewidth=1.5, label='1 nm')
axes[0].plot(x_d[3], dy_d0[3], '-y', linewidth=1.5, label='2 nm')
axes[0].plot(x_d[4], dy_d0[4], '-g', linewidth=1.5, label='4 nm')
axes[0].plot(x_d[5], dy_d0[5], '-c', linewidth=1.5, label='8 nm')
axes[0].tick_params(axis='both', which='major', labelsize=14)
axes[0].set_title('Homog. magn. disc', fontsize=18)
axes[0].set_xlabel('x [nm]', fontsize=15)
axes[0].set_ylabel('phase [mrad]', fontsize=15)
axes[0].set_xlim(0, 128)
print '--PLOT/SAVE PHASE SLICE ERRORS VORTEX STATE DISC PADDING = 0'
# Plot phase slices:
axes[1].plot(x_v[0], dy_v0[0], '-k', linewidth=1.5, label='analytic')
axes[1].plot(x_v[1], dy_v0[1], '-r', linewidth=1.5, label='0.5 nm')
axes[1].plot(x_v[2], dy_v0[2], '-m', linewidth=1.5, label='1 nm')
axes[1].plot(x_v[3], dy_v0[3], '-y', linewidth=1.5, label='2 nm')
axes[1].plot(x_v[4], dy_v0[4], '-g', linewidth=1.5, label='4 nm')
axes[1].plot(x_v[5], dy_v0[5], '-c', linewidth=1.5, label='8 nm')
axes[1].tick_params(axis='both', which='major', labelsize=14)
axes[1].set_title('Vortex state disc', fontsize=18)
axes[1].set_xlabel('x [nm]', fontsize=15)
axes[1].set_ylabel('phase [mrad]', fontsize=15)
axes[1].set_xlim(0, 128)
axes[1].legend(loc=4)
plt.show()
plt.figtext(0.15, 0.13, 'c)', fontsize=30)
plt.figtext(0.57, 0.13, 'd)', fontsize=30)
plt.savefig(directory + '/ch5-1-phase_slice_errors_padding_0.png', bbox_inches='tight')
# Create figure:
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Central phase slice errors (padding = 10)', fontsize=20)
print '--PLOT/SAVE PHASE SLICE ERRORS HOMOG. MAGN. DISC PADDING = 10'
# Plot phase slices:
axes[0].plot(x_d[0], dy_d10[0], '-k', linewidth=1.5, label='analytic')
axes[0].plot(x_d[1], dy_d10[1], '-r', linewidth=1.5, label='0.5 nm')
axes[0].plot(x_d[2], dy_d10[2], '-m', linewidth=1.5, label='1 nm')
axes[0].plot(x_d[3], dy_d10[3], '-y', linewidth=1.5, label='2 nm')
axes[0].plot(x_d[4], dy_d10[4], '-g', linewidth=1.5, label='4 nm')
axes[0].plot(x_d[5], dy_d10[5], '-c', linewidth=1.5, label='8 nm')
axes[0].tick_params(axis='both', which='major', labelsize=14)
axes[0].set_title('Homog. magn. disc', fontsize=18)
axes[0].set_xlabel('x [nm]', fontsize=15)
axes[0].set_ylabel('phase [mrad]', fontsize=15)
axes[0].set_xlim(0, 128)
print '--PLOT/SAVE PHASE SLICE ERRORS VORTEX STATE DISC PADDING = 10'
# Plot phase slices:
axes[1].plot(x_v[0], dy_v10[0], '-k', linewidth=1.5, label='analytic')
axes[1].plot(x_v[1], dy_v10[1], '-r', linewidth=1.5, label='0.5 nm')
axes[1].plot(x_v[2], dy_v10[2], '-m', linewidth=1.5, label='1 nm')
axes[1].plot(x_v[3], dy_v10[3], '-y', linewidth=1.5, label='2 nm')
axes[1].plot(x_v[4], dy_v10[4], '-g', linewidth=1.5, label='4 nm')
axes[1].plot(x_v[5], dy_v10[5], '-c', linewidth=1.5, label='8 nm')
axes[1].tick_params(axis='both', which='major', labelsize=14)
axes[1].set_title('Vortex state disc', fontsize=18)
axes[1].set_xlabel('x [nm]', fontsize=15)
axes[1].set_ylabel('phase [mrad]', fontsize=15)
axes[1].set_xlim(0, 128)
axes[1].legend(loc=4)
plt.show()
plt.figtext(0.15, 0.13, 'c)', fontsize=30)
plt.figtext(0.57, 0.13, 'd)', fontsize=30)
plt.savefig(directory + '/ch5-1-phase_slice_errors_padding_10.png', bbox_inches='tight')
###############################################################################################
print 'CH5-2 PHASE DIFFERENCES FOURIER SPACE'
# Input parameters:
a = 1.0 # in nm
phi = pi/2
dim = (16, 128, 128) # in px (z, y, x)
center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts with 0!
radius = dim[1]/4 # in px
height = dim[0]/2 # in px
key = 'ch5-1-phase_diff_mag_dist'
if key in data_shelve and not force_calculation:
print '--LOAD MAGNETIC DISTRIBUTIONS'
(mag_data_disc, mag_data_vort) = data_shelve[key]
else:
print '--CREATE MAGNETIC DISTRIBUTIONS'
# Create magnetic shape (4 times the size):
a_big = a / 2
dim_big = (dim[0]*2, dim[1]*2, dim[2]*2)
center_big = (dim_big[0]/2-0.5, dim_big[1]/2.-0.5, dim_big[2]/2.-0.5)
radius_big = dim_big[1]/4 # in px
height_big = dim_big[0]/2 # in px
mag_shape = mc.Shapes.disc(dim_big, center_big, radius_big, height_big)
# Create MagData (4 times the size):
mag_data_disc = MagData(a_big, mc.create_mag_dist_homog(mag_shape, phi))
mag_data_vort = MagData(a_big, mc.create_mag_dist_vortex(mag_shape, center_big))
# Scale mag_data, grid spacing and dimensions:
mag_data_disc.scale_down()
mag_data_vort.scale_down()
print '--SAVE MAGNETIC DISTRIBUTIONS'
# Shelve magnetic distributions:
data_shelve[key] = (mag_data_disc, mag_data_vort)
print '--PLOT/SAVE PHASE DIFFERENCES'
# Create projections along z-axis:
projector = SimpleProjector(dim)
# Get analytic solutions:
phase_map_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
phase_map_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
# Create phasemapper:
phasemapper = PhaseMapperFDFC(a, projector.dim_uv, padding=0)
# Create figure:
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Difference of the Fourier space approach from the analytical solution', fontsize=20)
# Create norm for the plots:
bounds = np.array([-100, -50, -25, -5, 0, 5, 25, 50, 100])
norm = BoundaryNorm(bounds, RdBu.N)
# Disc:
phase_map_num_disc = phasemapper(projector(mag_data_disc))
phase_map_num_disc.unit = 'mrad'
phase_diff_disc = phase_map_num_disc - phase_map_ana_disc
RMS_disc = np.sqrt(np.mean(phase_diff_disc.phase**2))
phase_diff_disc.display_phase('Homog. magn. disc, RMS = {:3.2f} mrad'.format(RMS_disc),
axis=axes[0], limit=np.max(bounds), norm=norm)
axes[0].set_aspect('equal')
# Create norm for the plots:
bounds = np.array([-3, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 3])
norm = BoundaryNorm(bounds, RdBu.N)
# Vortex:
phase_map_num_vort = phasemapper(projector(mag_data_vort))
phase_map_num_vort.unit = 'mrad'
phase_diff_vort = phase_map_num_vort - phase_map_ana_vort
phase_diff_vort -= np.mean(phase_diff_vort.phase)
RMS_vort = np.sqrt(np.mean(phase_diff_vort.phase**2))
phase_diff_vort.display_phase(u'Vortex state disc, RMS = {:3.2f} µrad'.format(RMS_vort*1000),
axis=axes[1], limit=np.max(bounds), norm=norm)
axes[1].set_aspect('equal')
plt.show()
plt.figtext(0.15, 0.2, 'a)', fontsize=30)
plt.figtext(0.52, 0.2, 'b)', fontsize=30)
plt.savefig(directory + '/ch5-2-fourier_phase_differe_no_padding.png', bbox_inches='tight')
# Create phasemapper:
phasemapper = PhaseMapperFDFC(a, projector.dim_uv, padding=10)
# Create figure:
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Difference of the Fourier space approach from the analytical solution', fontsize=20)
# Create norm for the plots:
bounds = np.array([-3, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 3])
norm = BoundaryNorm(bounds, RdBu.N)
# Disc:
phase_map_num_disc = phasemapper(projector(mag_data_disc))
phase_map_num_disc.unit = 'mrad'
phase_diff_disc = phase_map_num_disc - phase_map_ana_disc
RMS_disc = np.sqrt(np.mean(phase_diff_disc.phase**2))
phase_diff_disc.display_phase(u'Homog. magn. disc, RMS = {:3.2f} µrad'.format(RMS_disc*1000),
axis=axes[0], limit=np.max(bounds), norm=norm)
axes[0].set_aspect('equal')
# Create norm for the plots:
bounds = np.array([-3, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 3])
norm = BoundaryNorm(bounds, RdBu.N)
# Vortex:
phase_map_num_vort = phasemapper(projector(mag_data_vort))
phase_map_num_vort.unit = 'mrad'
phase_diff_vort = phase_map_num_vort - phase_map_ana_vort
phase_diff_vort -= np.mean(phase_diff_vort.phase)
RMS_vort = np.sqrt(np.mean(phase_diff_vort.phase**2))
phase_diff_vort.display_phase(u'Vortex state disc, RMS = {:3.2f} µrad'.format(RMS_vort*1000),
axis=axes[1], limit=np.max(bounds), norm=norm)
axes[1].set_aspect('equal')
plt.show()
plt.figtext(0.15, 0.2, 'c)', fontsize=30)
plt.figtext(0.52, 0.2, 'd)', fontsize=30)
plt.savefig(directory + '/ch5-2-fourier_phase_differe_padding_10.png', bbox_inches='tight')
###############################################################################################
print 'CH5-2 FOURIER PADDING VARIATION'
# Input parameters:
a = 1.0 # in nm
phi = pi/2
dim = (16, 128, 128) # in px (z, y, x)
center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts with 0!
radius = dim[1]/4 # in px
height = dim[0]/2 # in px
key = 'ch5-2-fourier_padding_mag_dist'
if key in data_shelve and not force_calculation:
print '--LOAD MAGNETIC DISTRIBUTIONS'
(mag_data_disc, mag_data_vort) = data_shelve[key]
else:
print '--CREATE MAGNETIC DISTRIBUTIONS'
# Create magnetic shape (4 times the size):
a_big = a / 2
dim_big = (dim[0]*2, dim[1]*2, dim[2]*2)
center_big = (dim_big[0]/2-0.5, dim_big[1]/2.-0.5, dim_big[2]/2.-0.5)
radius_big = dim_big[1]/4 # in px
height_big = dim_big[0]/2 # in px
mag_shape = mc.Shapes.disc(dim_big, center_big, radius_big, height_big)
# Create MagData (4 times the size):
mag_data_disc = MagData(a_big, mc.create_mag_dist_homog(mag_shape, phi))
mag_data_vort = MagData(a_big, mc.create_mag_dist_vortex(mag_shape, center_big))
# Scale mag_data, grid spacing and dimensions:
mag_data_disc.scale_down()
mag_data_vort.scale_down()
print '--SAVE MAGNETIC DISTRIBUTIONS'
# Shelve magnetic distributions:
data_shelve[key] = (mag_data_disc, mag_data_vort)
# Create projections along z-axis:
projector = SimpleProjector(dim)
# Get analytic solutions:
phase_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
phase_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
# List of applied paddings:
padding_list = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
print '--LOAD/CREATE PADDING SERIES OF HOMOG. MAGN. DISC'
data_disc = np.zeros((3, len(padding_list)))
data_disc[0, :] = padding_list
for (i, padding) in enumerate(padding_list):
key = ', '.join(['Padding->RMS|duration', 'Fourier', 'padding={}'.format(padding_list[i]),
'a={}'.format(a), 'dim={}'.format(dim), 'phi={}'.format(phi), 'disc'])
if key in data_shelve and not force_calculation:
data_disc[:, i] = data_shelve[key]
else:
print '----calculate and save padding =', padding_list[i]
phasemapper = PhaseMapperFDFC(a, projector.dim_uv, padding=padding_list[i])
start_time = time.time()
phase_num_disc = phasemapper(projector(mag_data_disc))
data_disc[2, i] = time.time() - start_time
phase_map_diff = phase_num_disc - phase_ana_disc
phase_map_diff.unit = 'mrad'
phase_map_diff.display_phase()
data_disc[1, i] = np.sqrt(np.mean(phase_map_diff.phase**2))*1E3 # *1E3: rad to mraddd
data_shelve[key] = data_disc[:, i]
print '--PLOT/SAVE PADDING SERIES OF HOMOG. MAGN. DISC'
# Create figure:
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Variation of the padding (homog. magn. disc)', fontsize=20)
# Plot RMS against padding:
axes[0].axhline(y=0.18, linestyle='--', color='g', label='RMS [mrad] (real space)')
axes[0].plot(data_disc[0], data_disc[1], 'go-', label='RMS [mrad] (Fourier space)')
axes[0].set_xlabel('padding', fontsize=15)
axes[0].set_ylabel('RMS [mrad]', fontsize=15)
axes[0].set_xlim(-0.5, 16.5)
axes[0].set_ylim(-5, 45)
axes[0].xaxis.set_major_locator(MaxNLocator(nbins=10, integer=True))
axes[0].tick_params(axis='both', which='major', labelsize=14)
axes[0].legend()
# Plot zoom inset:
ins_axis_d = plt.axes([0.2, 0.35, 0.25, 0.4])
ins_axis_d.axhline(y=0.18, linestyle='--', color='g')
ins_axis_d.plot(data_disc[0], data_disc[1], 'go-')
ins_axis_d.set_yscale('log')
ins_axis_d.set_xlim(5.5, 16.5)
ins_axis_d.set_ylim(0.1, 1.1)
ins_axis_d.tick_params(axis='both', which='major', labelsize=14)
# Plot duration against padding:
axes[1].plot(data_disc[0], data_disc[2], 'bo-')
axes[1].set_xlabel('padding', fontsize=15)
axes[1].set_ylabel('duration [s]', fontsize=15)
axes[1].set_xlim(-0.5, 16.5)
axes[1].set_ylim(-0.05, 2.5)
axes[1].xaxis.set_major_locator(MaxNLocator(nbins=10, integer=True))
axes[1].yaxis.set_major_locator(MaxNLocator(nbins=10))
axes[1].tick_params(axis='both', which='major', labelsize=14)
plt.show()
plt.figtext(0.15, 0.13, 'a)', fontsize=30)
plt.figtext(0.57, 0.17, 'b)', fontsize=30)
plt.savefig(directory + '/ch5-2-disc_padding_variation.png', bbox_inches='tight')
print '--LOAD/CREATE PADDING SERIES OF VORTEX STATE DISC'
data_vort = np.zeros((3, len(padding_list)))
data_vort[0, :] = padding_list
for (i, padding) in enumerate(padding_list):
key = ', '.join(['Padding->RMS|duration', 'Fourier', 'padding={}'.format(padding_list[i]),
'a={}'.format(a), 'dim={}'.format(dim), 'phi={}'.format(phi), 'vort'])
if key in data_shelve and not force_calculation:
data_vort[:, i] = data_shelve[key]
else:
print '----calculate and save padding =', padding_list[i]
phasemapper = PhaseMapperFDFC(a, projector.dim_uv, padding=padding_list[i])
start_time = time.time()
phase_num_vort = phasemapper(projector(mag_data_vort))
data_vort[2, i] = time.time() - start_time
phase_map_diff = phase_num_vort - phase_ana_vort
phase_map_diff -= np.mean(phase_map_diff.phase)
phase_map_diff.unit = 'mrad'
phase_map_diff.display_phase()
data_vort[1, i] = np.sqrt(np.mean(phase_map_diff.phase**2))*1E3 # *1E3: rad to mrad
data_shelve[key] = data_vort[:, i]
print '--PLOT/SAVE PADDING SERIES OF VORTEX STATE DISC'
# Create figure:
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Variation of the padding (Vortex state disc)', fontsize=20)
# Plot RMS against padding:
axes[0].axhline(y=0.22, linestyle='--', color='g', label='RMS [mrad] (real space)')
axes[0].plot(data_vort[0], data_vort[1], 'go-', label='RMS [mrad] (Fourier space)')
axes[0].set_xlabel('padding', fontsize=15)
axes[0].set_ylabel('RMS [mrad]', fontsize=15)
axes[0].set_xlim(-0.5, 16.5)
axes[0].set_ylim(-5, 45)
axes[0].tick_params(axis='both', which='major', labelsize=14)
axes[0].xaxis.set_major_locator(MaxNLocator(nbins=10, integer=True))
axes[0].legend()
# Plot zoom inset:
ins_axis_v = plt.axes([0.2, 0.35, 0.25, 0.4])
ins_axis_v.axhline(y=0.22, linestyle='--', color='g')
ins_axis_v.plot(data_vort[0], data_vort[1], 'go-')
ins_axis_v.set_yscale('log')
ins_axis_v.set_xlim(5.5, 16.5)
ins_axis_v.set_ylim(0.1, 1.1)
ins_axis_v.tick_params(axis='both', which='major', labelsize=14)
# Plot duration against padding:
axes[1].plot(data_vort[0], data_vort[2], 'bo-')
axes[1].set_xlabel('padding', fontsize=15)
axes[1].set_ylabel('duration [s]', fontsize=15)
axes[1].set_xlim(-0.5, 16.5)
axes[1].set_ylim(-0.05, 2.5)
axes[1].xaxis.set_major_locator(MaxNLocator(nbins=10, integer=True))
axes[1].yaxis.set_major_locator(MaxNLocator(nbins=10))
axes[1].tick_params(axis='both', which='major', labelsize=14)
plt.show()
plt.figtext(0.15, 0.13, 'c)', fontsize=30)
plt.figtext(0.57, 0.17, 'd)', fontsize=30)
plt.savefig(directory + '/ch5-2-vortex_padding_variation.png', bbox_inches='tight')
###############################################################################################
print 'CLOSING SHELVE\n'
# Close shelve:
data_shelve.close()
scripts/publications/paper 1/ch5-4-comparison_of_methods_new.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Created on Fri Jul 26 14:37:30 2013 @author: Jan"""
import os
import numpy as np
from numpy import pi
import matplotlib.pyplot as plt
from matplotlib.ticker import NullLocator, LogLocator, LogFormatter
import pyramid
import pyramid.magcreator as mc
import pyramid.analytic as an
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector
from pyramid.phasemapper import PhaseMapperRDFC, PhaseMapperFDFC
from pyramid.kernel import Kernel
import time
import shelve
import logging
import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
force_calculation = False
n = 50
def get_time(pm, proj, mag_data, n):
t = []
for i in range(n):
start = time.clock()
pm(proj(mag_data))
t.append(time.clock()-start)
t, dt = np.mean(t), np.std(t)
phase_map = pm(proj(mag_data))
return phase_map, t, dt
print '\nACCESS SHELVE'
# Create / Open databank:
directory = '../../../output/paper 1'
if not os.path.exists(directory):
os.makedirs(directory)
data_shelve = shelve.open(directory + '/paper_1_shelve')
###############################################################################################
print 'CH5-4 METHOD COMPARISON'
key = 'ch5-4-method_comparison_new'
if key in data_shelve and not force_calculation:
print '--LOAD METHOD DATA'
(data_disc_fourier0, data_vort_fourier0,
data_disc_fourier3, data_vort_fourier3,
data_disc_fourier10, data_vort_fourier10,
data_disc_real, data_vort_real) = data_shelve[key]
else:
# Input parameters:
steps = 4
a = 0.5 # in nm
phi = pi/2
dim = (32, 256, 256) # in px (z, y, x)
center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts at 0!
radius = dim[1]/4 # in px
height = dim[0]/2 # in px
print '--CREATE MAGNETIC SHAPE'
mag_shape = mc.Shapes.disc(dim, center, radius, height)
# Create MagData (4 times the size):
print '--CREATE MAG. DIST. HOMOG. MAGN. DISC'
mag_data_disc = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
print '--CREATE MAG. DIST. VORTEX STATE DISC'
mag_data_vort = MagData(a, mc.create_mag_dist_vortex(mag_shape, center))
# Create Data Arrays:
dim_list = [dim[2]/2**i for i in range(steps)]
data_disc_fourier0 = np.vstack((dim_list, np.zeros((3, steps))))
data_vort_fourier0 = np.vstack((dim_list, np.zeros((3, steps))))
data_disc_fourier3 = np.vstack((dim_list, np.zeros((3, steps))))
data_vort_fourier3 = np.vstack((dim_list, np.zeros((3, steps))))
data_disc_fourier10 = np.vstack((dim_list, np.zeros((3, steps))))
data_vort_fourier10 = np.vstack((dim_list, np.zeros((3, steps))))
data_disc_real = np.vstack((dim_list, np.zeros((3, steps))))
data_vort_real = np.vstack((dim_list, np.zeros((3, steps))))
for i in range(steps):
# Scale mag_data, grid spacing and dimensions:
dim = mag_data_disc.dim
a = mag_data_disc.a
center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px, index starts at 0!
radius = dim[1]/4 # in px
height = dim[0]/2 # in px
print '--a =', a, 'nm', 'dim =', dim
# Create projector along z-axis and phasemapper:
projector = SimpleProjector(dim)
pm_fourier0 = PhaseMapperFDFC(a, projector.dim_uv, padding=0)
pm_fourier3 = PhaseMapperFDFC(a, projector.dim_uv, padding=3)
pm_fourier10 = PhaseMapperFDFC(a, projector.dim_uv, padding=7)
pm_real = PhaseMapperRDFC(Kernel(a, projector.dim_uv))
print '----CALCULATE RMS/DURATION HOMOG. MAGN. DISC'
# Analytic solution:
phase_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
# Fourier unpadded:
phase_num_disc, t, dt = get_time(pm_fourier0, projector, mag_data_disc, n)
data_disc_fourier0[2, i], data_disc_fourier0[3, i] = t, dt
print '------time (disc, fourier0) =', data_disc_fourier0[2, i]
phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
print 'phase mean:', np.mean(phase_num_disc.phase)
data_disc_fourier0[1, i] = np.sqrt(np.mean(phase_diff_disc.phase**2))
# Fourier padding 3:
phase_num_disc, t, dt = get_time(pm_fourier3, projector, mag_data_disc, n)
data_disc_fourier3[2, i], data_disc_fourier3[3, i] = t, dt
print '------time (disc, fourier3) =', data_disc_fourier3[2, i]
phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
print 'phase mean:', np.mean(phase_num_disc.phase)
data_disc_fourier3[1, i] = np.sqrt(np.mean(phase_diff_disc.phase**2))
# Fourier padding 10:
phase_num_disc, t, dt = get_time(pm_fourier10, projector, mag_data_disc, n)
data_disc_fourier10[2, i], data_disc_fourier10[3, i] = t, dt
print '------time (disc, fourier10) =', data_disc_fourier10[2, i]
phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
print 'phase mean:', np.mean(phase_num_disc.phase)
data_disc_fourier10[1, i] = np.sqrt(np.mean(phase_diff_disc.phase**2))
# Real space disc:
phase_num_disc, t, dt = get_time(pm_real, projector, mag_data_disc, n)
data_disc_real[2, i], data_disc_real[3, i] = t, dt
print '------time (disc, real space) =', data_disc_real[2, i]
phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
print 'phase mean:', np.mean(phase_num_disc.phase)
data_disc_real[1, i] = np.sqrt(np.mean(phase_diff_disc.phase**2))
print 'TIME:', data_disc_real[2, i]
print 'RMS%:', np.sqrt(np.mean(((phase_ana_disc-phase_num_disc).phase /
phase_ana_disc.phase)**2))*100, '%'
print '----CALCULATE RMS/DURATION VORTEX STATE DISC'
# Analytic solution:
phase_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
# Fourier unpadded:
phase_num_vort, t, dt = get_time(pm_fourier0, projector, mag_data_vort, n)
phase_fft0 = phase_num_vort
data_vort_fourier0[2, i], data_vort_fourier0[3, i] = t, dt
print '------time (vortex, fourier0) =', data_vort_fourier0[2, i]
phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
print 'phase mean:', np.mean(phase_num_vort.phase)
phase_diff_vort -= np.mean(phase_diff_vort.phase)
data_vort_fourier0[1, i] = np.sqrt(np.mean(phase_diff_vort.phase**2))
# Fourier padding 3:
phase_num_vort, t, dt = get_time(pm_fourier3, projector, mag_data_vort, n)
phase_fft3 = phase_num_vort
data_vort_fourier3[2, i], data_vort_fourier3[3, i] = t, dt
print '------time (vortex, fourier3) =', data_vort_fourier3[2, i]
phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
print 'phase mean:', np.mean(phase_num_vort.phase)
phase_diff_vort -= np.mean(phase_diff_vort.phase)
data_vort_fourier3[1, i] = np.sqrt(np.mean(phase_diff_vort.phase**2))
# Fourier padding 10:
phase_num_vort, t, dt = get_time(pm_fourier10, projector, mag_data_vort, n)
phase_fft10 = phase_num_vort
data_vort_fourier10[2, i], data_vort_fourier10[3, i] = t, dt
print '------time (vortex, fourier10) =', data_vort_fourier10[2, i]
phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
print 'phase mean:', np.mean(phase_num_vort.phase)
phase_diff_vort -= np.mean(phase_diff_vort.phase)
data_vort_fourier10[1, i] = np.sqrt(np.mean(phase_diff_vort.phase**2))
# Real space disc:
phase_num_vort, t, dt = get_time(pm_real, projector, mag_data_vort, n)
phase_real = phase_num_vort
data_vort_real[2, i], data_vort_real[3, i] = t, dt
print '------time (vortex, real space) =', data_vort_real[2, i]
phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
print 'phase mean:', np.mean(phase_num_vort.phase)
# phase_diff_vort -= np.mean(phase_diff_vort.phase)
data_vort_real[1, i] = np.sqrt(np.mean(phase_diff_vort.phase**2))
print 'TIME:', data_disc_real[2, i]
print 'RMS%:', np.sqrt(np.mean(((phase_ana_vort-phase_num_vort).phase /
phase_ana_vort.phase)**2))*100, '%'
# Scale down for next iteration:
mag_data_disc.scale_down()
mag_data_vort.scale_down()
print '--SHELVE METHOD DATA'
data_shelve[key] = (data_disc_fourier0, data_vort_fourier0,
data_disc_fourier3, data_vort_fourier3,
data_disc_fourier10, data_vort_fourier10,
data_disc_real, data_vort_real)
print '--PLOT/SAVE METHOD DATA'
# Plot using shared rows and colums:
fig, axes = plt.subplots(2, 2, sharex='col', sharey='row', figsize=(12, 8))
fig.tight_layout(rect=(0.05, 0.05, 0.95, 0.95))
fig.suptitle('Method Comparison', fontsize=20)
# Plot RMS against a (disc) [top/left]:
axes[1, 0].set_xscale('log')
axes[1, 0].set_yscale('log')
axes[1, 0].plot(data_disc_fourier0[0], data_disc_fourier0[1], ':bs')
axes[1, 0].plot(data_disc_fourier3[0], data_disc_fourier3[1], ':bo')
axes[1, 0].plot(data_disc_fourier10[0], data_disc_fourier10[1], ':b^')
axes[1, 0].plot(data_disc_real[0], data_disc_real[1], '--ro')
axes[1, 0].set_xlabel('grid size [px]', fontsize=15)
axes[1, 0].set_ylabel('RMS [mrad]', fontsize=15)
axes[1, 0].set_xlim(25, 350)
axes[1, 0].tick_params(axis='both', which='major', labelsize=14)
axes[1, 0].xaxis.set_major_locator(LogLocator(base=2))
axes[1, 0].xaxis.set_major_formatter(LogFormatter(base=2))
axes[1, 0].xaxis.set_minor_locator(NullLocator())
axes[1, 0].grid()
# Plot duration against a (disc) [bottom/left]:
plt.tick_params(axis='both', which='major', labelsize=14)
axes[0, 0].set_xscale('log')
axes[0, 0].set_yscale('log')
axes[0, 0].errorbar(data_disc_fourier0[0], data_disc_fourier0[2],
yerr=0*data_vort_fourier3[3], fmt=':bs')
axes[0, 0].errorbar(data_disc_fourier3[0], data_disc_fourier3[2],
yerr=0*data_vort_fourier3[3], fmt=':bo')
axes[0, 0].errorbar(data_disc_fourier10[0], data_disc_fourier10[2],
yerr=0*data_vort_fourier3[3], fmt=':b^')
axes[0, 0].errorbar(data_disc_real[0], data_disc_real[2],
yerr=0*data_vort_fourier3[3], fmt='--ro')
axes[0, 0].set_title('Homog. magn. disc', fontsize=18)
axes[0, 0].set_ylabel('duration [s]', fontsize=15)
axes[0, 0].set_xlim(25, 350)
axes[0, 0].set_ylim(1E-4, 1E1)
axes[0, 0].tick_params(axis='both', which='major', labelsize=14)
axes[0, 0].xaxis.set_major_locator(LogLocator(base=2))
axes[0, 0].xaxis.set_major_formatter(LogFormatter(base=2))
axes[0, 0].xaxis.set_minor_locator(NullLocator())
axes[0, 0].grid()
# Plot RMS against a (vortex) [top/right]:
axes[1, 1].set_xscale('log')
axes[1, 1].set_yscale('log')
axes[1, 1].plot(data_vort_fourier0[0], data_vort_fourier0[1], ':bs',
label='Fourier padding=0')
axes[1, 1].plot(data_vort_fourier3[0], data_vort_fourier3[1], ':bo',
label='Fourier padding=3')
axes[1, 1].plot(data_vort_fourier10[0], data_vort_fourier10[1], ':b^',
label='Fourier padding=10')
axes[1, 1].plot(data_vort_real[0], data_vort_real[1], '--ro',
label='Real space')
axes[1, 1].set_xlabel('grid size [px]', fontsize=15)
axes[1, 1].set_xlim(25, 350)
axes[1, 1].tick_params(axis='both', which='major', labelsize=14)
axes[1, 1].xaxis.set_major_locator(LogLocator(base=2))
axes[1, 1].xaxis.set_major_formatter(LogFormatter(base=2))
axes[1, 1].xaxis.set_minor_locator(NullLocator())
axes[1, 1].grid()
# Plot duration against a (vortex) [bottom/right]:
axes[0, 1].set_xscale('log')
axes[0, 1].set_yscale('log')
axes[0, 1].errorbar(data_vort_fourier0[0], data_vort_fourier0[2], yerr=0*data_vort_fourier0[3],
fmt=':bs', label='Fourier padding=0')
axes[0, 1].errorbar(data_vort_fourier3[0], data_vort_fourier3[2], yerr=0*data_vort_fourier3[3],
fmt=':bo', label='Fourier padding=3')
axes[0, 1].errorbar(data_vort_fourier10[0], data_vort_fourier10[2], yerr=0*data_vort_fourier3[3],
fmt=':b^', label='Fourier padding=10')
axes[0, 1].errorbar(data_vort_real[0], data_vort_real[2], yerr=0*data_vort_fourier3[3],
fmt='--ro', label='Real space')
axes[0, 1].set_title('Vortex state disc', fontsize=18)
axes[0, 1].set_xlim(25, 350)
axes[0, 1].set_ylim(1E-4, 1E1)
axes[0, 1].tick_params(axis='both', which='major', labelsize=14)
axes[0, 1].xaxis.set_major_locator(LogLocator(base=2))
axes[0, 1].xaxis.set_major_formatter(LogFormatter(base=2))
axes[0, 1].xaxis.set_minor_locator(NullLocator())
axes[0, 1].grid()
# Add legend:
axes[1, 1].legend(bbox_to_anchor=(0, 0, 0.955, 0.615), bbox_transform=fig.transFigure,
prop={'size': 12})
# Save figure as .png:
plt.show()
plt.figtext(0.12, 0.85, 'a)', fontsize=30)
plt.figtext(0.57, 0.85, 'b)', fontsize=30)
plt.figtext(0.12, 0.15, 'c)', fontsize=30)
plt.figtext(0.57, 0.15, 'd)', fontsize=30)
plt.savefig(directory + '/ch5-3-method comparison.png', bbox_inches='tight')
###############################################################################################
print 'CLOSING SHELVE\n'
# Close shelve:
data_shelve.close()
scripts/publications/poster_pof.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""
Created on Wed Feb 05 19:19:19 2014
@author: Jan
"""
import os
from numpy import pi
import numpy as np
import pyramid
import pyramid.magcreator as mc
import pyramid.analytic as an
from pyramid.magdata import MagData
from pyramid.phasemapper import pm, PhaseMapperRDFC
from pyramid.projector import SimpleProjector
from pyramid.kernel import Kernel
from pyramid.phasemap import PhaseMap
from time import clock
import logging
import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
a = 10.0 # nm
dim = (1, 9, 9)
mag_shape = mc.Shapes.pixel(dim, (int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)))
mag_data_x = MagData(a, mc.create_mag_dist_homog(mag_shape, 0))
mag_data_y = MagData(a, mc.create_mag_dist_homog(mag_shape, pi/2))
pm(mag_data_x).display_phase()
pm(mag_data_y).display_phase()
a = 1
dim = (128, 128, 128)
center = (dim[0]/2, dim[1]/2, dim[2]/2)
radius = dim[1]/4 # in px
mag_shape = mc.Shapes.sphere(dim, center, radius)
mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape, pi/4))
projector = SimpleProjector(dim)
phasemapper = PhaseMapperRDFC(Kernel(a, projector.dim_uv))
start = clock()
phase_map = phasemapper(projector(mag_data))
print 'TIME:', clock() - start
phase_ana = an.phase_mag_sphere(dim, a, pi/4, center, radius)
phase_diff = (phase_ana - phase_map).phase
print 'RMS%:', np.sqrt(np.mean((phase_diff)**2))/phase_ana.phase.max()*100
print 'max:', phase_ana.phase.max()
print 'RMS:', np.sqrt(np.mean((phase_diff)**2))
PhaseMap(a, phase_diff).display_phase()
mag_data.scale_down(2)
mag_data.quiver_plot()
mag_data.quiver_plot3d()
phase_map.display_phase()
scripts/reconstruction/reconstruct_random_pixels.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Create random magnetic distributions."""
import os
import random as rnd
import numpy as np
from numpy import pi
import pyramid
import pyramid.magcreator as mc
from pyramid.magdata import MagData
from pyramid.phasemapper import pm
import pyramid.reconstruction as rc
import logging
import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
# Input parameters:
n_pixel = 5
dim = (1, 16, 16)
b_0 = 1 # in T
a = 10.0 # in nm
rnd.seed(18)
# Create empty MagData object and add random pixels:
mag_data = MagData(a, np.zeros((3,)+dim))
for i in range(n_pixel):
pixel = (rnd.randrange(dim[0]), rnd.randrange(dim[1]), rnd.randrange(dim[2]))
mag_shape = mc.Shapes.pixel(dim, pixel)
phi = 2 * pi * rnd.random()
magnitude = rnd.random()
mag_data += MagData(a, mc.create_mag_dist_homog(mag_shape, phi, magnitude))
# Plot magnetic distribution, phase map and holographic contour map:
mag_data.quiver_plot()
phase_map = pm(mag_data)
phase_map.display_combined('Generated Distribution', gain=10)
# Reconstruct the magnetic distribution:
mag_data_rec = rc.optimize_simple_leastsq(phase_map, mag_data.get_mask(), b_0, lam=1E-4, order=1)
# Display the reconstructed phase map and holography image:
phase_map_rec = pm(mag_data_rec)
phase_map_rec.display_combined('Reconstructed Distribution', gain=10)
scripts/reconstruction/reconstruction_sparse_cg_singular_values.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Created on Thu Apr 24 11:00:00 2014 @author: Jan"""
import os
import numpy as np
from numpy import pi
import pyramid
from pyramid.magdata import MagData
from pyramid.projector import YTiltProjector, XTiltProjector
from pyramid.dataset import DataSet
from pyramid.kernel import Kernel
from pyramid.phasemap import PhaseMap
from pyramid.forwardmodel import ForwardModel
import logging
import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
###################################################################################################
print('--Singular value decomposition')
a = 1.
b_0 = 1.
dim = (3, 3, 3)
dim_uv = dim[1:3]
count = 8
tilts_full = np.linspace(-pi/2, pi/2, num=count/2, endpoint=False)
tilts_miss = np.linspace(-pi/3, pi/3, num=count/2, endpoint=False)
phase_zero = PhaseMap(a, np.zeros(dim_uv))
projectors_xy_full, projectors_x_full, projectors_y_full = [], [], []
projectors_xy_miss, projectors_x_miss, projectors_y_miss = [], [], []
projectors_x_full.extend([XTiltProjector(dim, tilt) for tilt in tilts_full])
projectors_y_full.extend([YTiltProjector(dim, tilt) for tilt in tilts_full])
projectors_xy_full.extend(projectors_x_full)
projectors_xy_full.extend(projectors_y_full)
projectors_x_miss.extend([XTiltProjector(dim, tilt) for tilt in tilts_miss])
projectors_y_miss.extend([YTiltProjector(dim, tilt) for tilt in tilts_miss])
projectors_xy_miss.extend(projectors_x_miss)
projectors_xy_miss.extend(projectors_y_miss)
###################################################################################################
projectors = projectors_xy_full
###################################################################################################
size_2d = np.prod(dim_uv)
size_3d = np.prod(dim)
data = DataSet(a, dim, b_0)
[data.append(phase_zero, projectors[i]) for i in range(len(projectors))]
y = data.phase_vec
F = ForwardModel(data)
M = np.asmatrix([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
#MTM = M.T * M + lam * np.asmatrix(np.eye(3*size_3d))
U, s, V = np.linalg.svd(M) # TODO: M or MTM?
for i in range(len(s)):
print 'Singular value:', s[i]
title = 'Singular value = {:g}'.format(s[i])
MagData(data.a, np.array(V[i, :]).reshape((3,)+dim)).quiver_plot3d(title)
scripts/reconstruction/reconstruction_sparse_cg_test.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Created on Fri Feb 28 14:25:59 2014 @author: Jan"""
import os
import numpy as np
from numpy import pi
from time import clock
import pyramid
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
from pyramid.projector import YTiltProjector, XTiltProjector
from pyramid.dataset import DataSet
from pyramid.regularisator import ZeroOrderRegularisator
import pyramid.magcreator as mc
import pyramid.reconstruction as rc
import logging
import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
###################################################################################################
print('--Generating input phase_maps')
a = 10.
b_0 = 1.
dim = (8, 32, 32)
dim_uv = dim[1:3]
count = 16
center = (int(dim[0]/2)-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
mag_shape = mc.Shapes.disc(dim, center, dim[2]/4, dim[2]/2)
magnitude = mc.create_mag_dist_vortex(mag_shape)
mag_data = MagData(a, magnitude)
#vortex_shape = mc.Shapes.disc(dim, (3.5, 9.5, 9.5), 5, 4)
#sphere_shape = mc.Shapes.sphere(dim, (3.5, 22.5, 9.5), 3)
#slab_shape = mc.Shapes.slab(dim, (3.5, 15.5, 22.5), (5, 8, 8))
#mag_data = MagData(a, mc.create_mag_dist_vortex(vortex_shape, (3.5, 9.5, 9.5)))
#mag_data += MagData(a, mc.create_mag_dist_homog(sphere_shape, pi/4, pi/4))
#mag_data += MagData(a, mc.create_mag_dist_homog(slab_shape, -pi/6))
mag_data.quiver_plot3d()
tilts_full = np.linspace(-pi/2, pi/2, num=count/2, endpoint=False)
tilts_miss = np.linspace(-pi/3, pi/3, num=count/2, endpoint=False)
projectors_xy_full, projectors_x_full, projectors_y_full = [], [], []
projectors_xy_miss, projectors_x_miss, projectors_y_miss = [], [], []
projectors_x_full.extend([XTiltProjector(dim, tilt) for tilt in tilts_full])
projectors_y_full.extend([YTiltProjector(dim, tilt) for tilt in tilts_full])
projectors_xy_full.extend(projectors_x_full)
projectors_xy_full.extend(projectors_y_full)
projectors_x_miss.extend([XTiltProjector(dim, tilt) for tilt in tilts_miss])
projectors_y_miss.extend([YTiltProjector(dim, tilt) for tilt in tilts_miss])
projectors_xy_miss.extend(projectors_x_miss)
projectors_xy_miss.extend(projectors_y_miss)
###################################################################################################
projectors = projectors_xy_full
noise = 0.0
###################################################################################################
print('--Setting up data collection')
data = DataSet(a, dim, b_0)
data.projectors = projectors
data.phase_maps = data.create_phase_maps(mag_data)
if noise != 0:
for phase_map in data.phase_maps:
phase_map += PhaseMap(a, np.random.normal(0, noise, dim_uv))
###################################################################################################
print('--Test simple solver')
lam = 1E-4
reg = ZeroOrderRegularisator(lam)
start = clock()
mag_data_opt = rc.optimize_sparse_cg(data, regularisator=reg, maxiter=100, verbosity=2)
print 'Time:', str(clock()-start)
mag_data_opt.quiver_plot3d()
(mag_data_opt - mag_data).quiver_plot3d()
###################################################################################################
print('--Plot stuff')
phase_maps_opt = data.create_phase_maps(mag_data_opt)
#data.display_phase()
#data.display_phase(phase_maps_opt)
phase_diffs = [(data.phase_maps[i]-phase_maps_opt[i]) for i in range(len(data.phase_maps))]
[phase_diff.display_phase() for phase_diff in phase_diffs]
scripts/test methods/compare_methods.py deleted 100644 → 0
View file @ 04d9ed45
#! python
# -*- coding: utf-8 -*-
"""Compare the different methods to create phase maps."""
import os
import time
import numpy as np
from numpy import pi
import pyramid
import pyramid.magcreator as mc
import pyramid.analytic as an
from pyramid.projector import SimpleProjector
from pyramid.phasemapper import PhaseMapperRDRC, PhaseMapperRDFC, PhaseMapperFDFC
from pyramid.kernel import Kernel
from pyramid.magdata import MagData
import logging
import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
# Input parameters:
b_0 = 1.0 # in T
a = 1.0 # in nm
phi = pi/4
numcore = False
padding = 10
gain = 10
dim = (128, 128, 128) # in px (z, y, x)
# Create magnetic shape:
geometry = 'sphere'
if geometry == 'slab':
center = (dim[0]/2.-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z,y,x) index starts at 0!
width = (dim[0]/2, dim[1]/2, dim[2]/2) # in px (z, y, x)
mag_shape = mc.Shapes.slab(dim, center, width)
phase_map_ana = an.phase_mag_slab(dim, a, phi, center, width, b_0)
elif geometry == 'disc':
center = (dim[0]/2.-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z,y,x) index starts at 0!
radius = dim[1]/4 # in px
height = dim[0]/2 # in px
mag_shape = mc.Shapes.disc(dim, center, radius, height)
phase_map_ana = an.phase_mag_disc(dim, a, phi, center, radius, height, b_0)
elif geometry == 'sphere':
center = (dim[0]/2.-0.5, dim[1]/2.-0.5, dim[2]/2.-0.50) # in px (z, y, x) index starts with 0!
radius = dim[1]/4 # in px
mag_shape = mc.Shapes.sphere(dim, center, radius)
phase_map_ana = an.phase_mag_sphere(dim, a, phi, center, radius, b_0)
# Create MagData object and projector:
mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
projector = SimpleProjector(dim)
projection = projector(mag_data)
# Construct PhaseMapper objects:
pm_real = PhaseMapperRDRC(Kernel(a, projector.dim_uv, b_0), numcore=numcore)
pm_conv = PhaseMapperRDFC(Kernel(a, projector.dim_uv, b_0))
pm_four = PhaseMapperFDFC(a, projector.dim_uv, b_0, padding=padding)
# Get times for different approaches:
start_time = time.time()
phase_map_real = pm_real(projection)
print 'Time for RDRC: ', time.time() - start_time
start_time = time.time()
phase_map_conv = pm_conv(projection)
print 'Time for RDFC: ', time.time() - start_time
start_time = time.time()
phase_map_four = pm_four(projection)
print 'Time for FDFC: ', time.time() - start_time
# Display the combinated plots with phasemap and holography image:
phase_map_ana.display_combined('Analytic Solution', gain=gain)
phase_map_real.display_combined('RDRC', gain=gain)
phase_map_conv.display_combined('RDFC', gain=gain)
phase_map_four.display_combined('FDFC', gain=gain)
# Plot differences to the analytic solution:
phase_map_diff_real = phase_map_real - phase_map_ana
phase_map_diff_conv = phase_map_conv - phase_map_ana
phase_map_diff_four = phase_map_four - phase_map_ana
RMS_real = np.std(phase_map_diff_real.phase)
RMS_conv = np.std(phase_map_diff_conv.phase)
RMS_four = np.std(phase_map_diff_four.phase)
phase_map_diff_real.display_phase('RDRC difference (RMS = {:3.2e})'.format(RMS_real))
phase_map_diff_conv.display_phase('RDFC difference (RMS = {:3.2e})'.format(RMS_conv))
phase_map_diff_four.display_phase('FDFC difference (RMS = {:3.2e})'.format(RMS_four))
scripts/test methods/compare_pixel_fields.py deleted 100644 → 0
View file @ 04d9ed45
#! python
# -*- coding: utf-8 -*-
"""Compare the phase map of one pixel for different approaches."""
import os
import numpy as np
import pyramid
import pyramid.magcreator as mc
from pyramid.projector import SimpleProjector
from pyramid.phasemapper import PhaseMapperRDFC, PhaseMapperFDFC
from pyramid.kernel import Kernel
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
import matplotlib.pyplot as plt
from matplotlib.ticker import IndexLocator
import logging
import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
directory = '../../output/magnetic distributions'
if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
a = 1.0 # in nm
phi = 0 # in rad
dim = (1, 8, 8)
pixel = (0, int(dim[1]/2), int(dim[2]/2))
limit = 0.35
def get_fourier_kernel():
PHI_0 = 2067.83 # magnetic flux in T*nm²
b_0 = 1
coeff = - 1j * b_0 * a**2 / (2*PHI_0)
nyq = 1 / a # nyquist frequency
f_u = np.linspace(0, nyq/2, dim[2]/2.+1)
f_v = np.linspace(-nyq/2, nyq/2, dim[1], endpoint=False)
f_uu, f_vv = np.meshgrid(f_u, f_v)
phase_fft = coeff * f_vv / (f_uu**2 + f_vv**2 + 1e-30)
# Transform to real space and revert padding:
phase_fft = np.fft.ifftshift(phase_fft, axes=0)
phase_fft_kernel = np.fft.fftshift(np.fft.irfft2(phase_fft), axes=(0, 1))
return phase_fft_kernel
# Create magnetic data and projector:
mag_data = MagData(a, mc.create_mag_dist_homog(mc.Shapes.pixel(dim, pixel), phi))
mag_data.save_to_llg(directory + '/mag_dist_single_pixel.txt')
projector = SimpleProjector(dim)
# Kernel of a disc in real space:
phase_map_disc = PhaseMapperRDFC(Kernel(a, projector.dim_uv, geometry='disc'))(mag_data)
phase_map_disc.unit = 'mrad'
phase_map_disc.display_phase('Phase of one Pixel (Disc)', limit=limit)
# Kernel of a slab in real space:
phase_map_slab = PhaseMapperRDFC(Kernel(a, projector.dim_uv, geometry='slab'))(mag_data)
phase_map_slab.unit = 'mrad'
phase_map_slab.display_phase('Phase of one Pixel (Slab)', limit=limit)
# Kernel of the Fourier method:
phase_map_fft = PhaseMapperFDFC(a, projector.dim_uv, padding=0)(mag_data)
phase_map_fft.unit = 'mrad'
phase_map_fft.display_phase('Phase of one Pixel (Fourier)', limit=limit)
# Kernel of the Fourier method, calculated directly:
phase_map_fft_kernel = PhaseMap(a, get_fourier_kernel(), 'mrad')
phase_map_fft_kernel.display_phase('Phase of one Pixel (Fourier Kernel)', limit=limit)
# Kernel differences:
print 'Fourier Kernel, direct and indirect method are identical:', \
np.all(phase_map_fft_kernel.phase - phase_map_fft.phase) == 0
(phase_map_disc-phase_map_fft).display_phase('Phase difference of one Pixel (Disc - Fourier)')
# Cross section plots of real space kernels:
fig = plt.figure()
axis = fig.add_subplot(1, 1, 1)
x = np.linspace(-dim[1]/a/2, dim[1]/a/2-1, dim[1])
y_ft = phase_map_fft.phase[:, dim[1]/2]
y_re = phase_map_disc.phase[:, dim[1]/2]
axis.axhline(0, color='k')
axis.plot(x, y_re, label='Real space method')
axis.plot(x, y_ft, label='Fourier method')
axis.grid()
axis.legend()
axis.set_title('Real Space Kernel')
axis.set_xlim(-dim[1]/2, dim[1]/2-1)
axis.xaxis.set_major_locator(IndexLocator(base=dim[1]/8, offset=0))
# Cross section plots of Fourier space kernels:
fig = plt.figure()
axis = fig.add_subplot(1, 1, 1)
x = range(dim[1])
k_re = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))[:, 0]**2
k_ft = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0))[:, 0]**2
axis.axhline(0, color='k')
axis.plot(x, k_re, label='Real space method')
axis.plot(x, k_ft, label='Fourier method')
axis.grid()
axis.legend()
axis.set_title('Fourier Space Kernel')
axis.set_xlim(0, dim[1]-1)
axis.xaxis.set_major_locator(IndexLocator(base=dim[1]/8, offset=0))
scripts/test methods/projector_test.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Created on Fri Aug 15 08:09:10 2014 @author: Jan"""
import numpy as np
import pyramid.magcreator as mc
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector, XTiltProjector, YTiltProjector
from pyramid.phasemapper import PhaseMapperRDFC
from pyramid.kernel import Kernel
dim = (16, 24, 32)
center = (int(dim[0]/2), int(dim[1]/2), int(dim[2]/2))
width = (2, 8, 16)
dim_uv = (32, 48)
a = 2
b_0 = 1.5
shape = mc.Shapes.slab(dim, center, width)
mag_data = MagData(a, mc.create_mag_dist_homog(shape, phi=np.pi/6, theta=np.pi/3))
phase_mapper = PhaseMapperRDFC(Kernel(a, dim_uv, b_0))
# PROJECTOR TESTING
phase_mapper(SimpleProjector(dim, 'z', dim_uv)(mag_data)).display_phase('z simple')
phase_mapper(XTiltProjector(dim, 0, dim_uv)(mag_data)).display_phase('z (x-tilt)')
phase_mapper(SimpleProjector(dim, 'x', dim_uv)(mag_data)).display_phase('x simple')
phase_mapper(YTiltProjector(dim, np.pi/2, dim_uv)(mag_data)).display_phase('x (y-tilt)')
phase_mapper(XTiltProjector(dim, np.pi/2, dim_uv)(mag_data)).display_phase('y (x-tilt)')
phase_mapper(SimpleProjector(dim, 'y', dim_uv)(mag_data)).display_phase('y simple')
scripts/test methods/test_arb_projection.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Created on Tue Sep 03 12:55:40 2013 @author: Jan"""
import os
import numpy as np
from numpy import pi
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator, NullLocator
import pyramid
import itertools
import logging
import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
dim = (8, 8)
offset = (dim[0]/2., dim[1]/2.)
phi = 1.000000001*pi/4
field = np.zeros(dim)
def get_position(p, m, b, size):
x, y = np.array(p)[:, 0]+0.5, np.array(p)[:, 1]+0.5
return (y-m*x-b)/np.sqrt(m**2+1) + size/2.
def get_impact(pos, r, size):
return [x for x in np.arange(np.floor(pos-r), np.floor(pos+r)+1, dtype=int) if 0 <= x < size]
def get_weight(d, r):
return (1 - d/r) * (d < r)
direction = (-np.cos(phi), np.sin(phi))
xi = range(dim[0])
yj = range(dim[1])
ii, jj = np.meshgrid(xi, yj)
r = 1/np.sqrt(np.pi)
m = direction[0]/direction[1]
b = offset[0] - m * offset[1]
voxels = list(itertools.product(yj, xi))
positions = get_position(voxels, m, b, 2*dim[0])
weights = []
for i, voxel in enumerate(voxels):
voxel_weights = []
impacts = get_impact(positions[i], r, dim[0])
for impact in impacts:
distance = np.abs(impact+0.5 - positions[i])
voxel_weights.append((impact, get_weight(distance, r)))
weights.append((voxel, voxel_weights))
pixels = np.floor(positions).astype(int)
pixel_hits = zip(set(pixels), [list(pixels).count(i) for i in set(pixels)])
def Y(x):
return direction[0]/direction[1] * (x - offset[1]) + offset[0]
def Y_perp(x):
return - direction[1]/direction[0] * (x - offset[1]) + offset[0]
def X(y):
return direction[1]/direction[0] * (y - offset[0]) + offset[1]
x = np.linspace(-2, 2, 10)
fig = plt.figure()
axis = fig.add_subplot(1, 1, 1, aspect='equal')
axis.pcolormesh(field, cmap='PuRd')
axis.grid(which='both', color='k', linestyle='-')
x = np.linspace(0, dim[1])
y = Y(x)
y_perp = Y_perp(x)
axis.plot(x, y, '-r', linewidth=2)
axis.plot(x, y_perp, '-g', linewidth=2)
axis.set_xlim(0, dim[1])
axis.set_ylim(0, dim[0])
axis.xaxis.set_major_locator(MaxNLocator(nbins=dim[1], integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=dim[0], integer=True))
for i, p in enumerate(voxels):
if 0 <= positions[i] < dim[0]:
color = 'k'
else:
color = 'r'
plt.annotate('{:1.1f}'.format(positions[i]), p, size=8, color=color)
plt.show()
fig = plt.figure()
axis = fig.add_subplot(1, 1, 1, aspect='equal')
axis.scatter(positions, 0.5*np.ones_like(positions))
axis.grid(which='both', color='k', linestyle='-')
axis.vlines((0.5*dim[0], 1.5*dim[0]), 0, 1, colors='r', linewidth=3)
axis.set_xlim(0, 2*dim[0])
axis.set_ylim(0, 1)
axis.xaxis.set_major_locator(MaxNLocator(nbins=2*dim[0], integer=True))
axis.yaxis.set_major_locator(NullLocator())
for i, px in enumerate(pixel_hits):
plt.annotate('{:d}'.format(px[1]), (px[0], 0.6), size=8)
plt.show()
setup.cfg 0 → 100644
View file @ 6e3128b8
# -*- coding: utf-8 -*-
# setup.cfg
# :copyright: Copyright 2020 Jan Caron
# TODO: pip install .[io] is not working because hyperspy depends on trait which has no wheels on PyPI at the moment...
# TODO: See https://github.com/hyperspy/hyperspy/issues/2315 and https://github.com/enthought/traits/issues/357
# TODO: Add hyperspy back as soon (if?) this is resolved... Until then: install hyperspy with conda!
# TODO: mayavi has the same exact problem and tvtk is apparently also a problem...
# CONFIGURATION FOR TESTING:
[aliases]
test = pytest
[coverage:run]
branch = True
source = src/empyre
omit =
tests/*
[flake8]
max-line-length = 120
ignore =
# module import not at top of file
E402
# closing bracket does not match visual indentation
E124
# continuation line with same indent as next logical line
E125
# missing whitespace around arithmetic operator
E226
# line break before binary operator
W503
# line break after binary operator
W504
# do not use variables named ‘l’, ‘O’, or ‘I’
E741
per-file-ignores =
*/__init__.py: F401, F403, F405, F821
# F401: module imported but unused
# F403: 'from module import *' used; unable to detect undefined names
# F405: Name may be undefined, or defined from star imports: module
# F821: undefined name 'name'
setup.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Setup for testing, building, distributing and installing the 'Pyramid'-package"""
import numpy
import os
import sys
import sysconfig
import subprocess
from distutils.command.build import build
from Cython.Distutils import build_ext
from setuptools import setup
from setuptools import find_packages
from setuptools.extension import Extension
from pyramid._version import __version__
class custom_build(build):
'''Custom build command'''
def make_hgrevision(self, target):
output = subprocess.Popen(["hg", "id", "-i"], stdout=subprocess.PIPE).communicate()[0]
hgrevision_cc = file(str(target), "w")
hgrevision_cc.write('hg_revision = "{0}"'.format(output.strip()))
hgrevision_cc.close()
def run(self):
build.run(self)
print 'creating hg_revision.txt'
self.make_hgrevision(os.path.join('build', get_build_path('lib'), 'hg_revision.txt'))
def get_build_path(dname):
'''Returns the name of a distutils build directory'''
path = "{dirname}.{platform}-{version[0]}.{version[1]}"
return path.format(dirname=dname, platform=sysconfig.get_platform(), version=sys.version_info)
def get_files(rootdir):
'''Returns a list of .py-files inside rootdir'''
filepaths = []
for root, dirs, files in os.walk(rootdir):
for filename in files:
if filename.endswith('.py'):
filepaths.append(os.path.join(root, filename))
return filepaths
print '\n-------------------------------------------------------------------------------'
setup(
name = 'Pyramid',
version = __version__,
description = 'PYthon based Reconstruction Algorithm for MagnetIc Distributions',
author = 'Jan Caron',
author_email = 'j.caron@fz-juelich.de',
url = 'fz-juelich.de',
packages = find_packages(exclude=['tests']),
include_dirs = [numpy.get_include()],
requires = ['numpy', 'matplotlib', 'mayavi'],
scripts = get_files('scripts'),
test_suite = 'tests',
cmdclass = {'build_ext': build_ext, 'build': custom_build},
ext_package = 'pyramid/numcore',
ext_modules = [
Extension('phasemapper_core', ['pyramid/numcore/phasemapper_core.pyx'],
include_dirs = [numpy.get_include(), numpy.get_numarray_include()],
extra_compile_args=['-march=native', '-mtune=native']
)
]
)
print '-------------------------------------------------------------------------------\n'
src/empyre/__init__.py 0 → 100644
View file @ 6e3128b8
# -*- coding: utf-8 -*-
# Copyright 2020 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Subpackage containing functionality for visualisation of multidimensional fields."""
from . import fields
from . import io
from . import models
from . import reconstruct
from . import vis
from . import utils
__version__ = "0.3.4"
__git_revision__ = "undefined"
import logging
_log = logging.getLogger(__name__)
_log.info(f'Imported EMPyRe V-{__version__}')
del logging
__all__ = ['fields', 'io', 'models', 'reconstruct', 'vis', 'utils']
src/empyre/fields/__init__.py 0 → 100644
View file @ 6e3128b8
# -*- coding: utf-8 -*-
# Copyright 2020 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Subpackage containing container classes for multidimensional fields and ways to create them."""
from .field import *
from .shapes import *
from .vectors import *
__all__ = []
__all__.extend(field.__all__)
__all__.extend(shapes.__all__)
__all__.extend(vectors.__all__)
del field
del shapes
del vectors
src/empyre/fields/field.py 0 → 100644
View file @ 6e3128b8
# -*- coding: utf-8 -*-
# Copyright 2020 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module provides a container class for multidimensional scalar or vector fields."""
import logging
from numbers import Number
from numpy.lib.mixins import NDArrayOperatorsMixin
import numpy as np
try:
# numpy >= 2
from numpy.lib.array_utils import normalize_axis_tuple
except ModuleNotFoundError:
# numpy < 2
from numpy.core.numeric import normalize_axis_tuple
from scipy import ndimage
from ..utils import Quaternion
__all__ = ['Field']
class Field(NDArrayOperatorsMixin):
"""Container class for storing multidimensional scalar or vector fields.
The `Field` class is a sophisticated wrapper around a multidimensional numpy array. The user can access the
underlying data via the `data` attribute, or by using the `numpy.asarray` function.
`Field` defines the `ufunc` interface for numpys vast library of universal functions. This includes all operator
definitions, i.e., you can add, subtract, multiply and devide `Field` objects element-wise. Note that for vector
fields, `Field` will handle all vector components separately and for ufuncs with two inputs, the other object will
be broadcast accordingly. It is therefore also possible to add (or multiply) vector and scalar fields, assuming
their dimensions `dim` match (their `shape` can't match, because it includes the number of components `ncomp` for
the vector field). For ufuncs that reduce the output shape (e.g. `numpy.sum`), if the `axis` parameter is set to
`None` (default), component axes for vector fields are kept (they can however explicitely reduced by including
them in the `axis` paremeter, e.g. by setting `axis = -1`).
`Field` objects can be indexed like numpy arrays and scalar fields (`vector=False`) behave as expected. A vector
field (`vector=True`) can be indexed by only the non-component indices (see `dim`) and the according components of
the specified point are returned as a tuple. Slicing with integers will return a `Field` object with reduced
dimensionality (or a scalar/tuple if all spatial dimensions are reduced) that also drops the associated `scales`
as needed. New axes (indexing with `None` or `numpy.newaxis`) is disabled because that would mean that an unknown
scale would need to be added (you can always create a new `Field` object with an appropriately prepared ndarray).
Some callables exposed by numpy's public API (e.g. `numpy.mean(...)`) will treat `Field` objects as if they
were identical to their underlying numpy array (by envoking the `__array__` function which just returns `data`).
The result often is a simple `numpy.ndarray` and not a `Field` object. Some functions (e.g. `clip`) are also defined
on `Field` as methods. They might behave differently and will return `Field` objects instead, so take care and
refer to their docstrings for further information!
# TODO: Too verbose? In documentation instead?
Attributes
----------
data: np.ndarray
The underlying data array of the field.
scale: tuple of float
Scaling along the dimensions of the underlying data.
vector: bool
True if the field is a vector field, False if it is a scalar field.
shape: tuple of ints
Shape of the underlying data. Includes the number of components as the first dimension if the field is a
vector field.
dim: tuple of ints
Dimensions of the underlying data. Only includes spatial dimensions, NOT the number of vector components!
(Use `ncomp` to get the number of components, if any, `shape` for the full shape of the underlying array).
ncomp: None or tuple of ints
Number of components if the field is a vector field (can be 2 or 3), None if it is a scalar field. The component
axis is always the last axis (index `-1`) of the underlying data array!
Notes
-----
See https://docs.scipy.org/doc/numpy/reference/arrays.classes.html for information about ndarray subclassing!
See https://docs.scipy.org/doc/numpy-1.13.0/neps/ufunc-overrides.html for information about __array_ufunc__!
See https://numpy.org/neps/nep-0018-array-function-protocol.html for information about __array_function__!
"""
_log = logging.getLogger(__name__ + '.Field')
@property
def data(self):
return self.__data
@property
def vector(self):
return self.__vector
@vector.setter
def vector(self, vector):
assert isinstance(vector, bool), 'vector has to be a boolean!'
self.__vector = vector
@property
def scale(self):
return self.__scale
@scale.setter
def scale(self, scale):
if isinstance(scale, Number): # Scale is the same for each dimension!
self.__scale = (float(scale),) * len(self.dim)
elif isinstance(scale, tuple):
ndim = len(self.dim)
assert len(scale) == ndim, f'Each of the {ndim} dimensions needs a scale, but {scale} was given!'
self.__scale = scale
else:
raise AssertionError(f'Scaling has to be a number or a tuple of numbers, was {scale} instead!')
@property
def shape(self):
return self.data.shape
@property
def dim(self):
if self.vector:
return self.shape[:-1]
else:
return self.shape
@property
def ncomp(self):
return self.shape[-1] if self.vector else 0
@property
def comp(self): # TODO: Function or property? Philosophical question?
# Return empty list for scalars (ncomp is 0) and a list of components as scalar fields for a vector field:
return [Field(self.data[..., i], self.scale, vector=False) for i in range(self.ncomp)]
@property
def amp(self): # TODO: Function or property? Philosophical question?
if self.vector:
amp = np.sqrt(np.sum([comp.data**2 for comp in self.comp], axis=0))
else:
amp = np.abs(self.data)
return Field(amp, self.scale, vector=False)
@property
def mask(self): # TODO: Function or property? Philosophical question?
return Field(np.where(np.asarray(self.amp) > 0, True, False), self.scale, vector=False)
def __init__(self, data, scale=1.0, vector=False):
self._log.debug('Calling __init__')
self.__data = data
self.vector = vector # Set vector before scale, because scale needs to know if vector (via calling dim)!
self.scale = scale
if self.vector:
assert self.ncomp in (2, 3), 'Only 2- or 3-component vector fields are supported!'
self._log.debug('Created ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
data_string = str(self.data) # String of JUST the data array, without metadata (compared to repr)!
return f'{self.__class__.__name__}(data={data_string}, scale={self.scale}, vector={self.vector})'
def __str__(self):
self._log.debug('Calling __repr__')
ncomp_string = f', ncomp={self.ncomp}' if self.vector else ''
return f'{self.__class__.__name__}(dim={self.dim}, scale={self.scale}, vector={self.vector}{ncomp_string})'
def __array__(self):
self._log.debug('Calling __array__')
return self.data
def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
self._log.debug('Calling __array_ufunc__')
# outputs = kwargs.pop('out', ()) # TODO: Necessary?
outputs = kwargs.pop('out', (None,)*ufunc.nout) # Defaults to tuple of None (currently: nout=1 all the time)
outputs_arr = tuple([np.asarray(out) if isinstance(out, Field) else out for out in outputs])
if all(a is None for a in outputs_arr):
outputs_arr = None
# Cannot handle items that have __array_ufunc__ (other than our own).
for item in inputs + outputs:
if hasattr(item, '__array_ufunc__') and not isinstance(item, Field): # Something else with __array_ufunc__:
if type(item).__array_ufunc__ is not np.ndarray.__array_ufunc__: # Can't handle other overrides:
return NotImplemented
# TODO: BIGGEST NOTE HERE: Delegate work to ndarray.__array_ufunc__!
# TODO: For this to work, we have to make sure that we dispatch input as ndarrays, NOT Fields!
# Convert all Field inputs to simple ndarrays to avoid infinite recursion!
# ndarray.__array_ufunc__ delegates to the called ufunc method ONLY if all inputs and outputs are ndarrays
# (or have no __array_ufunc__ method, e.g. a scalar), so that's what we need to ensure here.:
# inputs = tuple([inp.view(np.ndarray) if isinstance(inp, Field) else inp for inp in inputs])
# TODO: possibly need to sort out which input (if two) is a Field (vector or scalar?), order matters?
# TODO: DOCUMENT: most things are same as in ndarrays, but if only one input is a vector field, try broadcast!
# TODO: for security, newaxis is not allowed (most other indexing works though), because scale would be unknown!
# 1 input (has to be a Field, otherwise we wouldn't be here):
self._log.debug(f'__array_ufunc__ inputs: {len(inputs)}')
self._log.info(f'ufunc: {ufunc}, method: {method}')
self._log.info(f'inputs: {inputs}')
self._log.info(f'outputs: {outputs}')
self._log.info(f'kwargs: {kwargs}')
if len(inputs) == 1:
field = inputs[0]
scale_new = field.scale
vector_new = field.vector
# Preprocess axis keyword if it exists:
axis = kwargs.get('axis', False) # Default must not be None, because None is a possible setting!
full_reduction = False
# for ufunc.outer and ufunc.at:
if axis is not False:
ax_full = tuple(range(len(field.dim))) # All axes (minus a possible component axis for vector Fields)!
ax_full_wc = tuple(range(len(field.dim) + 1)) # All axes WITH component axis (does not need to exist)!
axis = ax_full if axis is None else axis # This keeps possible components untouched if axis was None!
# axis=-1 reduces over the vector components, if they exist
# Takes care of pot. neg. indices, ensures tuple!
axis = normalize_axis_tuple(axis, len(field.shape))
kwargs['axis'] = axis # Put normalized axis back into kwargs!
if tuple(sorted(axis)) in (ax_full, ax_full_wc):
full_reduction = True # Full reduction (or reduction to just components) takes place:
if ax_full_wc[-1] in axis: # User explicitely wants component reduction (can only be true for vector):
vector_new = False # Force scalar field!
scale_new = tuple([s for i, s in enumerate(field.scale) if i not in axis]) # Drop axis from scale!
inputs_arr = np.asarray(field) # Convert inputs that are Fields to ndarrays to avoid recursion!
data_new = getattr(ufunc, method)(inputs_arr, out=outputs_arr, **kwargs)
if full_reduction: # Premature return because the result is no longer a Field:
return data_new
# More than 1 input (at least one has to be a Field, otherwise we wouldn't be here):
elif len(inputs) > 1:
is_field = [isinstance(inp, Field) for inp in inputs]
is_vector = [getattr(inp, 'vector', False) for inp in inputs]
# Determine scale:
if np.sum(is_field) > 1: # More than one input is a Field objects:
scales = [inp.scale for i, inp in enumerate(inputs) if is_field[i]] # Only takes scales of Field obj.!
scale_new = scales[0]
err_msg = f'Scales of all Field objects must match! Given scales: {scales}!'
if not all(scale == scale_new for scale in scales):
raise ValueError(err_msg)
else: # Only one input is a field, pick the scale of that one:
scale_new = inputs[np.argmax(is_field)].scale # argmax returns the index of first True!
# Determine vector:
vector_new = True if np.any(is_vector) else False # Output is vector field if any input is a vector field!
if np.sum(is_vector) > 1: # More than one input is a vector Field objects:
ncomps = [inp.ncomp for i, inp in enumerate(inputs) if is_vector[i]] # Only takes ncomp of v.-Fields!
err_msg = f'# of components of all Field objects must match! Given ncomps: {ncomps}!'
if not all(ncomp == ncomps[0] for ncomp in ncomps):
raise ValueError(err_msg)
# Append new axis at the end of non vector objects to broadcast to components:
if np.any(is_vector):
inputs = list(inputs)
for i, inp in enumerate(inputs):
if not is_vector[i] and not isinstance(inp, Number): # Numbers work for broadcasting anyway:
if len(np.asarray(inp).shape) == len(scale_new): # No. of dimensions w/o comp., have to match!
inputs[i] = np.asarray(inputs[i])[..., np.newaxis] # Broadcasting, try to cast as ndarray!
inputs = tuple(inputs)
# Convert inputs that are Fields to ndarrays to avoid recursion and determine data_new:
inputs_arr = tuple([np.asarray(inp) if isinstance(inp, Field) else inp for inp in inputs])
data_new = getattr(ufunc, method)(*inputs_arr, out=outputs_arr, **kwargs)
# Return results:
result = Field(data_new, scale_new, vector_new)
return result
def __getitem__(self, index):
self._log.debug('Calling __getitem__')
# Pre-processing index:
index = index if isinstance(index, tuple) else (index,) # Make sure item is a tuple!
if None in index:
raise Exception('Field does not support indexing with newaxis/None! Please cast as ndarray first!')
if len(index) < len(self.shape) and Ellipsis not in index: # Ellipsis is IMPLICIT at the end if index < dim:
index += (Ellipsis,)
if Ellipsis in index: # Expand Ellipsis (...) into slice(None) (:) to make iterating consistent:
index = list(index)
i = index.index(Ellipsis)
missing_dims = len(self.shape) - len(index) # Number of non-specified dimensions
index[i:i+1] = [slice(None)] * (missing_dims + 1) # +1 because at least Ellipsis is replaced!
index = tuple(index)
# Indexing with a scalar drops the dimension, indexing with slice always keeps the dimension:
index_scale = index[:-1] if self.vector else index # Disregard last index for vector fields (has no scale)!
scale_new = ()
for i, item in enumerate(index_scale):
if isinstance(item, slice): # Indexing with slice keeps the dimension and therefore the scale:
scale_new += (self.scale[i],)
# Get data with classic indexing from underlying data array:
data_new = self.data[index]
# Determine vector state of output:
vector_new = self.vector
if self.vector and isinstance(index[-1], Number): # For vector fields if last index (component) is dropped:
vector_new = False
# Return:
if not scale_new: # scale_new=(), i.e. full reduction of all dimensions (not regarding possible vector comp.):
if isinstance(data_new, Number): # full reduction:
return data_new
else: # Only important for vector fields:
return tuple(data_new) # return single vector as tuple
else: # Return Field object
return Field(data_new, scale_new, vector=vector_new)
@classmethod
def from_scalar_fields(cls, scalar_list):
"""Create a vector `Field` object from a list of scalar `Fields`.
Parameters
----------
scalar_list : list
List of scalar `Field` objects (`vector=False`).
Returns
-------
vector_field: `Field`
A vector `Field` object containing the input scalar fields as components.
"""
cls._log.debug('Calling from_scalar_fields')
assert len(scalar_list) in (2, 3), 'Needs two or three scalar fields as components for vector field!'
assert all(isinstance(item, Field) for item in scalar_list), 'All items have to be Field objects!'
assert all(not item.vector for item in scalar_list), 'All items have to be scalar fields!'
assert all(item.scale == scalar_list[0].scale for item in scalar_list), 'Scales of fields must match!'
return cls(np.stack(scalar_list, axis=-1), scalar_list[0].scale, vector=True)
@classmethod
def from_signal(cls, signal, scale=None, vector=None, comp_pos=-1):
"""Convert a :class:`~hyperspy.signals.Signal` object to a :class:`~.Field` object.
Parameters
----------
signal: :class:`~hyperspy.signals.Signal`
The :class:`~hyperspy.signals.Signal` object which should be converted to Field.
Returns
-------
field: :class:`~.Field`
A :class:`~.Field` object containing the loaded data.
scale: tuple of float, optional
Scaling along the dimensions of the underlying data. If not provided, will be read from the axes_manager.
vector: boolean, optional
If set to True, forces the signal to be interpreted as a vector field. EMPyRe will check if the first axis
is named 'vector components' (EMPyRe saves vector fields like this). If this is the case, vector will be
automatically set to True and the signal will also be interpreted as a vector field.
comp_pos: int, optoinal
The index of the axis containing the vector components (if `vector=True`). EMPyRe needs this to be the last
axis (index `-1`, which is the default). In case another position is given, the vector component will be
moved to the last axis. Old Pyramid files will have this axis at index `0`, so use this for backwards
compatibilty.
Notes
-----
Signals recquire the hyperspy package!
"""
cls._log.debug('Calling from_signal')
data = signal.data
if vector and comp_pos != -1: # component axis should be last, but is currently first -> roll to the end:
data = np.moveaxis(data, source=comp_pos, destination=-1)
if vector is None: # Automatic detection:
vector = True if signal.axes_manager[0].name == 'vector components' else False
if scale is None: # If not provided, try to read from axes_manager, one less axis if vector!:
scale = tuple([signal.axes_manager[i].scale for i in range(len(data.shape) - vector)])
return cls(data, scale, vector)
def to_signal(self):
"""Convert :class:`~.Field` data into a HyperSpy signal.
Returns
-------
signal: :class:`~hyperspy.signals.Signal`
Representation of the :class:`~.Field` object as a HyperSpy Signal.
Notes
-----
This method recquires the hyperspy package!
"""
self._log.debug('Calling to_signal')
try: # Try importing HyperSpy:
import hyperspy.api as hs
except ImportError:
self._log.error('This method recquires the hyperspy package!')
return
# Create signal:
signal = hs.signals.BaseSignal(self.data) # All axes are signal axes!
# Set axes:
if self.vector:
signal.axes_manager[0].name = 'vector components'
for i in range(len(self.dim)):
ax = i+1 if self.vector else i # take component axis into account if vector!
num = ['x', 'y', 'z'][i] if len(self.dim) <= 3 else i
signal.axes_manager[ax].name = f'axis {num}'
signal.axes_manager[ax].scale = self.scale[i]
signal.axes_manager[ax].units = 'nm'
signal.metadata.Signal.title = f"EMPyRe {'vector' if self.vector else 'scalar'} Field"
return signal
def copy(self):
"""Returns a copy of the :class:`~.Field` object.
Returns
-------
field: :class:`~.Field`
A copy of the :class:`~.Field`.
"""
self._log.debug('Calling copy')
return Field(self.data.copy(), self.scale, self.vector)
def rotate(self, angle, axis='z', **kwargs):
"""Rotate the :class:`~.Field`, based on :meth:`~scipy.ndimage.rotate`.
Rotation direction is from the first towards the second axis. Works for 2D and 3D Fields.
Parameters
----------
angle : float
The rotation angle in degrees.
axis: {'x', 'y', 'z'}, optional
The axis around which the vector field is rotated. Default is 'z'. Ignored for 2D Fields.
Returns
-------
field: :class:`~.Field`
The rotated :class:`~.Field`.
Notes
-----
All additional kwargs are passed through to :meth:`~scipy.ndimage.rotate`.
The `reshape` parameter, controlling if the output shape is adapted so that the input array is contained
completely in the output is False per default, contrary to :meth:`~scipy.ndimage.rotate`,
where it is True.
"""
self._log.debug('Calling rotate')
assert len(self.dim) in (2, 3), 'rotate is currently only defined for 2D and 3D Fields!'
kwargs.setdefault('reshape', False) # Default here is no reshaping!
if len(self.dim) == 2: # For 2D, there are only 2 axes:
axis = 'z' # Overwrite potential argument if 2D!
axes = (0, 1) # y and x
else: # 3D case:
# order of axes is important for permutation, to get positive levi_civita
axes = {'x': (0, 1), 'y': (2, 0), 'z': (1, 2)}[axis]
sc_0, sc_1 = self.scale[axes[0]], self.scale[axes[1]]
assert sc_0 == sc_1, f'rotate needs the scales in the rotation plane to be equal (they are {sc_0} & {sc_1})!'
np_angle = angle
if axis in ('x', 'z'):
np_angle *= -1
if not self.vector: # Scalar field:
data_new = ndimage.rotate(self.data, np_angle, axes=axes, **kwargs)
else: # Vector field:
# Rotate coordinate system:
comps = [np.asarray(comp) for comp in self.comp]
if self.ncomp == 3:
data_new = np.stack([ndimage.rotate(c, np_angle, axes=axes, **kwargs) for c in comps], axis=-1)
# Up till now, only the coordinates are rotated, now we need to rotate the vectors inside the voxels:
rot_axis = {'x': 2, 'y': 1, 'z': 0}[axis]
i, j, k = axes[0], axes[1], rot_axis # next line only works if i != j != k
levi_civita = int((j-i) * (k-i) * (k-j) / (np.abs(j-i) * np.abs(k-i) * np.abs(k-j)))
# Create Quaternion, note that they have (x, y, z) order instead of (z, y, x):
quat_axis = tuple([levi_civita if i == rot_axis else 0 for i in (2, 1, 0)])
quat = Quaternion.from_axisangle(quat_axis, np.deg2rad(angle))
# T needed b.c. ordering!
data_new = quat.matrix.dot(data_new.reshape((-1, 3)).T).T.reshape(data_new.shape)
elif self.ncomp == 2:
u_comp, v_comp = comps
u_rot = ndimage.rotate(u_comp, np_angle, axes=axes, **kwargs)
v_rot = ndimage.rotate(v_comp, np_angle, axes=axes, **kwargs)
# Up till now, only the coordinates are rotated, now we need to rotate the vectors inside the voxels:
ang_rad = np.deg2rad(angle)
u_mix = np.cos(ang_rad)*u_rot - np.sin(ang_rad)*v_rot
v_mix = np.sin(ang_rad)*u_rot + np.cos(ang_rad)*v_rot
data_new = np.stack((u_mix, v_mix), axis=-1)
else:
raise ValueError('rotate currently only works for 2- or 3-component vector fields!')
return Field(data_new, self.scale, self.vector)
def rot90(self, k=1, axis='z'):
"""Rotate the :class:`~.Field` 90° around the specified axis (right hand rotation).
Parameters
----------
k : integer
Number of times the array is rotated by 90 degrees.
axis: {'x', 'y', 'z'}, optional
The axis around which the vector field is rotated. Default is 'z'. Ignored for 2D Fields.
Returns
-------
field: :class:`~.Field`
The rotated :class:`~.Field`.
"""
self._log.debug('Calling rot90')
assert axis in ('x', 'y', 'z'), "Wrong input! 'x', 'y', 'z' allowed!"
assert len(self.dim) in (2, 3), 'rotate is currently only defined for 2D and 3D Fields!'
if len(self.dim) == 2: # For 2D, there are only 2 axes:
axis = 'z' # Overwrite potential argument if 2D!
axes = (0, 1) # y and x
else: # 3D case:
axes = {'x': (0, 1), 'y': (0, 2), 'z': (1, 2)}[axis]
sc_0, sc_1 = self.scale[axes[0]], self.scale[axes[1]]
assert sc_0 == sc_1, f'rot90 needs the scales in the rotation plane to be equal (they are {sc_0} & {sc_1})!'
# TODO: Later on, rotation could also flip the scale (not implemented here, yet)!
if axis != 'y':
k = -k # rotation is inverted if around x or z due to y flip compared to numpy
if not self.vector: # Scalar Field:
data_new = np.rot90(self.data, k=k, axes=axes).copy()
else: # Vector Field:
if len(self.dim) == 3: # 3D:
assert self.ncomp == 3, 'rot90 currently only works for vector fields with 3 components in 3D!'
comp_x, comp_y, comp_z = self.comp
# reference:
# https://en.wikipedia.org/wiki/Rotation_matrix#In_three_dimensions
if axis == 'x': # RotMatrix for 90°: [[1, 0, 0], [0, 0, -1], [0, 1, 0]]
comp_x_rot = np.rot90(comp_x, k=k, axes=axes)
comp_y_rot = -np.rot90(comp_z, k=k, axes=axes)
comp_z_rot = np.rot90(comp_y, k=k, axes=axes)
elif axis == 'y': # RotMatrix for 90°: [[0, 0, 1], [0, 1, 0], [-1, 0, 0]]
comp_x_rot = np.rot90(comp_z, k=k, axes=axes)
comp_y_rot = np.rot90(comp_y, k=k, axes=axes)
comp_z_rot = -np.rot90(comp_x, k=k, axes=axes)
elif axis == 'z': # RotMatrix for 90°: [[0, -1, 0], [1, 0, 0], [0, 0, 1]]
comp_x_rot = -np.rot90(comp_y, k=k, axes=axes)
comp_y_rot = np.rot90(comp_x, k=k, axes=axes)
comp_z_rot = np.rot90(comp_z, k=k, axes=axes)
data_new = np.stack((comp_x_rot, comp_y_rot, comp_z_rot), axis=-1)
if len(self.dim) == 2: # 2D:
assert self.ncomp == 2, 'rot90 currently only works for vector fields with 2 components in 2D!'
comp_x, comp_y = self.comp
comp_x_rot = -np.rot90(comp_y, k=k, axes=axes)
comp_y_rot = np.rot90(comp_x, k=k, axes=axes)
data_new = np.stack((comp_x_rot, comp_y_rot), axis=-1)
# Return result:
return Field(data_new, self.scale, self.vector)
def get_vector(self, mask=None):
"""Returns the field as a vector, specified by a mask.
Parameters
----------
mask : :class:`~numpy.ndarray` (boolean)
Masks the pixels from which the entries should be taken. Must be a numpy array for indexing to work!
Returns
-------
vector : :class:`~numpy.ndarray` (N=1)
The vector containing the field of the specified pixels. If the field is a vector field, components are
first vectorised, then concatenated!
"""
self._log.debug('Calling get_vector')
if mask is None: # If not given, take full data:
mask = np.full(self.dim, True)
if self.vector: # Vector field:
return np.ravel([comp.data[mask] for comp in self.comp])
else: # Scalar field:
return np.ravel(self.data[mask])
def set_vector(self, vector, mask=None):
"""Set the field of the masked pixels to the values specified by `vector`.
Parameters
----------
mask : :class:`~numpy.ndarray` (boolean), optional
Masks the pixels from which the field should be taken.
vector : :class:`~numpy.ndarray` (N=1)
The vector containing the field of the specified pixels.
Returns
-------
None
"""
self._log.debug('Calling set_vector')
if mask is None: # If not given, set full data:
mask = np.full(self.dim, True)
if self.vector: # Vector field:
assert np.size(vector) % self.ncomp == 0, 'Vector has to contain all components for every pixel!'
count = np.size(vector) // self.ncomp
for i in range(self.ncomp):
sl = slice(i*count, (i+1)*count)
self.data[..., i][mask] = vector[sl]
else: # Scalar field:
self.data[mask] = vector
def squeeze(self):
"""Squeeze the `Field` object to remove single-dimensional entries in the shape.
Returns
-------
field: `Field`
Squeezed `Field` object. Note that scales associated with squeezed dimensions are also dropped.
Notes
-----
Also works via `numpy.squeeze(field)`, because `numpy.squeeze` searches for a local implementation first and
then uses `_wrapit` to envoke this function here!
"""
self._log.debug('Calling squeeze')
squeezed_indices = np.flatnonzero(np.asarray(self.dim) == 1)
squeezed_data = np.squeeze(self.data)
if squeezed_indices:
self._log.info(f'The following indices were squeezed: {squeezed_indices}')
squeezed_scale = tuple([self.scale[i] for i in range(len(self.dim)) if i not in squeezed_indices])
return Field(squeezed_data, squeezed_scale, self.vector)
def pad(self, pad_width, mode='constant', **kwargs):
"""Pad the `Field` object and increase the size of the underlying array.
Parameters
----------
pad_width : {sequence, array_like, int}
Number of values padded to the edges of each axis. Can be a single number for all, one number per axis or
a tuple `(before, after)` for each axis for finer control.
mode : str, optional
Padding mode (see `numpy.pad` for all options), by default 'constant', which pads with zeros.
Returns
-------
field: `Field`
The padded `Field` object.
"""
if isinstance(pad_width, Number): # Paddding is the same for each dimension (make sure it is a tuple)!
pad_width = (pad_width,) * len(self.dim)
pad_width = [(p, p) if isinstance(p, Number) else p for p in pad_width]
if self.vector: # Append zeros to padding, so component axis stays as is:
pad_width = pad_width + [(0, 0)]
data_new = np.pad(self.data, pad_width, mode, **kwargs)
return Field(data_new, self.scale, self.vector)
def bin(self, n):
"""Bins data of the `Field` to decrease the size of the underlying array by averaging over a number of values.
Parameters
----------
n : {sequence, array_like, int}
Number of entries along each axis over which the average is taken. Can be a single integer for all axes or a
list, specifying the number of entries over which to average for each individual axis.
Returns
-------
field: `Field`
The binned `Field` object.
Notes
-----
Padding takes place before binning to ensure the dimensions are multiples of `n`. The padding mode is `edge`.
"""
self._log.debug('Calling bin')
assert isinstance(n, (int, tuple)), 'n must be a positive integer or a tuple of positive integers!'
if isinstance(n, Number): # Binning is the same for each dimension (make sure it is a tuple)!
n = (n,) * len(self.dim)
assert all([n_i >= 0 for n_i in n]), 'All entries of n must be positive integers!'
# Pad if necessary (use padded 'field' from here on), formula for clarity: (n - dim % n) % n
pad_width = [(0, (n[i] - self.dim[i] % n[i]) % n[i]) for i in range(len(self.dim))]
field = self.pad(pad_width, mode='edge')
# Create new shape used for binning (mean over every second axis will be taken):
bin_shape = list(np.ravel([(field.dim[i]//n[i], n[i]) for i in range(len(field.dim))]))
mean_axes = np.arange(1, 2*len(field.dim), 2) # every 2nd axis!
if self.vector: # Vector field:
bin_shape += [field.ncomp] # Append component axis (they stay unchanged)
# Bin data and scale accordingly:
data_new = field.data.reshape(bin_shape).mean(axis=tuple(mean_axes))
scale_new = tuple([field.scale[i] * n[i] for i in range(len(field.dim))])
return Field(data_new, scale_new, self.vector)
def zoom(self, zoom, **kwargs):
"""Wrapper for the `scipy.ndimage.zoom` function.
Parameters
----------
zoom : float or sequence
Zoom factor along the axes. If a float, `zoom` is the same for each axis. If a sequence, `zoom` needs to
contain one value for each axis.
Returns
-------
field: `Field`
The zoomed in `Field` object.
Notes
-----
As in `scipy.ndimage.zoom`, a spline order can be specified, which defaults to 3.
"""
self._log.debug('Calling zoom')
if not self.vector: # Scalar field:
data_new = ndimage.zoom(self.data, zoom, **kwargs)
else: # Vector field:
comps = [np.asarray(comp) for comp in self.comp]
data_new = np.stack([ndimage.zoom(comp, zoom, **kwargs) for comp in comps], axis=-1)
if isinstance(zoom, Number): # Zoom is the same for each dimension!
zoom = (zoom,) * len(self.dim)
scale_new = tuple([self.scale[i]/z for i, z in enumerate(zoom)])
return Field(data_new, scale_new, self.vector)
def flip(self, axis=None, **kwargs):
"""Returns a `Field` with entries flipped along specified axes.
Parameters
----------
axis : tuple or None, optional
Axes whose entries will be flipped, by default None.
"""
self._log.debug('Calling flip')
if self.vector and axis is None:
axis = tuple(range(len(self.dim)))
axis = normalize_axis_tuple(axis, len(self.shape)) # Make sure, axis is a tuple!
data_new = np.flip(self.data, axis, **kwargs).copy() # Only flips space, not component direction!
if self.vector:
flip_vec = [
-1 if i in axis else 1
for i in reversed(range(self.ncomp))
]
data_new *= flip_vec
return Field(data_new, self.scale, self.vector)
def gradient(self):
"""Returns the gradient of an N-dimensional scalar `Field`. Wrapper around the according numpy function.
Returns
-------
gradients: ndarray or list of ndarray
A set of ndarrays (or a single ndarray for 1D input), corresponding to the derivatives of the field with
respect to each dimension.
Notes
-----
The field is implicitely squeezed before the gradient is calculated!
"""
self._log.debug('Calling gradient')
assert not self.vector, 'Gradient can only be computed from scalar fields!'
squeezed_field = self.squeeze() # Gradient along dimension of length 1 does not work!
gradients = np.gradient(np.asarray(squeezed_field), *squeezed_field.scale)
if len(squeezed_field.dim) == 1: # Only one gradient!
return Field(gradients, squeezed_field.scale, vector=False)
else: # Multidimensional gradient (flip component order with [::-1], so that e.g x/y/z instead of z/y/x):
return Field(np.stack(gradients[::-1], axis=-1), squeezed_field.scale, vector=True)
def curl(self):
"""Returns the curl of an N-dimensional `Field`.
Returns
-------
field: `Field`
The curl/rotation.
Notes
-----
The calculation depends on the input:
3 dimensions, 3 components: Calculates the full 3D rotational vector field!
2 dimensions, 2 components: Calculates the out-of-plane component of the curl as a 2D scalar field!
2 dimensions, scalar: Calculates the planar rotation as a 2D vector field!
"""
self._log.debug('Calling curl')
squeezed_field = self.squeeze()
if len(squeezed_field.dim) == 3: # 3 dimensions:
if squeezed_field.ncomp == 3: # 3 component vector field (standard case):
self._log.debug('input: 3 dimensions, 3 components!')
field_x, field_y, field_z = squeezed_field.comp
gradx_x, grady_x, gradz_x = field_x.gradient().comp
gradx_y, grady_y, gradz_y = field_y.gradient().comp
gradx_z, grady_z, gradz_z = field_z.gradient().comp
curl_x = grady_z - gradz_y
curl_y = gradz_x - gradx_z
curl_z = gradx_y - grady_x
return Field.from_scalar_fields([curl_x, curl_y, curl_z])
else:
raise AssertionError('Can only handle 3 component vector fields in 3D!')
elif len(squeezed_field.dim) == 2: # 2 dimensions (tricky, usually not hardly defined):
if squeezed_field.ncomp == 2: # 2 component vector field (return perpendicular component as scalar field):
self._log.debug('input: 2 dimensions, 2 components!')
field_x, field_y = squeezed_field.comp
gradx_x, grady_x = field_x.gradient().comp
gradx_y, grady_y = field_y.gradient().comp
return gradx_y - grady_x
elif not squeezed_field.vector: # scalar field (return planar components as 2D vector field):
self._log.debug('input: 3 dimensions, scalar field!')
gradx, grady = squeezed_field.gradient().comp
curl_x = grady
curl_y = -gradx
return Field.from_scalar_fields([curl_x, curl_y])
else:
raise AssertionError('Can only handle 3 component vector or scalar fields in 2D!')
else:
raise AssertionError('Can only handle 3D or 2D cases (see documentation for specifics)!')
def clip(self, vmin=None, vmax=None, sigma=None, mask=None):
"""Clip (limit) the values in an array. For vector fields, this will take the amplitude into account.
Parameters
----------
vmin : float, optional
Mimimum value, by default None. Is ignored for vector fields. Will overwrite values found via `sigma` or
`mask`.
vmax : float, optional
Maximum value, by default None. Will overwrite values found via `sigma` or `mask`.
sigma : float, optional
Defines a range in standard deviations, that results in a boolean mask that excludes outliers, by default
None. E.g. `sigma=2` will mark points within the 5% highest amplitude as outliers. `vmin` and `vmax` will be
searched in the remaining region. Is logically combined with `mask` if both are set.
mask : ndarray, optional
Boolean array that directly describes where to search for `vmin` and `vmax`, by default None. Is logically
combined with the mask from `sigma` if both are set.
Returns
-------
field: `Field`
The clipped `Field`.
Notes
-----
In contrast to the corresponding numpy function, `vmin` and `vmax` can both be `None` due to the alternative
clipping strategies employed here.
"""
self._log.debug('Calling clip')
# Get a scalar indicator array for where clipping needs to happen:
if self.vector: # For vector fields, it is the amplitude of the data:
indicator = self.amp.data
else: # For scalar fields this is the data itself:
indicator = self.data
if mask is None: # If no mask is set yet, default to True everywhere:
mask = np.full(self.dim, True)
if sigma is not None: # Mark outliers that are outside `sigma` standard deviations:
sigma_mask = (indicator - indicator.mean()) < (sigma * indicator.std())
mask = np.logical_and(mask, sigma_mask)
# Determine vmin and vmax if they are not set by the user:
indicator_masked = np.where(mask, indicator, np.nan)
if vmin is None:
vmin = np.nanmin(indicator_masked)
if vmax is None:
vmax = np.nanmax(indicator_masked)
if self.vector: # Vector fields need to scale components according to masked amplitude:
# Only vmax is important for vectors! mask_vec broadcast to components!
mask_vec = (indicator <= vmax)[..., None]
data_new = np.where(mask_vec, self.data, vmax * self.data/indicator[..., None]) # Scale outliers to vmax!
else: # For scalar fields, just delegate to the numpy function:
data_new = np.clip(self.data, vmin, vmax)
return Field(data_new, self.scale, self.vector)
src/empyre/fields/shapes.py 0 → 100644
View file @ 6e3128b8
# -*- coding: utf-8 -*-
# Copyright 2020 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
import logging
import numpy as np
from .field import Field
__all__ = ['create_shape_slab', 'create_shape_disc', 'create_shape_ellipse', 'create_shape_ellipsoid',
'create_shape_sphere', 'create_shape_filament', 'create_shape_voxel']
_log = logging.getLogger(__name__)
# TODO: TEST!!!
def create_shape_slab(dim, center=None, width=None, scale=1.0):
"""Creates a Field object with the shape of a slab as a scalar field in arbitrary dimensions.
Attributes
----------
dim : tuple
The dimensions of the grid.
center : tuple, optional
The center of the slab in pixel coordinates.
width : tuple, optional
The width of the slab in pixel coordinates.
scale: tuple of float
Scaling along the dimensions of the underlying data.
"""
_log.debug('Calling __init__')
if center is None:
center = tuple([d/2 for d in dim])
if width is None:
width = tuple([d/2 for d in dim])
assert len(dim) == len(center) == len(width), 'Parameters dim, center and width must have the same dimensions!'
data = np.zeros(dim)
bounds = ()
for i in range(len(dim)):
start = int(np.floor(center[i] - width[i]/2))
stop = int(np.ceil(center[i] + width[i]/2))
bounds += (slice(start, stop),)
data[bounds] = 1
return Field(data=data, scale=scale, vector=False)
def create_shape_disc(dim, center=None, radius=None, height=None, axis=0, scale=1.0):
"""Creates a Field object with the shape of a cylindrical disc in 2D or 3D.
Attributes
----------
dim : tuple
The dimensions of the grid.
center : tuple, optional
The center of the disc in pixel coordinates.
radius : float, optional
The radius of the disc in pixel coordinates.
height : float, optional
The height of the disc in pixel coordinates. Unused if only 2D.
axis : int, optional
The orientation of the discs orthogonal axis. Only used in 3D case with z-axis as default.
scale: tuple of float
Scaling along the dimensions of the underlying data.
"""
_log.debug('Calling __init__')
assert len(dim) in (2, 3), 'Disc can only be build in 2 or 3 dimensions!'
# Find indices of the disc plane axes:
idx_uv = [0, 1, 2]
if len(dim) == 3: # 3D:
idx_uv.remove(axis)
else: # 2D:
idx_uv.remove(2)
# Find default values:
if center is None:
center = tuple([d/2 for d in dim])
if radius is None:
radius = np.max((dim[idx_uv[0]], dim[idx_uv[1]])) / 4
if height is None and len(dim) == 3: # only used for 3D!
height = dim[axis] / 2
assert height > 0, 'Height has to be a positive scalar value!'
assert len(dim) == len(center), 'Parameters dim and center must have the same dimensions!'
assert radius > 0, 'Radius has to be a positive scalar value!'
coords = np.indices(dim) + 0.5 # 0.5 to get to pixel/voxel center!
coords = coords - np.asarray(center)[(slice(None),) + (None,)*len(dim)] # [:, None, None, None]
data = np.where(np.hypot(coords[idx_uv[0]], coords[idx_uv[1]]) <= radius, 1, 0)
if len(dim) == 3: # 3D: Implement bounds above and below the disc:
height_shape = np.zeros_like(data)
bounds = [slice(None)] * 3 # i.e.: [:, :, :]
start = int(np.floor(center[axis] - height/2))
stop = int(np.ceil(center[axis] + height/2))
bounds[axis] = slice(start, stop) # replace with actual bounds along disc symmetry axis
height_shape[tuple(bounds)] = 1
data *= height_shape
return Field(data=data, scale=scale, vector=False)
def create_shape_ellipse(dim, center=None, width=None, height=None, axis=0, scale=1.0):
"""Creates a Field object with the shape of an ellipse in 2D or 3D.
Attributes
----------
dim : tuple
The dimensions of the grid.
center : tuple, optional
The center of the ellipse in pixel coordinates.
width : tuple, optional
The two half axes of the ellipse in pixel coordinates.
height : float, optional
The height of the ellipse in pixel coordinates. Unused if only 2D.
axis : int, optional
The orientation of the ellipses orthogonal axis. Only used in 3D case with z-axis as default.
scale: tuple of float
Scaling along the dimensions of the underlying data.
"""
_log.debug('Calling __init__')
assert len(dim) in (2, 3), 'Ellipse can only be build in 2 or 3 dimensions!'
# Find indices of the disc plane axes:
idx_uv = [0, 1, 2]
if len(dim) == 3: # 3D:
idx_uv.remove(axis)
else: # 2D:
idx_uv.remove(2)
# Find default values:
if center is None:
center = tuple([d/2 for d in dim])
if width is None:
dim_uv = (dim[idx_uv[0]], dim[idx_uv[1]])
width = (np.max(dim_uv)*1/3, np.min(dim_uv)*2/3)
if height is None and len(dim) == 3: # only used for 3D!
height = dim[axis] / 2
assert len(dim) == len(center), 'Parameters dim and center must have the same dimensions!'
assert len(width) == 2, 'Width has to contain the two half axes!'
coords = np.indices(dim) + 0.5 # 0.5 to get to pixel/voxel center!
coords = coords - np.asarray(center)[(slice(None),) + (None,)*len(dim)] # [:, None, None, None]
distance = np.hypot(coords[idx_uv[1]] / (width[1]/2), coords[idx_uv[0]] / (width[0]/2))
data = np.where(distance <= 1, 1, 0)
if len(dim) == 3: # 3D: Implement bounds above and below the disc:
height_shape = np.zeros_like(data)
bounds = [slice(None)] * 3 # i.e.: [:, :, :]
start = int(np.floor(center[axis] - height/2))
stop = int(np.ceil(center[axis] + height/2))
bounds[axis] = slice(start, stop) # replace with actual bounds along disc symmetry axis
height_shape[tuple(bounds)] = 1
data *= height_shape
return Field(data=data, scale=scale, vector=False)
def create_shape_ellipsoid(dim, center=None, width=None, scale=1.0):
"""Creates a Field object with the shape of an ellipsoid in arbitrary dimensions.
Attributes
----------
dim : tuple
The dimensions of the grid.
center : tuple, optional
The center of the ellipsoid in pixel coordinates.
width : tuple, optional
The half axes of the ellipsoid in pixel coordinates.
scale: tuple of float
Scaling along the dimensions of the underlying data.
"""
_log.debug('Calling __init__')
if center is None:
center = tuple([d/2 for d in dim])
if width is None:
width = tuple([d/2 for d in dim])
assert len(dim) == len(center) == len(width), 'Parameters dim, center and width must have the same dimensions!'
coords = np.indices(dim) + 0.5 # 0.5 to get to pixel/voxel center!
coords = coords - np.asarray(center)[(slice(None),) + (None,)*len(dim)] # [:, None, None, None]
distance = np.sqrt([(coords[i] / (width[i]/2))**2 for i in range(len(dim))])
data = np.where(distance <= 1, 1, 0)
return Field(data=data, scale=scale, vector=False)
def create_shape_sphere(dim, center=None, radius=None, scale=1.0):
"""Creates a Field object with the shape of a sphere in arbitrary dimensions.
Attributes
----------
dim : tuple
The dimensions of the grid.
center : tuple, optional
The center of the sphere in pixel coordinates.
width : tuple, optional
The half axes of the sphere in pixel coordinates.
scale: tuple of float
Scaling along the dimensions of the underlying data.
"""
_log.debug('Calling __init__')
if center is None:
center = tuple([d/2 for d in dim])
if radius is None:
radius = np.max(dim) / 4
assert len(dim) == len(center), 'Parameters dim and center must have the same dimensions!'
assert radius > 0, 'Radius has to be a positive scalar value!'
coords = np.indices(dim) + 0.5 # 0.5 to get to pixel/voxel center!
coords = coords - np.asarray(center)[(slice(None),) + (None,)*len(dim)] # [:, None, None, None]
distance = np.sqrt([coords[i]**2 for i in range(len(dim))])
data = np.where(distance <= radius, 1, 0)
return Field(data=data, scale=scale, vector=False)
def create_shape_filament(dim, pos=None, axis=0, scale=1.0):
"""Creates a Field object with the shape of a filament in arbitrary dimension.
Parameters
----------
dim : tuple
The dimensions of the grid.
pos : tuple, optional
The position of the filament in pixel coordinates. Has to be a tuple of len(dim) - 1 and denotes the
index positions along all axes aside from the specified `axis` along which the filament is placed.
axis : int, optional
The orientation of the filament axis. Defaults to the first axis.
scale: tuple of float
Scaling along the dimensions of the underlying data.
"""
_log.debug('Calling __init__')
assert len(dim) > 1, 'Only usable for multidimensional Fields (at least 2 dimensions)!'
if pos is None:
pos = tuple([d/2 for d in dim])
assert len(dim) == len(pos)-1, 'Position has to specify one less coordinates than dimensions!'
data = np.zeros(dim)
index = list(pos)
index.insert(axis, slice(None)) # [:] for the filament axis! pos everywhere else!
data[tuple(index)] = 1
return Field(data=data, scale=scale, vector=False)
def create_shape_voxel(dim, pos=None, scale=1.0):
"""Creates a Field object with the shape of a single voxel in arbitrary dimension.
Parameters
----------
dim : tuple
The dimensions of the grid.
pos : tuple, optional
The position of the voxel.
scale: tuple of float
Scaling along the dimensions of the underlying data.
"""
_log.debug('Calling __init__')
assert len(dim) > 1, 'Only usable for multidimensional Fields (at least 2 dimensions)!'
if pos is None:
pos = tuple([d/2 for d in dim])
assert len(dim) == len(pos), 'Parameters dim and pos must have the same dimensions!'
data = np.zeros(dim)
data[pos] = 1
return Field(data=data, scale=scale, vector=False)
src/empyre/fields/vectors.py 0 → 100644
View file @ 6e3128b8
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
import logging
from numbers import Number
import numpy as np
from .field import Field
__all__ = ['create_vector_homog', 'create_vector_vortex', 'create_vector_skyrmion', 'create_vector_singularity']
_log = logging.getLogger(__name__)
def create_vector_homog(dim, phi=0, theta=None, scale=1.0):
"""Field subclass implementing a homogeneous vector field with 2 or 3 components in 2 or 3 dimensions.
Attributes
----------
dim : tuple
The dimensions of the grid.
phi : float
The azimuthal angle. The default is 0, i.e., the vectors point in x direction.
theta : float, optional
The polar angle. If None (default), only two components will be created (corresponds to pi/2 in 3D, i.e., the
vectors are in the xy plane).
scale: tuple of float
Scaling along the dimensions of the underlying data.
"""
_log.debug('Calling __init__')
assert len(dim) in (2, 3), 'Disc can only be build in 2 or 3 dimensions!'
assert isinstance(phi, Number), 'phi has to be an angle in radians!'
assert isinstance(theta, Number) or theta is None, 'theta has to be an angle in radians or None!'
if len(dim) == 2 and theta is None: # 2D and 3rd component not explicitely wanted:
y_comp = np.ones(dim) * np.sin(phi)
x_comp = np.ones(dim) * np.cos(phi)
data = np.stack([x_comp, y_comp], axis=-1)
else: # 3D, always have a third component:
if theta is None:
theta = np.pi/2 # xy-plane if not specified otherwise!
z_comp = np.ones(dim) * np.cos(theta)
y_comp = np.ones(dim) * np.sin(theta) * np.sin(phi)
x_comp = np.ones(dim) * np.sin(theta) * np.cos(phi)
data = np.stack([x_comp, y_comp, z_comp], axis=-1)
return Field(data=data, scale=scale, vector=True)
def create_vector_vortex(dim, center=None, phi_0=np.pi/2, oop_r=None, oop_sign=1, core_r=0, axis=0, scale=1.0):
"""Field subclass implementing a vortex vector field with 3 components in 2 or 3 dimensions.
Attributes
----------
dim : tuple
The dimensions of the grid.
center : tuple (N=2 or N=3), optional
The vortex center, given in 2D `(v, u)` or 3D `(z, y, x)`, where the perpendicular axis is discarded
(determined by the `axis` parameter). Is set to the center of the field of view if not specified.
The vortex center should be between two pixels to avoid singularities.
axis : int, optional
The orientation of the vortex axis. The default is 0 and corresponds to the z-axis. Is ignored if dim is
only 2D.
phi_0 : float, optional
Angular offset that allows all arrows to be rotated simultaneously. The default value is `pi/2` and corresponds
to an anti-clockwise rotation, while a value of `-pi/2` corresponds to a clockwise rotation. `0` and `pi` lead
to arrows that all point in/outward.
oop_r : float, optional
Radius of a potential out-of-plane ("oop") component in the vortex center. If `None`, the vortex is purely
in-plane, which is the default. Set this to a number (in pixels) that determines the radius in which the
oop-component is tilted into the plane. A `0` leads to an immediate/sharp tilt (not visible without setting a
core radius, see `core_r`), while setting this to the radius of your vortex disc, will let it tilt smoothly
until reaching the edge. If this is `None`, only two components will be created for 2-dimensional vortices!
oop_sign : {1, -1}, optional
Only used if `oop_r` is not None (i.e. there is an out-of-plane component). Can be `1` (default) or `-1` and
determines if the core (if it exists) is aligned parallel (`+1`) or antiparallel (`-1`) to the chosen symmetry
axis.
core_r : float, optional
Radius of a potential vortex core that's homogeneously oriented out-of-plane. Is not used when `oop_r` is not
set and defaults to `0`, which means the vortex core is infinitely small.
scale: tuple of float
Scaling along the dimensions of the underlying data.
Returns
-------
field : `~.Field` object
The resulting vector field.
"""
_log.debug('Calling create_vector_vortex')
assert len(dim) in (2, 3), 'Vortex can only be build in 2 or 3 dimensions!'
# Find indices of the vortex plane axes:
idx_uv = [0, 1, 2]
if len(dim) == 3: # 3D:
idx_uv.remove(axis)
else: # 2D:
idx_uv.remove(2)
# 2D dimensions:
dim_uv = tuple([dim[i] for i in idx_uv])
# Find default values:
if center is None:
center = tuple([dim[i] / 2 for i in idx_uv])
elif len(center) == 3: # if a 3D-center is given, just take the relevant coordinates:
center = list(center)
del center[axis]
center = tuple(center)
# Create vortex plane (2D):
coords_uv = np.indices(dim_uv) + 0.5 # 0.5 to get to pixel/voxel center!
coords_uv = coords_uv - np.asarray(center, dtype=float)[:, None, None] # Shift by center!
phi = np.arctan2(coords_uv[0], coords_uv[1]) - phi_0
rr = np.hypot(coords_uv[0], coords_uv[1])
rr_clip = np.clip(rr - core_r, a_min=0, a_max=None) # negative inside core_r (clipped to 0), positive outside!
# Check if a core should be constructed (oop_r != 0):
if oop_r is None:
w_comp = np.zeros(dim_uv)
else:
assert oop_r >= 0, 'oop_r has to be a positive number!'
assert oop_sign in (1, -1), 'Sign of the out-of-plane component has to be either +1 or -1!'
w_comp = 1 - 2/np.pi * np.arcsin(np.tanh(np.pi*rr_clip/(oop_r+1E-30))) # orthogonal: 1 inside, to 0 outside!
w_comp *= oop_sign * w_comp
v_comp = np.ones(dim_uv) * np.sin(phi) * np.sqrt(1 - np.abs(w_comp))
u_comp = np.ones(dim_uv) * np.cos(phi) * np.sqrt(1 - np.abs(w_comp))
if len(dim) == 3: # Expand to 3D:
w_comp = np.expand_dims(w_comp, axis=axis)
v_comp = np.expand_dims(v_comp, axis=axis)
u_comp = np.expand_dims(u_comp, axis=axis)
reps = [1, 1, 1] # repetitions for tiling
reps[axis] = dim[axis] # repeat along chosen axis
w_comp = np.tile(w_comp, reps)
v_comp = np.tile(v_comp, reps)
u_comp = np.tile(u_comp, reps)
if axis == 0: # z-axis
z_comp = w_comp
y_comp = -v_comp
x_comp = -u_comp
elif axis == 1: # y-axis
z_comp = v_comp
y_comp = w_comp
x_comp = u_comp
elif axis == 2: # x-axis
z_comp = -v_comp
y_comp = -u_comp
x_comp = w_comp
else:
raise ValueError(f'{axis} is not a valid argument for axis (has to be 0, 1 or 2)')
data = np.stack([x_comp, y_comp, z_comp], axis=-1)
return Field(data=data, scale=scale, vector=True)
def create_vector_skyrmion(dim, center=None, phi_0=0, core_sign=1, skyrm_d=None, wall_d=None, axis=0, scale=1.0):
"""Create a 3-dimensional magnetic Bloch or Neel type skyrmion distribution.
Parameters
----------
dim : tuple
The dimensions of the grid.
center : tuple (N=2 or N=3), optional
The source center, given in 2D `(v, u)` or 3D `(z, y, x)`, where the perpendicular axis
is discarded. Is set to the center of the field of view if not specified.
The center has to be between two pixels.
phi_0 : float, optional
Angular offset switching between Neel type (0 [default] or pi) or Bloch type (+/- pi/2)
skyrmions.
core_sign : {1, -1}, optional
Can be `1` (default) or `-1` and determines if the skyrmion core is aligned parallel (`+1`) or antiparallel
(`-1`) to the chosen symmetry axis.
skyrm_d : float, optional
Diameter of the skyrmion. Defaults to half of the smaller dimension perpendicular to the
skyrmion axis if not specified.
wall_d : float, optional
Diameter of the domain wall of the skyrmion. Defaults to `skyrm_d / 4` if not specified.
axis : int, optional
The orientation of the vortex axis. The default is 0 and corresponds to the z-axis. Is ignored if dim is
only 2D.
scale: tuple of float
Scaling along the dimensions of the underlying data.
Returns
-------
field : `~.Field` object
The resulting vector field.
Notes
-----
To avoid singularities, the source center should lie between the pixel centers (which
reside at coordinates with _.5 at the end), i.e. integer values should be used as center
coordinates (e.g. coordinate 1 lies between the first and the second pixel).
Skyrmion wall width is dependant on exchange stiffness A [J/m] and anisotropy K [J/m³]
The out-of-plane magnetization at the domain wall is described in a paper by Romming et al
(see DOI: 10.1103/PhysRevLett.114.177203s)
"""
def _theta(r):
theta_1 = + np.arcsin(np.tanh((r + skyrm_d/2)/(wall_d/2)))
theta_2 = - np.arcsin(np.tanh((r - skyrm_d/2)/(wall_d/2)))
theta = theta_1 + theta_2
theta /= np.abs(theta).max() / np.pi
return theta
_log.debug('Calling create_vector_skyrmion')
assert len(dim) in (2, 3), 'Skyrmion can only be build in 2 or 3 dimensions!'
assert core_sign in (1, -1), 'Sign of the out-of-plane component has to be either +1 or -1!'
# Find indices of the skyrmion plane axes:
idx_uv = [0, 1, 2]
if len(dim) == 3: # 3D:
idx_uv.remove(axis)
else: # 2D:
idx_uv.remove(2)
# 2D dimensions:
dim_uv = tuple([dim[i] for i in idx_uv])
# Find default values:
if skyrm_d is None:
skyrm_d = np.min(dim_uv) / 2
if wall_d is None:
wall_d = skyrm_d / 4
if center is None:
center = tuple([dim[i] / 2 for i in idx_uv])
elif len(center) == 3: # if a 3D-center is given, just take the relevant coordinates:
center = list(center)
del center[axis]
center = tuple(center)
# Create skyrmion plane (2D):
coords_uv = np.indices(dim_uv) + 0.5 # 0.5 to get to pixel/voxel center!
coords_uv = coords_uv - np.asarray(center, dtype=float)[:, None, None] # Shift by center!
rr = np.hypot(coords_uv[0], coords_uv[1])
phi = np.arctan2(coords_uv[0], coords_uv[1]) - phi_0
theta = _theta(rr)
w_comp = core_sign * np.cos(theta)
v_comp = np.sin(theta) * np.sin(phi)
u_comp = np.sin(theta) * np.cos(phi)
# Expansion to 3D if necessary and component shuffling:
if len(dim) == 3: # Expand to 3D:
w_comp = np.expand_dims(w_comp, axis=axis)
v_comp = np.expand_dims(v_comp, axis=axis)
u_comp = np.expand_dims(u_comp, axis=axis)
reps = [1, 1, 1] # repetitions for tiling
reps[axis] = dim[axis] # repeat along chosen axis
w_comp = np.tile(w_comp, reps)
v_comp = np.tile(v_comp, reps)
u_comp = np.tile(u_comp, reps)
if axis == 0: # z-axis
z_comp = w_comp
y_comp = -v_comp
x_comp = -u_comp
elif axis == 1: # y-axis
z_comp = v_comp
y_comp = w_comp
x_comp = u_comp
elif axis == 2: # x-axis
z_comp = -v_comp
y_comp = -u_comp
x_comp = w_comp
else:
raise ValueError(f'{axis} is not a valid argument for axis (has to be 0, 1 or 2)')
data = np.stack([x_comp, y_comp, z_comp], axis=-1)
return Field(data=data, scale=scale, vector=True)
def create_vector_singularity(dim, center=None, scale=1.0):
"""Create a 3-dimensional magnetic distribution of a homogeneous magnetized object.
Parameters
----------
dim : tuple
The dimensions of the grid.
center : tuple (N=2 or N=3), optional
The source center, given in 2D `(v, u)` or 3D `(z, y, x)`, where the perpendicular axis
is discarded. Is set to the center of the field of view if not specified.
The source center has to be between two pixels.
scale: tuple of float
Scaling along the dimensions of the underlying data.
Returns
-------
field : `~.Field` object
The resulting vector field.
Notes
-----
To avoid singularities, the source center should lie between the pixel centers (which
reside at coordinates with _.5 at the end), i.e. integer values should be used as center
coordinates (e.g. coordinate 1 lies between the first and the second pixel).
Per default, all arrows point outwards, negate the resulting `Field` object with a minus after creation to let
them point inwards.
"""
_log.debug('Calling create_vector_singularity')
# Find default values:
if center is None:
center = tuple([d / 2 for d in dim])
assert len(dim) == len(center), f"Length of dim ({len(dim)}) and center ({len(center)}) don't match!"
# Setup coordinates, shape is (c, z, y, x), if 3D, or (c, y, x), if 2D (c: components):
coords = np.indices(dim) + 0.5 # 0.5 to get to pixel/voxel center!
bc_shape = (len(dim,),) + (1,)*len(dim) # Shape for broadcasting, (3,1,1,1) for 3D, (2,1,1) for 2D!
coords = coords - np.reshape(center, bc_shape) # Shift by center (append 1s for broadcasting)!
rr = np.sqrt(np.sum([coords[i]**2 for i in range(len(dim))], axis=0))
coords /= rr + 1E-30 # Normalise amplitude (keep direction), rr (z,y,x) is broadcasted to data (c,z,y,x)!
# coords has components listed in wrong order (z, y, x) in the first dimension, we need (x, y, z) in the last:
coords = np.flip(coords, axis=0)
# Finally, the data has the coordinate axis at the end, not at the front:
data = np.moveaxis(coords, 0, -1) # (c,z,y,x) -> (z,y,x,c)
return Field(data=data, scale=scale, vector=True)
src/empyre/io/__init__.py 0 → 100644
View file @ 6e3128b8
# -*- coding: utf-8 -*-
# Copyright 2020 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Subpackage containing EMPyRe IO functionality for several EMPyRe classes."""
from .io_field import *
__all__ = []
__all__.extend(io_field.__all__)
del io_field