From 6461da37c1da1b0ecff2271b6164d1cedc5876a4 Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Tue, 30 Dec 2014 16:58:12 +0100
Subject: [PATCH] Added diagnostics class!
---
pyramid/__init__.py | 4 +
pyramid/costfunction.py | 19 +++
pyramid/diagnostics.py | 96 ++++++++++++++
pyramid/fft.py | 14 +-
pyramid/forwardmodel.py | 2 +-
pyramid/magdata.py | 4 +-
pyramid/phasemap.py | 4 +-
pyramid/reconstruction.py | 6 +-
pyramid/regularisator.py | 4 -
...rick_nanowire_reconstruction_simulation.py | 121 ++++++++++++++++++
.../vtk_conversion_nanowire_full.py | 3 +-
.../vtk_conversion_nanowire_full_90deg.py | 83 ++++++------
.../vtk_conversion_nanowire_segment.py | 2 +-
.../create_multiple_samples.py | 2 +-
.../reconstruction_linear_test.py | 32 ++++-
15 files changed, 329 insertions(+), 67 deletions(-)
create mode 100644 pyramid/diagnostics.py
create mode 100644 scripts/collaborations/patrick_nanowire_reconstruction_simulation.py
diff --git a/pyramid/__init__.py b/pyramid/__init__.py
index ad35b5a..ad7b174 100644
--- a/pyramid/__init__.py
+++ b/pyramid/__init__.py
@@ -28,6 +28,10 @@ costfunction
Class for the evaluation of the cost of a function.
reconstruction
Reconstruct magnetic distributions from given phasemaps.
+regularisator
+ Class to instantiate different regularisation strategies.
+fft
+ Class for custom FFT functions using numpy or FFTW.
Subpackages
-----------
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index b177e05..50c41ff 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -132,8 +132,27 @@ class Costfunction(object):
def hess_diag(self, _):
# TODO: Docstring!
+ # TODO: What for again?
return np.ones(self.n)
+ def estimate_lambda(self):
+ # TODO: Docstring!
+ # TODO: Not very efficient? Is this even correct?
+
+ def unit_vec(length, index):
+ result = np.zeros(length)
+ result[index] = 1
+ return result
+
+ fwd, reg = self.fwd_model, self.regularisator
+ trace_fwd = np.sum([fwd.jac_T_dot(None, fwd.jac_dot(None, unit_vec(self.n, i)))
+ for i in range(self.n)])
+ trace_reg = np.sum([reg(unit_vec(self.n, i)) for i in range(self.n)])
+ print 'fwd:', trace_fwd
+ print 'reg:', trace_reg
+ import pdb; pdb.set_trace()
+ return trace_fwd / trace_reg
+
class CFAdapterScipyCG(LinearOperator):
diff --git a/pyramid/diagnostics.py b/pyramid/diagnostics.py
new file mode 100644
index 0000000..fbec1ff
--- /dev/null
+++ b/pyramid/diagnostics.py
@@ -0,0 +1,96 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Dec 11 17:20:27 2014
+
+@author: Jan
+"""
+
+
+import numpy as np
+
+import jutil
+
+from pyramid import fft
+from pyramid.magdata import MagData
+from pyramid.phasemap import PhaseMap
+
+class Diagnostics(object):
+
+ # TODO: Docstrings and position of properties!
+
+ def __init__(self, x_rec, cost, max_iter=100):
+ self.x_rec = x_rec
+ self.cost = cost
+ self.max_iter = max_iter
+ self.fwd_model = cost.fwd_model
+ self.Se_inv = self.cost.Se_inv
+ self.dim = self.cost.data_set.dim
+ self.row_idx = None
+ self.set_position(0)#(0, self.dim[0]//2, self.dim[1]//2, self.dim[2]//2))
+ self._A = jutil.operator.CostFunctionOperator(self.cost, self.x_rec)
+ self._P = jutil.preconditioner.CostFunctionPreconditioner(self.cost, self.x_rec)
+
+ def set_position(self, pos):
+ # TODO: Does not know about the mask... thus gives wrong results or errors
+# m, z, y, x = pos
+# row_idx = m*np.prod(self.dim) + z*self.dim[1]*self.dim[2] + y*self.dim[2] + x
+ row_idx = pos
+ if row_idx != self.row_idx:
+ self.row_idx = row_idx
+ self._updated_std = False
+ self._updated_gain_row = False
+ self._updated_avrg_kern_row = False
+ self._updated_measure_contribution = False
+
+ @property
+ def std(self):
+ if not self._updated_std:
+ e_i = fft.zeros(self.cost.n, dtype=fft.FLOAT)
+ e_i[self.row_idx] = 1
+ row = jutil.cg.conj_grad_solve(self._A, e_i, P=self._P, max_iter=self.max_iter)
+ self._m_inv_row = row
+ self._std = np.sqrt(self._m_inv_row[self.row_idx])
+ self._updated_std = True
+ return self._std
+
+ @property
+ def gain_row(self):
+ if not self._updated_gain_row:
+ self.std # evoke to update self._m_inv_row if necessary # TODO: make _m_inv_row checked!
+ self._gain_row = self.Se_inv.dot(self.fwd_model.jac_dot(self.x_rec, self._m_inv_row))
+ self._updated_gain_row = True
+ return self._gain_row
+
+ @property
+ def avrg_kern_row(self):
+ if not self._updated_avrg_kern_row:
+ self._avrg_kern_row = self.fwd_model.jac_T_dot(self.x_rec, self.gain_row)
+ self._updated_avrg_kern_row = True
+ return self._avrg_kern_row
+
+ @property
+ def measure_contribution(self):
+ if not self._updated_measure_contribution:
+ cache = self.fwd_model.jac_dot(self.x_rec, fft.ones(self.cost.n, fft.FLOAT))
+ cache = self.fwd_model.jac_T_dot(self.x_rec, self.Se_inv.dot(cache))
+ mc = jutil.cg.conj_grad_solve(self._A, cache, P=self._P, max_iter=self.max_iter)
+ self._measure_contribution = mc
+ self._updated_measure_contribution = True
+ return self._measure_contribution
+
+ def get_avg_kernel(self, pos=None):
+ if pos is not None:
+ self.set_position(pos)
+ mag_data_avg_kern = MagData(self.cost.data_set.a, fft.zeros((3,)+self.dim))
+ mag_data_avg_kern.set_vector(self.avrg_kern_row, mask=self.cost.data_set.mask)
+ return mag_data_avg_kern
+
+ def get_gain_maps(self, pos=None):
+ if pos is not None:
+ self.set_position(pos)
+ hp = self.cost.data_set.hook_points
+ result = []
+ for i, projector in enumerate(self.cost.data_set.projectors):
+ gain = self.gain_row[hp[i]:hp[i+1]].reshape(projector.dim_uv)
+ result.append(PhaseMap(self.cost.data_set.a, gain))
+ return result
diff --git a/pyramid/fft.py b/pyramid/fft.py
index e595437..1c7733d 100644
--- a/pyramid/fft.py
+++ b/pyramid/fft.py
@@ -171,7 +171,15 @@ def load_wisdom(fname):
# Array setups:
-def zeros(shape, dtype):
+def empty(shape, dtype=FLOAT):
+ # TODO: Docstring!
+ result = np.empty(shape, dtype)
+ if pyfftw is not None:
+ result = pyfftw.n_byte_align(result, pyfftw.simd_alignment)
+ return result
+
+
+def zeros(shape, dtype=FLOAT):
# TODO: Docstring!
result = np.zeros(shape, dtype)
if pyfftw is not None:
@@ -179,9 +187,9 @@ def zeros(shape, dtype):
return result
-def empty(shape, dtype):
+def ones(shape, dtype=FLOAT):
# TODO: Docstring!
- result = np.empty(shape, dtype)
+ result = np.ones(shape, dtype)
if pyfftw is not None:
result = pyfftw.n_byte_align(result, pyfftw.simd_alignment)
return result
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index eacc079..f2758b9 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -93,7 +93,7 @@ class ForwardModel(object):
`vector`.
'''
- self.mag_data.magnitude[:] = 0
+ self.mag_data.magnitude[...] = 0
self.mag_data.set_vector(vector, self.data_set.mask)
result = np.zeros(self.m)
hp = self.hook_points
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index 9b64f6d..7fe3e6c 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -498,9 +498,9 @@ class MagData(object):
# Setup quiver:
dim_uv = plot_u.shape
ad = ar_dens
- xx, yy = np.meshgrid(np.arange(dim_uv[0]), np.arange(dim_uv[1]))
+ xx, yy = np.meshgrid(np.arange(dim_uv[1]), np.arange(dim_uv[0]))
axis.quiver(xx[::ad, ::ad], yy[::ad, ::ad], plot_u[::ad, ::ad], plot_v[::ad, ::ad],
- pivot='middle', units='xy', angles=angles[::ad, ::ad],
+ pivot='middle', units='xy', angles=angles[::ad, ::ad], minlength=0.25,
scale_units='xy', scale=1./ad, headwidth=6, headlength=7)
axis.set_xlim(-1, dim_uv[1])
axis.set_ylim(-1, dim_uv[0])
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index 4b962dd..d60b410 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -344,8 +344,8 @@ class PhaseMap(object):
axis.set_ylim(0, self.dim_uv[0])
axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
- axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:.3g}'.format(x*self.a)))
- axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:.3g}'.format(x*self.a)))
+ axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*self.a)))
+ axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*self.a)))
axis.tick_params(axis='both', which='major', labelsize=14)
axis.set_title(title, fontsize=18)
axis.set_xlabel('u-axis [nm]', fontsize=15)
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index 4da85c7..d9386ee 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -83,7 +83,7 @@ class PrintIterator(object):
return self.iteration
-def optimize_linear(data, regularisator=None, max_iter=None, info=None):
+def optimize_linear(data, regularisator=None, max_iter=None):
'''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
Blazingly fast for l2-based cost functions.
@@ -114,9 +114,7 @@ def optimize_linear(data, regularisator=None, max_iter=None, info=None):
# Create and return fitting MagData object:
mag_opt = MagData(data.a, np.zeros((3,) + data.dim))
mag_opt.set_vector(x_opt, data.mask)
- if info is not None:
- info[:] = cost.chisq, cost.chisq_m, cost.chisq_a
- return mag_opt
+ return mag_opt, cost
def optimize_nonlin(data, first_guess=None, regularisator=None):
diff --git a/pyramid/regularisator.py b/pyramid/regularisator.py
index 03cdaaa..5c115bb 100644
--- a/pyramid/regularisator.py
+++ b/pyramid/regularisator.py
@@ -52,12 +52,10 @@ class Regularisator(object):
def jac(self, x):
# TODO: Docstring!
- self._log.debug('Calling jac')
return self.lam * self.norm.jac(x)
def hess_dot(self, x, vector):
# TODO: Docstring!
-# self._log.debug('Calling hess_dot') # TODO: Profiler says this was slow...
return self.lam * self.norm.hess_dot(x, vector)
def hess_diag(self, x, vector):
@@ -85,12 +83,10 @@ class NoneRegularisator(Regularisator):
def jac(self, x):
# TODO: Docstring!
- self._log.debug('Calling jac')
return np.zeros_like(x)
def hess_dot(self, x, vector):
# TODO: Docstring!
- self._log.debug('Calling hess_dot')
return np.zeros_like(vector)
def hess_diag(self, x, vector):
diff --git a/scripts/collaborations/patrick_nanowire_reconstruction_simulation.py b/scripts/collaborations/patrick_nanowire_reconstruction_simulation.py
new file mode 100644
index 0000000..3c1edba
--- /dev/null
+++ b/scripts/collaborations/patrick_nanowire_reconstruction_simulation.py
@@ -0,0 +1,121 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Nov 14 14:17:56 2014
+
+@author: Jan
+"""
+
+from PIL import Image
+import numpy as np
+
+from jutil.taketime import TakeTime
+from pyramid import * # analysis:ignore
+
+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)
+
+###################################################################################################
+gain = 5
+b_0 = 1
+lam = 1E-3
+length = 300
+cutoff = 50
+trans = 40
+pads = [100, 250, 500, 800]
+INPATH = '../../output/vtk data\longtube_withcap/'
+OUTPATH = '../../output/patrick/'
+FILE = 'CoFeB_tube_cap_4nm_lying_down'
+###################################################################################################
+
+# Read in big 3D files:
+mag_data_3d = MagData.load_from_netcdf4(INPATH+FILE+'.nc')
+a = mag_data_3d.a
+dim_3d = mag_data_3d.dim
+
+# Make 2D projection:
+dim_uv_big = (dim_3d[1]+200, dim_3d[2]+200)
+projector = SimpleProjector(dim_3d, dim_uv=dim_uv_big)
+mag_data_big = projector(mag_data_3d)
+phase_map_big = pm(mag_data_big, dim_uv=dim_uv_big)
+phase_map_big.display_combined(gain=gain)
+mask_big = mag_data_big.get_mask()
+dim_big = mag_data_big.dim
+dim_uv_big = phase_map_big.dim_uv
+
+# Reconstruction of complete wire:
+data_set = DataSet(a, dim_big, b_0, mask_big)
+data_set.append(phase_map_big, SimpleProjector(dim_big))
+regularisator = FirstOrderRegularisator(mask_big, lam, p=2)
+with TakeTime('reconstruction time'):
+ mag_data_big_rec = rc.optimize_linear(data_set, regularisator=regularisator, max_iter=50)[0]
+phase_map_big_rec = pm(mag_data_big_rec)
+axis_p, axis_h = phase_map_big.display_combined('Reconstruction (complete wire)', gain=gain)
+axis_p.axhline(y=length, linewidth=1, color='b', linestyle='-')
+axis_h.axhline(y=length, linewidth=1, color='b', linestyle='-')
+plt.savefig(OUTPATH+FILE+'RECON_COMPLETE.png'.format(pad), dpi=500)
+plt.close('all')
+
+for pad in pads:
+ # Make smaller cutout with optional padding:
+ mask = mask_big[:, :length+pad, :].copy()
+ mask[:, mask.shape[1]-cutoff:, :] = False
+ dim = mask.shape
+ dim_uv = (dim[1], dim[2])
+ phase = np.zeros(dim_uv)
+ phase[:length, :] = phase_map_big.phase[:length, :]
+ phase_map = PhaseMap(a, phase)
+ phase_map_orig = PhaseMap(a, phase[:length, :])
+
+ # Create DataSet and regularisator:
+ data_set = DataSet(a, dim, b_0, mask)
+ data_set.append(phase_map, SimpleProjector(dim))
+
+ # Create Se_inv
+ mask_Se = np.ones(dim_uv)
+ mask_Se[length:length+pad, :] = 0
+ ramp = np.tile(np.linspace(1, 0, trans), mask_Se.shape[1]).reshape((mask_Se.shape[1], trans)).T
+ mask_Se[length-trans:length, :] = ramp
+ data_set.set_Se_inv_diag_with_masks([mask_Se])
+
+ # Create regularisator:
+ regularisator = FirstOrderRegularisator(mask, lam, p=2)
+
+ # Reconstruct:
+ with TakeTime('reconstruction time'):
+ mag_data_rec = rc.optimize_linear(data_set, regularisator=regularisator, max_iter=50)[0]
+
+ # Calculate additional PhaseMaps:
+ phase_map_rec = pm(mag_data_rec)
+ phase_map_cut = PhaseMap(a, phase_map_rec.phase[:length, :])
+ phase_map_diff = phase_map_orig - phase_map_cut
+
+ # Plotting:
+ phase_map_orig.display_combined('Original distribution', gain=gain)
+ plt.savefig(OUTPATH+FILE+'_{}pxpadding_ORIG.png'.format(pad), dpi=500)
+ axis_p, axis_h = phase_map.display_combined('Original Distribution (padded)', gain=gain)
+ axis_p.axhline(y=length, linewidth=1, color='b', linestyle='-')
+ axis_h.axhline(y=length, linewidth=1, color='b', linestyle='-')
+ axis_p.axhline(y=length-trans, linewidth=1, color='r', linestyle='--')
+ axis_h.axhline(y=length-trans, linewidth=1, color='r', linestyle='--')
+ plt.savefig(OUTPATH+FILE+'_{}ramp_{}pxpadding_ORIG_PAD.png'.format(trans, pad), dpi=500)
+ axis_p, axis_h = phase_map_rec.display_combined('Reconstr. Distribution (padded)', gain=gain)
+ axis_p.axhline(y=length, linewidth=1, color='b', linestyle='-')
+ axis_h.axhline(y=length, linewidth=1, color='b', linestyle='-')
+ axis_p.axhline(y=length-trans, linewidth=1, color='r', linestyle='--')
+ axis_h.axhline(y=length-trans, linewidth=1, color='r', linestyle='--')
+ plt.savefig(OUTPATH+FILE+'_{}ramp_{}pxpadding_RECON_PAD.png'.format(trans, pad), dpi=500)
+ phase_map_cut.display_combined('Reconstr. Distribution', gain=gain)
+ plt.savefig(OUTPATH+FILE+'_{}ramp_{}pxpadding_RECON.png'.format(trans, pad))
+ phase_map_diff.display_combined('Difference')
+ plt.savefig(OUTPATH+FILE+'_{}ramp_{}pxpadding_DIFF.png'.format(trans, pad))
+ axis = mag_data_rec.quiver_plot(ar_dens=4, log=True)
+ axis.axhline(y=length, linewidth=1, color='b', linestyle='-')
+ axis.axhline(y=length-trans, linewidth=1, color='r', linestyle='--')
+ plt.savefig(OUTPATH+FILE+'_{}ramp_{}pxpadding_MAG.png'.format(trans, pad), dpi=500)
+ plt.close('all')
diff --git a/scripts/collaborations/vtk_conversion_nanowire_full.py b/scripts/collaborations/vtk_conversion_nanowire_full.py
index 67f32b3..a6f3edb 100644
--- a/scripts/collaborations/vtk_conversion_nanowire_full.py
+++ b/scripts/collaborations/vtk_conversion_nanowire_full.py
@@ -170,6 +170,7 @@ for i in range(mag_data.magnitude.shape[1]):
if i == 100:
mag_data.quiver_plot(ax_slice=i)
plt.savefig(PATH+'_quiver_500_post.png')
+mag_data.save_to_netcdf4(PATH+'_90deg_around_z.nc')
# Iterate over all angles:
for angle in angles:
angle_rad = angle * np.pi/180
@@ -199,7 +200,7 @@ for angle in angles:
PM = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))
phase_map_tip = PhaseMap(mag_data.a, PM(projector(mag_data_bot)).phase[50+shift:225+shift, :])
phase_map_tip.display_combined('Phase Map Nanowire Bottom', gain=gain,
- interpolation='bilinear')
+ b interpolation='bilinear')
plt.savefig(PATH+'_nanowire_bot_xtilt_{}_no_frame.png'.format(angle))
mag_data_bot.scale_down()
mag_proj_bot = projector_scaled(mag_data_bot)
diff --git a/scripts/collaborations/vtk_conversion_nanowire_full_90deg.py b/scripts/collaborations/vtk_conversion_nanowire_full_90deg.py
index acf17b6..d1ab275 100644
--- a/scripts/collaborations/vtk_conversion_nanowire_full_90deg.py
+++ b/scripts/collaborations/vtk_conversion_nanowire_full_90deg.py
@@ -136,45 +136,46 @@ for i in range(mag_data.magnitude.shape[2]):
magnitude_new[1, ..., i] = -z_rot
magnitude_new[2, ..., i] = y_rot
mag_data.magnitude = magnitude_new
+mag_data.save_to_netcdf4(PATH+'_lying_down.nc')
# Iterate over all angles:
-for angle in angles:
- angle_rad = angle * np.pi/180
- projector = YTiltProjector(dim, angle_rad, dim_uv)
- projector_scaled = YTiltProjector((dim[0]/2, dim[1]/2, dim[2]/2), angle_rad,
- (dim_uv[0]/2, dim_uv[1]/2))
- # Tip:
- mag_data_tip = MagData(mag_data.a, mag_data.magnitude[:, :, 608:, :])
- PM = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))
- phase_map_tip = PhaseMap(mag_data.a, PM(projector(mag_data_tip)).phase[350:530, :])
- phase_map_tip.display_combined('Phase Map Nanowire Tip', gain=gain,
- interpolation='bilinear')
- plt.savefig(PATH+'_nanowire_tip_xtilt_{}.png'.format(angle))
- mag_data_tip.scale_down()
- mag_proj_tip = projector_scaled(mag_data_tip)
- axis = mag_proj_tip.quiver_plot()
-# axis.set_xlim(17, 55)
-# axis.set_ylim(180-shift/2, 240-shift/2)
- plt.savefig(PATH+'_nanowire_tip_mag_xtilt_{}.png'.format(angle))
- axis = mag_proj_tip.quiver_plot(log=True)
-# axis.set_xlim(17, 55)
-# axis.set_ylim(180-shift/2, 240-shift/2)
- plt.savefig(PATH+'_nanowire_tip_mag_log_xtilt_{}.png'.format(angle))
- # Bottom:
- mag_data_bot = MagData(mag_data.a, mag_data.magnitude[:, :, :300, :])
- PM = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))
- phase_map_tip = PhaseMap(mag_data.a, PM(projector(mag_data_bot)).phase[50:225, :])
- phase_map_tip.display_combined('Phase Map Nanowire Bottom', gain=gain,
- interpolation='bilinear')
- plt.savefig(PATH+'_nanowire_bot_xtilt_{}_no_frame.png'.format(angle))
- mag_data_bot.scale_down()
- mag_proj_bot = projector_scaled(mag_data_bot)
- axis = mag_proj_bot.quiver_plot()
-# axis.set_xlim(17, 55)
-# axis.set_ylim(50+shift/2, 110+shift/2)
- plt.savefig(PATH+'_nanowire_bot_mag_xtilt_{}.png'.format(angle))
- axis = mag_proj_bot.quiver_plot(log=True)
-# axis.set_xlim(17, 55)
-# axis.set_ylim(50+shift/2, 110+shift/2)
- plt.savefig(PATH+'_nanowire_bot_mag_log_xtilt_{}.png'.format(angle))
- # Close plots:
- plt.close('all')
+#for angle in angles:
+# angle_rad = angle * np.pi/180
+# projector = YTiltProjector(dim, angle_rad, dim_uv)
+# projector_scaled = YTiltProjector((dim[0]/2, dim[1]/2, dim[2]/2), angle_rad,
+# (dim_uv[0]/2, dim_uv[1]/2))
+# # Tip:
+# mag_data_tip = MagData(mag_data.a, mag_data.magnitude[:, :, 608:, :])
+# PM = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))
+# phase_map_tip = PhaseMap(mag_data.a, PM(projector(mag_data_tip)).phase[350:530, :])
+# phase_map_tip.display_combined('Phase Map Nanowire Tip', gain=gain,
+# interpolation='bilinear')
+# plt.savefig(PATH+'_nanowire_tip_xtilt_{}.png'.format(angle))
+# mag_data_tip.scale_down()
+# mag_proj_tip = projector_scaled(mag_data_tip)
+# axis = mag_proj_tip.quiver_plot()
+## axis.set_xlim(17, 55)
+## axis.set_ylim(180-shift/2, 240-shift/2)
+# plt.savefig(PATH+'_nanowire_tip_mag_xtilt_{}.png'.format(angle))
+# axis = mag_proj_tip.quiver_plot(log=True)
+## axis.set_xlim(17, 55)
+## axis.set_ylim(180-shift/2, 240-shift/2)
+# plt.savefig(PATH+'_nanowire_tip_mag_log_xtilt_{}.png'.format(angle))
+# # Bottom:
+# mag_data_bot = MagData(mag_data.a, mag_data.magnitude[:, :, :300, :])
+# PM = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))
+# phase_map_tip = PhaseMap(mag_data.a, PM(projector(mag_data_bot)).phase[50:225, :])
+# phase_map_tip.display_combined('Phase Map Nanowire Bottom', gain=gain,
+# interpolation='bilinear')
+# plt.savefig(PATH+'_nanowire_bot_xtilt_{}_no_frame.png'.format(angle))
+# mag_data_bot.scale_down()
+# mag_proj_bot = projector_scaled(mag_data_bot)
+# axis = mag_proj_bot.quiver_plot()
+## axis.set_xlim(17, 55)
+## axis.set_ylim(50+shift/2, 110+shift/2)
+# plt.savefig(PATH+'_nanowire_bot_mag_xtilt_{}.png'.format(angle))
+# axis = mag_proj_bot.quiver_plot(log=True)
+## axis.set_xlim(17, 55)
+## axis.set_ylim(50+shift/2, 110+shift/2)
+# plt.savefig(PATH+'_nanowire_bot_mag_log_xtilt_{}.png'.format(angle))
+# # Close plots:
+# plt.close('all')
diff --git a/scripts/collaborations/vtk_conversion_nanowire_segment.py b/scripts/collaborations/vtk_conversion_nanowire_segment.py
index 1f00810..9fac2fb 100644
--- a/scripts/collaborations/vtk_conversion_nanowire_segment.py
+++ b/scripts/collaborations/vtk_conversion_nanowire_segment.py
@@ -142,7 +142,7 @@ else:
mag_data.quiver_plot()
###################################################################################################
# Phasemapping:
-phase_map = pm(mag_data)
+phase_map = pm(mag_data, dim_uv=(300, 300))
(-phase_map).display_combined(title=r'Combined Plot (B$_0$={} T, Cos x {})'.format(b_0, gain),
gain=gain)
plt.savefig(PATH+'_{}T_cosx{}.png'.format(b_0, gain))
diff --git a/scripts/create distributions/create_multiple_samples.py b/scripts/create distributions/create_multiple_samples.py
index 1c09a0e..e978e51 100644
--- a/scripts/create distributions/create_multiple_samples.py
+++ b/scripts/create distributions/create_multiple_samples.py
@@ -26,7 +26,7 @@ if not os.path.exists(directory):
# Input parameters:
filename = directory + '/mag_dist_multiple_samples.txt'
a = 10.0 # nm
-dim = (64, 128, 128)
+dim = (64, 128, 256)
# Slab:
center = (32, 32, 32) # in px (z, y, x), index starts with 0!
width = (48, 48, 48) # in px (z, y, x)
diff --git a/scripts/reconstruction/reconstruction_linear_test.py b/scripts/reconstruction/reconstruction_linear_test.py
index b26c728..7da8b7a 100644
--- a/scripts/reconstruction/reconstruction_linear_test.py
+++ b/scripts/reconstruction/reconstruction_linear_test.py
@@ -52,9 +52,9 @@ shape = mc.Shapes.disc(dim, center, radius_shell, height)
magnitude = mc.create_mag_dist_vortex(shape)
mag_data = MagData(a, magnitude)
-#mag_data.quiver_plot('z-projection', proj_axis='z')
-#mag_data.quiver_plot('y-projection', proj_axis='y')
-#mag_data.quiver_plot3d('Original distribution')
+mag_data.quiver_plot('z-projection', proj_axis='z')
+mag_data.quiver_plot('y-projection', proj_axis='y')
+mag_data.quiver_plot3d('Original distribution')
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)
@@ -90,15 +90,33 @@ if noise != 0:
print('--Reconstruction')
reg = FirstOrderRegularisator(mask, lam, p=2)
-info = []
with TakeTime('reconstruction'):
- mag_data_opt = rc.optimize_linear(data, regularisator=reg, max_iter=100, info=info)
+ mag_data_opt, cost = rc.optimize_linear(data, regularisator=reg, max_iter=100)
+print 'Cost:', cost.chisq
###################################################################################################
print('--Plot stuff')
mag_data_opt.quiver_plot3d('Reconstructed distribution')
-#(mag_data_opt - mag_data).quiver_plot3d('Difference')
-#phase_maps_opt = data.create_phase_maps(mag_data_opt)
+(mag_data_opt - mag_data).quiver_plot3d('Difference')
+phase_maps_opt = data.create_phase_maps(mag_data_opt)
# TODO: iterations in jutil is one to many!
+
+from pyramid.diagnostics import Diagnostics
+
+diag = Diagnostics(mag_data_opt.mag_vec, cost)
+
+print 'std:', diag.std
+gain_maps = diag.get_gain_maps()
+gain_maps[count//2].display_phase()
+diag.get_avg_kernel().quiver_plot3d()
+measure_contribution = diag.measure_contribution
+
+diag.set_position(cost.data_set.mask.sum()//2)
+
+print 'std:', diag.std
+gain_maps = diag.get_gain_maps()
+gain_maps[count//2].display_phase()
+diag.get_avg_kernel().quiver_plot3d()
+measure_contribution = diag.measure_contribution
--
GitLab