Added diagnostics class!

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

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

pyramid/diagnostics.py

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

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

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

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

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)

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

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

scripts/collaborations/patrick_nanowire_reconstruction_simulation.py

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

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)

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

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

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)

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