Updated long overdue scripts to work with new layout and refactoring!

scripts/phasemap/phasemap_tilt_series.py

-# -*- coding: utf-8 -*-
-"""Create PhaseMaps from MagDatas via simple axis projection."""
-
-
-import numpy as np
-import pyramid as py
-import logging.config
-
-
-logging.config.fileConfig(py.LOGGING_CONFIG, disable_existing_loggers=False)
-
-###################################################################################################
-filename = 'magdata_mc_homog_sphere_dim=(32, 32, 32).nc'
-b_0 = 1.
-axis = 'z'
-dim_uv = None
-
-angles = np.linspace()
-###################################################################################################
-
-mag_data = py.VectorData.load_from_netcdf4(filename)
-phase_map = py.pm(mag_data, mode=axis, dim_uv=dim_uv, b_0=b_0)
-phase_map.save_to_netcdf4('phasemap_{}_axis={}'.format(filename.replace('magdata_', ''), axis))
-phase_map.display_combined()
-
-
-dim = mag_data.dim
-dim_uv = (500, 200)
-angles = np.arange(-60, 61, 5)#[0, 20, 40, 60]
-
-#mag_data_xy = mag_data.copy()
-# mag_data_xy.field[2] = 0
-#
-#mag_data_z = mag_data.copy()
-# mag_data_z.field[0] = 0
-# mag_data_z.field[1] = 0
-
-# Iterate over all angles:
-for angle in angles:
-    angle_rad = np.pi/2 + angle*np.pi/180
-    projector = XTiltProjector(dim, angle_rad, dim_uv)
-    mag_proj = projector(mag_data)
-    phase_map = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))(mag_proj)
-    phase_map.display_phase('Phase Map Nanowire Tip', cmap='gray')
-    plt.savefig(PATH+'_nanowire_xtilt_{}.png'.format(angle), dpi=500)
-#    phase_map = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))(mag_proj)
-#    phase_map.display_combined('Phase Map Nanowire Tip', gain=gain,
-#                               interpolation='bilinear')
-#    plt.savefig(PATH+'_nanowire_xtilt_{}.png'.format(angle), dpi=500)
-#    mag_proj.scale_down(2)
-#    axis = mag_proj.quiver_plot()
-#    plt.savefig(PATH+'_nanowire_mag_xtilt_{}.png'.format(angle), dpi=500)
-#    axis = mag_proj.quiver_plot(log=True)
-#    plt.savefig(PATH+'_nanowire_mag_log_xtilt_{}.png'.format(angle), dpi=500)
-    # Close plots:
-    plt.close('all')
-    gc.collect()
-    print 'RSS = {:.2f} MB'.format(proc.memory_info().rss/1024.**2)
-    print angle, 'deg. done!'

scripts/reconstruction/reconstruction_2d_from_phasemap.py

 # -*- coding: utf-8 -*-
 """Reconstruct a magnetization distributions from a single phase map."""
 
+from __future__ import print_function
 
 import numpy as np
 import pyramid as pr
 from jutil.taketime import TakeTime
 import logging.config
 
-
 logging.config.fileConfig(pr.LOGGING_CONFIG, disable_existing_loggers=False)
 
 ###################################################################################################
-phase_name = 'phasemap_gui_Pha1'
+phase_name = 'phasemap_gui_150917_skyrm_M_100mT'
 b_0 = 1  # in T
 lam = 1E-1
 max_iter = 100
@@ -21,7 +21,7 @@ order = 1
 
 
 # Load phase map:
-phase_map = pr.PhaseMap.load_from_netcdf4(phase_name+'.nc')
+phase_map = pr.PhaseMap.load_from_hdf5(phase_name + '.hdf5')
 phase_map.pad((buffer_pixel, buffer_pixel))
 dim = (1,) + phase_map.dim_uv
 
@@ -41,14 +41,17 @@ with TakeTime('reconstruction time'):
 if order >= 1:
     offset, ramp = param_cache[0][0], (param_cache[1][0], param_cache[2][0])
 elif order == 0:
-    offset = param_cache[0][0]
+    offset, ramp = param_cache[0][0], (0, 0)
+else:
+    offset, ramp = 0, (0, 0)
 mag_data_buffer = mag_data_rec.copy()
 mag_data_rec.crop((0, buffer_pixel, buffer_pixel))
 mag_name = '{}_lam={}'.format(phase_name.replace('phasemap', 'magdata_rec'), lam)
-mag_data_rec.save_to_netcdf4(mag_name+'.nc')
+mag_data_rec = mag_data_rec.flip(axis='y')
+mag_data_rec.save_to_hdf5(mag_name + '.hdf5', overwrite=True)
 
 # Plot stuff:
-mag_data_rec.quiver_plot('Reconstructed Distribution', ar_dens=int(np.ceil(np.max(dim)/128.)))
+mag_data_rec.quiver_plot('Reconstructed Distribution', ar_dens=int(np.ceil(np.max(dim) / 128.)))
 phase_map.crop((buffer_pixel, buffer_pixel))
 phase_map.display_combined('Input Phase')
 phase_map -= fwd_model.ramp(index=0)
@@ -56,13 +59,13 @@ phase_map.display_combined('Input Phase (ramp corrected)')
 phase_map_rec = pr.pm(mag_data_rec)
 title = 'Reconstructed Phase'
 if order >= 0:
-    print 'offset:', offset
+    print('offset:', offset)
     title += ', fitted Offset: {:.2g} [rad]'.format(offset)
 if order >= 1:
