Minor changes in some modules, splitted "compare_method_erros" script

holoimage: densitiy default set to 1

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
 

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()
 
     

scripts/compare_method_errors_res.py

+# -*- 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

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':

scripts/create_multiple_samples.py

+# -*- 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

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))