From 1430e331a7d8ddb2e0a449a29aff49fe4a8cc7eb Mon Sep 17 00:00:00 2001 From: Jan Caron <j.caron@fz-juelich.de> Date: Mon, 1 Jul 2013 10:04:57 +0200 Subject: [PATCH] Minor changes in some modules, splitted "compare_method_erros" script holoimage: densitiy default set to 1 --- pyramid/holoimage.py | 4 +- scripts/compare_method_errors.py | 192 ++++++++---- scripts/compare_method_errors_res.py | 423 +++++++++++++++++++++++++++ scripts/compare_methods.py | 4 +- scripts/create_multiple_samples.py | 62 ++++ scripts/reconstruct_random_pixels.py | 7 +- 6 files changed, 632 insertions(+), 60 deletions(-) create mode 100644 scripts/compare_method_errors_res.py create mode 100644 scripts/create_multiple_samples.py diff --git a/pyramid/holoimage.py b/pyramid/holoimage.py index 0819527..2aeb2c0 100644 --- a/pyramid/holoimage.py +++ b/pyramid/holoimage.py @@ -31,11 +31,11 @@ CDICT = {'red': [(0.00, 1.0, 0.0), HOLO_CMAP = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT, 256) -def holo_image(phase_map, density): +def holo_image(phase_map, density=1): '''Returns a holography image with color-encoded gradient direction. Arguments: phase_map - a PhaseMap object storing the phase informations - density - the factor for determining the number of contour lines + density - the gain factor for determining the number of contour lines (default: 1) Returns: holography image diff --git a/scripts/compare_method_errors.py b/scripts/compare_method_errors.py index ce967cd..0c451b1 100644 --- a/scripts/compare_method_errors.py +++ b/scripts/compare_method_errors.py @@ -11,7 +11,6 @@ import matplotlib.pyplot as plt import pyramid.magcreator as mc import pyramid.projector as pj import pyramid.phasemapper as pm -import pyramid.holoimage as hi import pyramid.analytic as an from pyramid.magdata import MagData from pyramid.phasemap import PhaseMap @@ -88,58 +87,145 @@ def phase_from_mag(): axis.set_ylabel('duration [s]') -# fig = plt.figure() -# axis = fig.add_subplot(1, 1, 1, aspect='equal') -# im = plt.pcolormesh(self.phase, cmap=cmap) -# # Set the axes ticks and labels: -# ticks = axis.get_xticks()*self.res -# axis.set_xticklabels(ticks) -# ticks = axis.get_yticks()*self.res -# axis.set_yticklabels(ticks) -# axis.set_title(title) -# axis.set_xlabel('x-axis [nm]') -# axis.set_ylabel('y-axis [nm]') -# # Plot the phase map: -# fig = plt.gcf() -# fig.subplots_adjust(right=0.85) -# cbar_ax = fig.add_axes([0.9, 0.15, 0.02, 0.7]) -# fig.colorbar(im, cax=cbar_ax) -# plt.show() - -# '''REAL SLAB ''' -# phi = np.linspace(0, 2*pi, endpoint=False, num=72) -# -# RMS = np.zeros(len(phi)) -# for i in range(len(phi)): -# print 'phi =', round(360*phi[i]/(2*pi)) -# mag_data = MagData(res, mc.create_mag_dist(mag_shape, phi[i])) -# projection = pj.simple_axis_projection(mag_data) -# phase_num = pm.phase_mag_real(res, projection, 'slab', b_0) -# phase_ana = an.phase_mag_slab(dim, res, phi[i], center, width, b_0) -# phase_diff = phase_ana - phase_num -# RMS[i] = np.std(phase_diff) -# -# fig = plt.figure() -# fig.add_subplot(111) -# plt.plot(np.round(360*phi/(2*pi)), RMS) - - -# phase_map_slab = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab', b_0)) -# phase_map_disc = PhaseMap(res, pm.phase_mag_real(res, projection, 'disc', b_0)) -# # Display the combinated plots with phasemap and holography image: -# hi.display_combined(phase_map_ana, density, 'Analytic Solution') -# hi.display_combined(phase_map_fft, density, 'Fourier Space') -# hi.display_combined(phase_map_slab, density, 'Real Space (Slab)') -# hi.display_combined(phase_map_disc, density, 'Real Space (Disc)') -# -# # Plot differences to the analytic solution: -# -# phase_map_diff_slab = PhaseMap(res, phase_map_ana.phase-phase_map_slab.phase) -# phase_map_diff_disc = PhaseMap(res, phase_map_ana.phase-phase_map_disc.phase) -# -# RMS_slab = phase_map_diff_slab.phase -# RMS_disc = phase_map_diff_disc.phase - + + + + + + '''VARY DIMENSIONS FOR ALL APPROACHES''' + + b_0 = 1 # in T + phi = -pi/4 + dim_list = [(1, 4, 4), (1, 8, 8), (1, 16, 16), (1, 32, 32), (1, 64, 64), (1, 128, 128)] + res_list = [64., 32., 16., 8., 4., 2., 1.] # in nm + + + + data_sl_p_fourier0 = np.zeros((3, len(res_list))) + data_sl_w_fourier0 = np.zeros((3, len(res_list))) + data_disc_fourier0 = np.zeros((3, len(res_list))) + + data_sl_p_fourier1 = np.zeros((3, len(res_list))) + data_sl_w_fourier1 = np.zeros((3, len(res_list))) + data_disc_fourier1 = np.zeros((3, len(res_list))) + + data_sl_p_fourier20 = np.zeros((3, len(res_list))) + data_sl_w_fourier20 = np.zeros((3, len(res_list))) + data_disc_fourier20 = np.zeros((3, len(res_list))) + + data_sl_p_real_s = np.zeros((3, len(res_list))) + data_sl_w_real_s = np.zeros((3, len(res_list))) + data_disc_real_s = np.zeros((3, len(res_list))) + + data_sl_p_real_d= np.zeros((3, len(res_list))) + data_sl_w_real_d = np.zeros((3, len(res_list))) + data_disc_real_d = np.zeros((3, len(res_list))) + + data_slab_perfect[0, :] = res_list + data_slab_worst[0, :] = res_list + data_disc[0, :] = res_list + + + + '''FOURIER UNPADDED''' + + for i, (dim, res) in enumerate(zip(dim_list, res_list)): + + print 'dim =', str(dim) + + # ANALYTIC SOLUTIONS: + # Slab (perfectly aligned): + center = (0, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts with 0! + width = (1, dim[1], dim[2]) # in px (z, y, x) + mag_shape_slab_perfect = mc.Shapes.slab(dim, center, width) + phase_ana_slab_perfect = an.phase_mag_slab(dim, res, phi, center, width, b_0) + mag_data_slab_perfect = MagData(res, mc.create_mag_dist(mag_shape_slab_perfect, phi)) + projection_slab_perfect = pj.simple_axis_projection(mag_data_slab_perfect) + # Slab (worst case): + center = (0, dim[1]/2, dim[2]/2) # in px (z, y, x) index starts with 0! + width = (1, dim[1], dim[2]) # in px (z, y, x) + mag_shape_slab_worst = mc.Shapes.slab(dim, center, width) + phase_ana_slab_worst = an.phase_mag_slab(dim, res, phi, center, width, b_0) + mag_data_slab_worst = MagData(res, mc.create_mag_dist(mag_shape_slab_worst, phi)) + projection_slab_worst = pj.simple_axis_projection(mag_data_slab_worst) + # Disc: + center = (0, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts with 0! + radius = dim[1]/2 # in px + height = 1 # in px + mag_shape_disc = mc.Shapes.disc(dim, center, radius, height) + phase_ana_disc = an.phase_mag_disc(dim, res, phi, center, radius, height, b_0) + mag_data_disc = MagData(res, mc.create_mag_dist(mag_shape_disc, phi)) + projection_disc = pj.simple_axis_projection(mag_data_disc) + + # NUMERICAL SOLUTIONS: + # Slab (perfectly aligned): + key = ', '.join(['Resolution->RMS|duration', 'Fourier', 'padding=0', + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=slab_perfect']) + if data_shelve.has_key(key): + data_slab_perfect[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_fourier(res, projection_slab_perfect, b_0, 0) + data_slab_perfect[2, i] = time.time() - start_time + phase_diff = phase_ana_slab_perfect - phase_num + data_slab_perfect[1, i] = np.std(phase_diff) + data_shelve[key] = data_slab_perfect[:, i] + # Slab (worst case): + key = ', '.join(['Resolution->RMS|duration', 'Fourier', 'padding=0', + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=slab_worst']) + if data_shelve.has_key(key): + data_slab_worst[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_fourier(res, projection_slab_worst, b_0, 0) + data_slab_worst[2, i] = time.time() - start_time + phase_diff = phase_ana_slab_worst - phase_num + data_slab_worst[1, i] = np.std(phase_diff) + data_shelve[key] = data_slab_worst[:, i] + # Disc: + key = ', '.join(['Resolution->RMS|duration', 'Fourier', 'padding=0', + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=disc']) + if data_shelve.has_key(key): + data_disc[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_fourier(res, projection_disc, b_0, 0) + data_disc[2, i] = time.time() - start_time + phase_diff = phase_ana_disc - phase_num + data_disc[1, i] = np.std(phase_diff) + data_shelve[key] = data_disc[:, i] + + # Plot duration against res: + fig = plt.figure() + axis = fig.add_subplot(1, 1, 1) + axis.set_yscale('log') + axis.plot(data_slab_perfect[0], data_slab_perfect[1], + data_slab_worst[0], data_slab_worst[1], + data_disc[0], data_disc[1]) + axis.set_title('Variation of the resolution (Fourier without padding)') + axis.set_xlabel('res [nm]') + axis.set_ylabel('RMS') + # Plot RMS against res: + fig = plt.figure() + axis = fig.add_subplot(1, 1, 1) + axis.plot(data_slab_perfect[0], data_slab_perfect[1], + data_slab_worst[0], data_slab_worst[1], + data_disc[0], data_disc[1]) + axis.set_title('Variation of the resolution (Fourier without padding)') + axis.set_xlabel('res [nm]') + axis.set_ylabel('duration [s]') + + + + + + + + + data_shelve.close() diff --git a/scripts/compare_method_errors_res.py b/scripts/compare_method_errors_res.py new file mode 100644 index 0000000..6fcc589 --- /dev/null +++ b/scripts/compare_method_errors_res.py @@ -0,0 +1,423 @@ +# -*- coding: utf-8 -*- +"""Compare the different methods to create phase maps.""" + + +import time +import pdb, traceback, sys +import numpy as np +from numpy import pi +import matplotlib.pyplot as plt + +import pyramid.magcreator as mc +import pyramid.projector as pj +import pyramid.phasemapper as pm +import pyramid.analytic as an +from pyramid.magdata import MagData +import shelve + + +def compare_method_errors_res(): + '''Calculate and display the phase map from a given magnetization. + Arguments: + None + Returns: + None + + ''' + # Create / Open databank: + data_shelve = shelve.open('../output/method_errors_shelve') + + + + + + + + '''VARY DIMENSIONS FOR ALL APPROACHES''' + + b_0 = 1 # in T + phi = -pi/4 + dim_list = [(1, 4, 4), (1, 8, 8), (1, 16, 16), (1, 32, 32), (1, 64, 64), (1, 128, 128), (1, 256, 256), (1, 512, 512)] + res_list = [64., 32., 16., 8., 4., 2., 1., 0.5, 0.25] # in nm + + + '''CREATE DATA ARRAYS''' + + data_sl_p_fourier0 = np.zeros((3, len(res_list))) + data_sl_w_fourier0 = np.zeros((3, len(res_list))) + data_disc_fourier0 = np.zeros((3, len(res_list))) + + data_sl_p_fourier1 = np.zeros((3, len(res_list))) + data_sl_w_fourier1 = np.zeros((3, len(res_list))) + data_disc_fourier1 = np.zeros((3, len(res_list))) + + data_sl_p_fourier20 = np.zeros((3, len(res_list))) + data_sl_w_fourier20 = np.zeros((3, len(res_list))) + data_disc_fourier20 = np.zeros((3, len(res_list))) + + data_sl_p_real_s = np.zeros((3, len(res_list))) + data_sl_w_real_s = np.zeros((3, len(res_list))) + data_disc_real_s = np.zeros((3, len(res_list))) + + data_sl_p_real_d= np.zeros((3, len(res_list))) + data_sl_w_real_d = np.zeros((3, len(res_list))) + data_disc_real_d = np.zeros((3, len(res_list))) + + + '''CREATE DATA ARRAYS''' + + data_sl_p_fourier0[0, :] = res_list + data_sl_w_fourier0[0, :] = res_list + data_disc_fourier0[0, :] = res_list + + data_sl_p_fourier1[0, :] = res_list + data_sl_w_fourier1[0, :] = res_list + data_disc_fourier1[0, :] = res_list + + data_sl_p_fourier20[0, :] = res_list + data_sl_w_fourier20[0, :] = res_list + data_disc_fourier20[0, :] = res_list + + data_sl_p_real_s[0, :] = res_list + data_sl_w_real_s[0, :] = res_list + data_disc_real_s[0, :] = res_list + + data_sl_p_real_d[0, :]= res_list + data_sl_w_real_d[0, :] = res_list + data_disc_real_d[0, :] = res_list + + + + + for i, (dim, res) in enumerate(zip(dim_list, res_list)): + + print 'i =', i, ' dim =', str(dim) + + '''ANALYTIC SOLUTIONS''' + + # Slab (perfectly aligned): + center = (0, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts with 0! + width = (1, dim[1]/2, dim[2]/2) # in px (z, y, x) + mag_shape_sl_p = mc.Shapes.slab(dim, center, width) + phase_ana_sl_p = an.phase_mag_slab(dim, res, phi, center, width, b_0) + mag_data_sl_p = MagData(res, mc.create_mag_dist(mag_shape_sl_p, phi)) + projection_sl_p = pj.simple_axis_projection(mag_data_sl_p) + # Slab (worst case): + center = (0, dim[1]/2, dim[2]/2) # in px (z, y, x) index starts with 0! + width = (1, dim[1]/2, dim[2]/2) # in px (z, y, x) + mag_shape_sl_w = mc.Shapes.slab(dim, center, width) + phase_ana_sl_w = an.phase_mag_slab(dim, res, phi, center, width, b_0) + mag_data_sl_w = MagData(res, mc.create_mag_dist(mag_shape_sl_w, phi)) + projection_sl_w = pj.simple_axis_projection(mag_data_sl_w) + # Disc: + center = (0, 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 = 1 # in px + mag_shape_disc = mc.Shapes.disc(dim, center, radius, height) + phase_ana_disc = an.phase_mag_disc(dim, res, phi, center, radius, height, b_0) + mag_data_disc = MagData(res, mc.create_mag_dist(mag_shape_disc, phi)) + projection_disc = pj.simple_axis_projection(mag_data_disc) + + '''FOURIER UNPADDED''' + padding = 0 + # Slab (perfectly aligned): + key = ', '.join(['Resolution->RMS|duration', 'Fourier', 'padding={}'.format(padding), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=sl_p']) + if data_shelve.has_key(key): + data_sl_p_fourier0[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_fourier(res, projection_sl_p, b_0, padding) + data_sl_p_fourier0[2, i] = time.time() - start_time + phase_diff = phase_ana_sl_p - phase_num + data_sl_p_fourier0[1, i] = np.std(phase_diff) + data_shelve[key] = data_sl_p_fourier0[:, i] + # Slab (worst case): + key = ', '.join(['Resolution->RMS|duration', 'Fourier', 'padding={}'.format(padding), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=sl_w']) + if data_shelve.has_key(key): + data_sl_w_fourier0[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_fourier(res, projection_sl_w, b_0, padding) + data_sl_w_fourier0[2, i] = time.time() - start_time + phase_diff = phase_ana_sl_w - phase_num + data_sl_w_fourier0[1, i] = np.std(phase_diff) + data_shelve[key] = data_sl_w_fourier0[:, i] + # Disc: + key = ', '.join(['Resolution->RMS|duration', 'Fourier', 'padding={}'.format(padding), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=disc']) + if data_shelve.has_key(key): + data_disc_fourier0[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_fourier(res, projection_disc, b_0, padding) + data_disc_fourier0[2, i] = time.time() - start_time + phase_diff = phase_ana_disc - phase_num + data_disc_fourier0[1, i] = np.std(phase_diff) + data_shelve[key] = data_disc_fourier0[:, i] + + + '''FOURIER PADDED ONCE''' + padding = 1 + # Slab (perfectly aligned): + key = ', '.join(['Resolution->RMS|duration', 'Fourier', 'padding={}'.format(padding), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=sl_p']) + if data_shelve.has_key(key): + data_sl_p_fourier1[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_fourier(res, projection_sl_p, b_0, padding) + data_sl_p_fourier1[2, i] = time.time() - start_time + phase_diff = phase_ana_sl_p - phase_num + data_sl_p_fourier1[1, i] = np.std(phase_diff) + data_shelve[key] = data_sl_p_fourier1[:, i] + # Slab (worst case): + key = ', '.join(['Resolution->RMS|duration', 'Fourier', 'padding={}'.format(padding), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=sl_w']) + if data_shelve.has_key(key): + data_sl_w_fourier1[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_fourier(res, projection_sl_w, b_0, padding) + data_sl_w_fourier1[2, i] = time.time() - start_time + phase_diff = phase_ana_sl_w - phase_num + data_sl_w_fourier1[1, i] = np.std(phase_diff) + data_shelve[key] = data_sl_w_fourier1[:, i] + # Disc: + key = ', '.join(['Resolution->RMS|duration', 'Fourier', 'padding={}'.format(padding), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=disc']) + if data_shelve.has_key(key): + data_disc_fourier1[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_fourier(res, projection_disc, b_0, padding) + data_disc_fourier1[2, i] = time.time() - start_time + phase_diff = phase_ana_disc - phase_num + data_disc_fourier1[1, i] = np.std(phase_diff) + data_shelve[key] = data_disc_fourier1[:, i] + + + '''FOURIER PADDED 20''' + if dim[1] <= 128: + padding = 20 + # Slab (perfectly aligned): + key = ', '.join(['Resolution->RMS|duration', 'Fourier', 'padding={}'.format(padding), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=sl_p']) + if data_shelve.has_key(key): + data_sl_p_fourier20[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_fourier(res, projection_sl_p, b_0, padding) + data_sl_p_fourier20[2, i] = time.time() - start_time + phase_diff = phase_ana_sl_p - phase_num + data_sl_p_fourier20[1, i] = np.std(phase_diff) + data_shelve[key] = data_sl_p_fourier20[:, i] + # Slab (worst case): + key = ', '.join(['Resolution->RMS|duration', 'Fourier', 'padding={}'.format(padding), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=sl_w']) + if data_shelve.has_key(key): + data_sl_w_fourier20[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_fourier(res, projection_sl_w, b_0, padding) + data_sl_w_fourier20[2, i] = time.time() - start_time + phase_diff = phase_ana_sl_w - phase_num + data_sl_w_fourier20[1, i] = np.std(phase_diff) + data_shelve[key] = data_sl_w_fourier20[:, i] + # Disc: + key = ', '.join(['Resolution->RMS|duration', 'Fourier', 'padding={}'.format(padding), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=disc']) + if data_shelve.has_key(key): + data_disc_fourier20[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_fourier(res, projection_disc, b_0, padding) + data_disc_fourier20[2, i] = time.time() - start_time + phase_diff = phase_ana_disc - phase_num + data_disc_fourier20[1, i] = np.std(phase_diff) + data_shelve[key] = data_disc_fourier20[:, i] + + + '''REAL SLAB''' + method = 'slab' + # Slab (perfectly aligned): + key = ', '.join(['Resolution->RMS|duration', 'Real', 'method={}'.format(method), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=sl_p']) + if data_shelve.has_key(key): + data_sl_p_real_s[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_real(res, projection_sl_p, method, b_0) + data_sl_p_real_s[2, i] = time.time() - start_time + phase_diff = phase_ana_sl_p - phase_num + data_sl_p_real_s[1, i] = np.std(phase_diff) + data_shelve[key] = data_sl_p_real_s[:, i] + # Slab (worst case): + key = ', '.join(['Resolution->RMS|duration', 'Real', 'method={}'.format(method), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=sl_w']) + if data_shelve.has_key(key): + data_sl_w_real_s[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_real(res, projection_sl_w, method, b_0) + data_sl_w_real_s[2, i] = time.time() - start_time + phase_diff = phase_ana_sl_w - phase_num + data_sl_w_real_s[1, i] = np.std(phase_diff) + data_shelve[key] = data_sl_w_real_s[:, i] + # Disc: + key = ', '.join(['Resolution->RMS|duration', 'Real', 'method={}'.format(method), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=disc']) + if data_shelve.has_key(key): + data_disc_real_s[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_real(res, projection_disc, method, b_0) + data_disc_real_s[2, i] = time.time() - start_time + phase_diff = phase_ana_disc - phase_num + data_disc_real_s[1, i] = np.std(phase_diff) + data_shelve[key] = data_disc_real_s[:, i] + + + '''REAL DISC''' + method = 'disc' + # Slab (perfectly aligned): + key = ', '.join(['Resolution->RMS|duration', 'Real', 'method={}'.format(method), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=sl_p']) + if data_shelve.has_key(key): + data_sl_p_real_d[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_real(res, projection_sl_p, method, b_0) + data_sl_p_real_d[2, i] = time.time() - start_time + phase_diff = phase_ana_sl_p - phase_num + data_sl_p_real_d[1, i] = np.std(phase_diff) + data_shelve[key] = data_sl_p_real_d[:, i] + # Slab (worst case): + key = ', '.join(['Resolution->RMS|duration', 'Real', 'method={}'.format(method), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=sl_w']) + if data_shelve.has_key(key): + data_sl_w_real_d[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_real(res, projection_sl_w, method, b_0) + data_sl_w_real_d[2, i] = time.time() - start_time + phase_diff = phase_ana_sl_w - phase_num + data_sl_w_real_d[1, i] = np.std(phase_diff) + data_shelve[key] = data_sl_w_real_d[:, i] + # Disc: + key = ', '.join(['Resolution->RMS|duration', 'Real', 'method={}'.format(method), + 'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim), + 'phi={}'.format(phi), 'geometry=disc']) + if data_shelve.has_key(key): + data_disc_real_d[:, i] = data_shelve[key] + else: + start_time = time.time() + phase_num = pm.phase_mag_real(res, projection_disc, method, b_0) + data_disc_real_d[2, i] = time.time() - start_time + phase_diff = phase_ana_disc - phase_num + data_disc_real_d[1, i] = np.std(phase_diff) + data_shelve[key] = data_disc_real_d[:, i] + + + # Plot duration against res (perfect slab): + fig = plt.figure() + axis = fig.add_subplot(1, 1, 1) + axis.set_yscale('log') + axis.plot(data_sl_p_fourier0[0], data_sl_p_fourier0[1], 'b+', + data_sl_p_fourier1[0], data_sl_p_fourier1[1], 'bx', + data_sl_p_fourier20[0], data_sl_p_fourier20[1], 'b*', + data_sl_p_real_s[0], data_sl_p_real_s[1], 'rs', + data_sl_p_real_d[0], data_sl_p_real_d[1], 'ro') + axis.set_title('Variation of the resolution (perfectly adjusted slab)') + axis.set_xlabel('res [nm]') + axis.set_ylabel('RMS') + # Plot RMS against res (perfect slab): + fig = plt.figure() + axis = fig.add_subplot(1, 1, 1) + axis.set_yscale('log') + axis.plot(data_sl_p_fourier0[0], data_sl_p_fourier0[2], 'b+', + data_sl_p_fourier1[0], data_sl_p_fourier1[2], 'bx', + data_sl_p_fourier20[0], data_sl_p_fourier20[2], 'b*', + data_sl_p_real_s[0], data_sl_p_real_s[2], 'rs', + data_sl_p_real_d[0], data_sl_p_real_d[2], 'ro') + axis.set_title('Variation of the resolution (perfectly adjusted slab)') + axis.set_xlabel('res [nm]') + axis.set_ylabel('duration [s]') + + # Plot duration against res (worst case slab): + fig = plt.figure() + axis = fig.add_subplot(1, 1, 1) + axis.set_yscale('log') + axis.plot(data_sl_w_fourier0[0], data_sl_w_fourier0[1], 'b+', + data_sl_w_fourier1[0], data_sl_w_fourier1[1], 'bx', + data_sl_w_fourier20[0], data_sl_w_fourier20[1], 'b*', + data_sl_w_real_s[0], data_sl_w_real_s[1], 'rs', + data_sl_w_real_d[0], data_sl_w_real_d[1], 'ro') + axis.set_title('Variation of the resolution (worst case slab)') + axis.set_xlabel('res [nm]') + axis.set_ylabel('RMS') + # Plot RMS against res (worst case slab): + fig = plt.figure() + axis = fig.add_subplot(1, 1, 1) + axis.set_yscale('log') + axis.plot(data_sl_w_fourier0[0], data_sl_w_fourier0[2], 'b+', + data_sl_w_fourier1[0], data_sl_w_fourier1[2], 'bx', + data_sl_w_fourier20[0], data_sl_w_fourier20[2], 'b*', + data_sl_w_real_s[0], data_sl_w_real_s[2], 'rs', + data_sl_w_real_d[0], data_sl_w_real_d[2], 'ro') + axis.set_title('Variation of the resolution (worst case slab)') + axis.set_xlabel('res [nm]') + axis.set_ylabel('duration [s]') + + # Plot duration against res (disc<): + fig = plt.figure() + axis = fig.add_subplot(1, 1, 1) + axis.set_yscale('log') + axis.plot(data_disc_fourier0[0], data_disc_fourier0[1], 'b+', + data_disc_fourier1[0], data_disc_fourier1[1], 'bx', + data_disc_fourier20[0], data_disc_fourier20[1], 'b*', + data_disc_real_s[0], data_disc_real_s[1], 'rs', + data_disc_real_d[0], data_disc_real_d[1], 'ro') + axis.set_title('Variation of the resolution (disc)') + axis.set_xlabel('res [nm]') + axis.set_ylabel('RMS') + # Plot RMS against res (disc): + fig = plt.figure() + axis = fig.add_subplot(1, 1, 1) + axis.set_yscale('log') + axis.plot(data_disc_fourier0[0], data_disc_fourier0[2], 'b+', + data_disc_fourier1[0], data_disc_fourier1[2], 'bx', + data_disc_fourier20[0], data_disc_fourier20[2], 'b*', + data_disc_real_s[0], data_disc_real_s[2], 'rs', + data_disc_real_d[0], data_disc_real_d[2], 'ro') + axis.set_title('Variation of the resolution (disc)') + axis.set_xlabel('res [nm]') + axis.set_ylabel('duration [s]') + + + data_shelve.close() + + + +if __name__ == "__main__": + try: + compare_method_errors_res() + except: + type, value, tb = sys.exc_info() + traceback.print_exc() + pdb.post_mortem(tb) \ No newline at end of file diff --git a/scripts/compare_methods.py b/scripts/compare_methods.py index 27262a4..7327b7f 100644 --- a/scripts/compare_methods.py +++ b/scripts/compare_methods.py @@ -34,8 +34,8 @@ def compare_methods(): # Create magnetic shape: geometry = 'slab' if geometry == 'slab': - center = (0, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts with 0! - width = (1, dim[1]/2., dim[2]/2.) # in px (z, y, x) + center = (0, dim[1]/2., dim[2]/2.) # in px (z, y, x) index starts with 0! + width = (1, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) mag_shape = mc.Shapes.slab(dim, center, width) phase_ana = an.phase_mag_slab(dim, res, phi, center, width, b_0) elif geometry == 'disc': diff --git a/scripts/create_multiple_samples.py b/scripts/create_multiple_samples.py new file mode 100644 index 0000000..bddaa81 --- /dev/null +++ b/scripts/create_multiple_samples.py @@ -0,0 +1,62 @@ +# -*- coding: utf-8 -*- +"""Create random magnetic distributions.""" + +import random as rnd +import pdb, traceback, sys +import numpy as np +from numpy import pi + +import pyramid.magcreator as mc +from pyramid.magdata import MagData +import pyramid.phasemapper as pm +import pyramid.projector as pj +import pyramid.holoimage as hi +from pyramid.phasemap import PhaseMap + + +def create_multiple_samples(): + '''Calculate, display and save a random magnetic distribution to file. + Arguments: + None + Returns: + None + + ''' + filename = '../output/mag_dist_multiple_samples.txt' + res = 10.0 # nm + dim = (64, 128, 128) + # Slab: + center = (32, 32, 32) # in px (z, y, x), index starts with 0! + width = (48, 48, 48) # in px (z, y, x) + mag_shape_slab = mc.Shapes.slab(dim, center, width) + # Disc: + center = (32, 32, 96) # in px (z, y, x), index starts with 0! + radius = 24 # in px + height = 24 # in px + mag_shape_disc = mc.Shapes.disc(dim, center, radius, height) + # Sphere: + center = (32, 96, 64) # in px (z, y, x), index starts with 0! + radius = 24 # in px + mag_shape_sphere = mc.Shapes.sphere(dim, center, radius) + # Create lists for magnetic objects: + mag_shape_list = [mag_shape_slab, mag_shape_disc, mag_shape_sphere] + phi_list = [pi/4, pi/2, pi] + # Create magnetic distribution: + magnitude = mc.create_mag_dist_comb(mag_shape_list, phi_list) + mag_data = MagData(res, magnitude) + mag_data.quiver_plot('z', dim[0]/2) + mag_data.save_to_llg(filename) + projection = pj.simple_axis_projection(mag_data) + phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab')) + phase_map.display() + hi.display(hi.holo_image(phase_map, 0.5)) + + + +if __name__ == "__main__": + try: + create_multiple_samples() + except: + type, value, tb = sys.exc_info() + traceback.print_exc() + pdb.post_mortem(tb) \ No newline at end of file diff --git a/scripts/reconstruct_random_pixels.py b/scripts/reconstruct_random_pixels.py index a705273..1cf7c70 100644 --- a/scripts/reconstruct_random_pixels.py +++ b/scripts/reconstruct_random_pixels.py @@ -24,11 +24,11 @@ def reconstruct_random_distribution(): ''' # Input parameters: - n_pixel = 10 - dim = (32, 32, 32) + n_pixel = 5 + dim = (1, 16, 16) b_0 = 1 # in T res = 10.0 # in nm - rnd.seed(12) + rnd.seed(18) threshold = 0 # Create lists for magnetic objects: @@ -43,6 +43,7 @@ def reconstruct_random_distribution(): # Create magnetic distribution: magnitude = mc.create_mag_dist_comb(mag_shape_list, phi_list, magnitude_list) mag_data = MagData(res, magnitude) + mag_data.quiver_plot() # Display phase map and holography image: projection = pj.simple_axis_projection(mag_data) phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab', b_0)) -- GitLab