-    print 'ramp:', ramp
+    print('ramp:', ramp)
     title += ', (Fitted Ramp: (u:{:.2g}, v:{:.2g}) [rad/nm]'.format(*ramp)
 phase_map_rec.display_combined(title)
-difference = (phase_map_rec.phase-phase_map.phase).mean()
-(phase_map_rec-phase_map).display_phase('Difference (mean: {:.2g})'.format(difference))
+difference = (phase_map_rec.phase - phase_map.phase).mean()
+(phase_map_rec - phase_map).display_phase('Difference (mean: {:.2g})'.format(difference))
 if order is not None:
     fwd_model.ramp(0).display_combined('Fitted Ramp')

scripts/reconstruction/reconstruction_3d_from_magdata.py

@@ -12,7 +12,7 @@ from jutil.taketime import TakeTime
 logging.config.fileConfig(pr.LOGGING_CONFIG, disable_existing_loggers=False)
 
 ###################################################################################################
-mag_name = 'magdata_mc_vortex_disc'
+mag_name = 'magdata_mc_horseshoe'
 dim_uv = None
 angles = np.linspace(-90, 90, num=7)
 axes = {'x': True, 'y': True}
@@ -25,73 +25,67 @@ offset_max = 2
 ramp_max = 0.01
 order = 1
 show_input = True
-# TODO: toggle multiprocessing if nprocs > 1
+nprocs = mp.cpu_count()  # or 1 for non-multiprocessing
 ###################################################################################################
 
 
 if __name__ == '__main__':
-
     mp.freeze_support()
-
     # Load magnetization distribution:
     mag_data = pr.VectorData.load_from_hdf5(mag_name + '.hdf5')
     dim = mag_data.dim
-
     # Construct data set and regularisator:
     data = pr.DataSet(mag_data.a, mag_data.dim, b_0)
-
     # Construct projectors:
     projectors = []
     for angle in angles:
-        angle_rad = angle*np.pi/180
+        angle_rad = angle * np.pi / 180
         if axes['x']:
             projectors.append(pr.XTiltProjector(mag_data.dim, angle_rad, dim_uv))
         if axes['y']:
             projectors.append(pr.YTiltProjector(mag_data.dim, angle_rad, dim_uv))
-
     # Add projectors and construct according phase maps:
     data.projectors = projectors
     data.phase_maps = data.create_phase_maps(mag_data)
-
     for i, phase_map in enumerate(data.phase_maps):
         offset = np.random.uniform(-offset_max, offset_max)
         ramp_u = np.random.uniform(-ramp_max, ramp_max)
         ramp_v = np.random.uniform(-ramp_max, ramp_max)
         phase_map += pr.Ramp.create_ramp(phase_map.a, phase_map.dim_uv, (offset, ramp_u, ramp_v))
-
     # Add noise if necessary:
     if noise != 0:
         for i, phase_map in enumerate(data.phase_maps):
             phase_map.phase += np.random.normal(0, noise, phase_map.dim_uv)
             data.phase_maps[i] = phase_map
-
     # Plot input:
     if show_input:
         for i, phase_map in enumerate(data.phase_maps):
             phase_map.display_phase()
-
     # Construct mask:
-#    if use_internal_mask:
-#        data.mask = mag_data.get_mask()  # Use perfect mask from mag_data!
-#    else:
-#        data.set_3d_mask()  # Construct mask from 2D phase masks!
-
+    if use_internal_mask:
+        data.mask = mag_data.get_mask()  # Use perfect mask from mag_data!
+    else:
+        data.set_3d_mask()  # Construct mask from 2D phase masks!
     # Construct regularisator, forward model and costfunction:
-    fwd_model = pr.DistributedForwardModel(data, ramp_order=order, nprocs=4)
+    if nprocs > 1:
+        fwd_model = pr.DistributedForwardModel(data, ramp_order=order, nprocs=nprocs)
+    else:
+        fwd_model = pr.ForwardModel(data, ramp_order=order)
     reg = pr.FirstOrderRegularisator(data.mask, lam, add_params=fwd_model.ramp.n)
     cost = pr.Costfunction(fwd_model, reg)
-
     # Reconstruct and save:
     with TakeTime('reconstruction time'):
         mag_data_rec = pr.reconstruction.optimize_linear(cost, max_iter=max_iter)
         param_cache = cost.fwd_model.ramp.param_cache
     fwd_model.finalize()
     mag_name_rec = mag_name.replace('magdata', 'magdata_rec')
-    mag_data_rec.save_to_netcdf4(mag_name_rec+'.nc')
-
+    mag_data_rec.save_to_hdf5(mag_name_rec + '.hdf5', overwrite=True)
     # Plot stuff:
-    data.display_mask(ar_dens=np.ceil(np.max(dim)/64.))
-    mag_data.quiver_plot3d('Original Distribution', ar_dens=np.ceil(np.max(dim)/64.))
-    mag_data_rec.quiver_plot3d('Reconstructed Distribution', ar_dens=np.ceil(np.max(dim)/64.))
-    mag_data_rec.quiver_plot3d('Reconstructed Distribution', ar_dens=np.ceil(np.max(dim)/64.),
+    data.display_mask(ar_dens=np.ceil(np.max(dim) / 64.))
+    mag_data.quiver_plot3d('Original Distribution', ar_dens=np.ceil(np.max(dim) / 64.))
+    mag_data_rec.quiver_plot3d('Reconstructed Distribution', ar_dens=np.ceil(np.max(dim) / 64.))
+    mag_data_rec.quiver_plot3d('Reconstructed Distribution', ar_dens=np.ceil(np.max(dim) / 64.),
                                coloring='amplitude')
+    import matplotlib.pyplot as plt
+
+    plt.show()