fixed some minor compilation bugs and updated to current jutil. prepared jutil...

fixed some minor compilation bugs and updated to current jutil. prepared jutil reconstruction routines.

pyramid/costfunction.py

@@ -53,6 +53,8 @@ class Costfunction(object):
         # Extract important information:
         self.y = data_set.phase_vec
         self.Se_inv = data_set.Se_inv
+        self.n = data_set.n * 3
+        self.m = len(self.y)
         self.LOG.debug('Created '+str(self))
 
     def __repr__(self):
@@ -65,10 +67,14 @@ class Costfunction(object):
         return 'Costfunction(data_set=%s, fwd_model=%s, regularisator=%s)' % \
             (self.data_set, self.fwd_model, self.regularisator)
 
+    def init(self, x):
+        self(x)
+
     def __call__(self, x):
         self.LOG.debug('Calling __call__')
-        delta_y = self.fwd_model(x)-self.y
-        return delta_y.dot(self.Se_inv.dot(delta_y)) + self.regularisator(x)
+        delta_y = self.fwd_model(x) - self.y
+        self.chisq = delta_y.dot(self.Se_inv.dot(delta_y)) + self.regularisator(x)
+        return self.chisq
 
     def jac(self, x):
         '''Calculate the derivative of the costfunction for a given magnetization distribution.
@@ -85,7 +91,8 @@ class Costfunction(object):
 
         '''
         self.LOG.debug('Calling jac')
-        return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model(x)-self.y))
+        assert len(x) == self.n
+        return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model(x) - self.y))
                 + self.regularisator.jac(x))
 
     def hess_dot(self, x, vector):
@@ -111,6 +118,9 @@ class Costfunction(object):
         return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model.jac_dot(x, vector)))
                 + self.regularisator.hess_dot(x, vector))
 
+    def hess_diag(self, _):
+        return np.ones(self.n)
+
 
 class CFAdapterScipyCG(LinearOperator):
 

pyramid/reconstruction.py

@@ -83,7 +83,7 @@ class PrintIterator(object):
         return self.iteration
 
 
-def optimize_sparse_cg(data, regularisator=None, maxiter=1000, verbosity=0):
+def optimize_linear(data, regularisator=None, maxiter=1000, verbosity=0):
     '''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
     conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
 
@@ -113,22 +113,21 @@ def optimize_sparse_cg(data, regularisator=None, maxiter=1000, verbosity=0):
         The reconstructed magnetic distribution as a :class:`~.MagData` object.
 
     '''
+    import jutil
     LOG.debug('Calling optimize_sparse_cg')
     # Set up necessary objects:
-    y = data.phase_vec
     fwd_model = ForwardModel(data)
     cost = Costfunction(data, regularisator)
-    # Optimize:
-    A = CFAdapterScipyCG(cost)
-    b = fwd_model.jac_T_dot(None, y)
-    x_opt, info = cg(A, b, callback=PrintIterator(cost, verbosity), maxiter=maxiter)
+    print cost(np.zeros(cost.n))
+    x_opt = jutil.cg.conj_grad_minimize(cost)
+    print cost(x_opt)
     # Create and return fitting MagData object:
     mag_opt = MagData(data.a, np.zeros((3,)+data.dim))
     mag_opt.mag_vec = x_opt
     return mag_opt
 
 
-def optimize_cg(data, first_guess):
+def optimize_nonlin(data, first_guess=None, regularisator=None):
     '''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
     conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
 
@@ -149,19 +148,18 @@ def optimize_cg(data, first_guess):
         The reconstructed magnetic distribution as a :class:`~.MagData` object.
 
     '''
+    import jutil
     LOG.debug('Calling optimize_cg')
-    mag_0 = first_guess
+    if first_guess is None:
+        first_guess = MagData(data.a, np.zeros((3,)+data.dim))
     x_0 = first_guess.mag_vec
-    y = data.phase_vec
-    kernel = Kernel(data.a, data.dim_uv)
-    fwd_model = ForwardModel(data.projectors, kernel, data.b_0)
-    cost = Costfunction(y, fwd_model)
-    # Optimize:
-    result = minimize(cost, x_0, method='Newton-CG', jac=cost.jac, hessp=cost.hess_dot,
-                      options={'maxiter': 1000, 'disp': True})
-    # Create optimized MagData object:
+    fwd_model = ForwardModel(data)
+    cost = Costfunction(data, regularisator)
+    assert len(x_0) == cost.n, (len(x_0), cost.m, cost.n)
+    result = jutil.minimizer.minimize(cost, x_0, options={"conv_rel":1e-20}, tol={"max_iteration":1})
     x_opt = result.x
-    mag_opt = MagData(mag_0.a, np.zeros((3,)+mag_0.dim))
+    print cost(x_opt)
+    mag_opt = MagData(data.a, np.zeros((3,)+data.dim))
     mag_opt.mag_vec = x_opt
     return mag_opt
 
@@ -236,6 +234,8 @@ def optimize_simple_leastsq(phase_map, mask, b_0=1, lam=1E-4, order=0):
 
     J_DICT = [J_0, J_1]  # list of cost-functions with different regularisations
     # Reconstruct the magnetization components:
+    # TODO Use jutil.minimizer.minimize(jutil.costfunction.LeastSquaresCostFunction(J_DICT[order],
+    # ...) or a simpler frontend.
     x_rec, _ = leastsq(J_DICT[order], np.zeros(3*count))
     mag_data_rec.set_vector(x_rec, mask)
     return mag_data_rec

pyramid/regularisator.py

@@ -52,7 +52,7 @@ class Regularisator(object):
     def jac(self, x):
         # TODO: Docstring!
         self.LOG.debug('Calling jac')
-        return self.lam * self.norm.jac_dot(x)
+        return self.lam * self.norm.jac(x)
 
     def hess_dot(self, x, vector):
         # TODO: Docstring!
@@ -105,7 +105,7 @@ class ZeroOrderRegularisator(Regularisator):
 
     def __init__(self, lam):
         self.LOG.debug('Calling __init__')
-        norm = jnorm.NormL2Square()
+        norm = jnorm.L2Square()
         super(ZeroOrderRegularisator, self).__init__(norm, lam)
         self.LOG.debug('Created '+str(self))
 

scripts/collaborations/example_joern.py

@@ -28,6 +28,7 @@ LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)),
 
 
 logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+logging.basicConfig(level=logging.INFO)
 ###################################################################################################
 threshold = 1
 a = 1.0  # in nm
@@ -40,12 +41,13 @@ smoothed_pictures = True
 lam = 1E-4
 order = 1
 log = True
-PATH = '../../output/joern/'
+PATH = './'
+dirname = PATH
 ###################################################################################################
 # Read in files:
 phase_map = PhaseMap.load_from_netcdf4(PATH+'phase_map.nc')
 #mask = np.genfromtxt('mask.txt', dtype=bool)
-with open(PATH+'mask.pickle') as pf:
+with open(PATH + 'mask.pickle') as pf:
     mask = pickle.load(pf)
 # Setup:
 data_set = DataSet(a, dim, b_0)
@@ -53,29 +55,40 @@ data_set.append(phase_map, SimpleProjector(dim))
 regularisator = ZeroOrderRegularisator(lam)
 # Reconstruct the magnetic distribution:
 tic = clock()
-mag_data_rec = rc.optimize_sparse_cg(data_set, regularisator, maxiter=100, verbosity=2)
+mag_data_rec1 = rc.optimize_linear(data_set, regularisator=regularisator)
+mag_data_rec = rc.optimize_nonlin(data_set, regularisator=regularisator)
 #  .optimize_simple_leastsq(phase_map, mask, b_0, lam=lam, order=order)
+
 print 'reconstruction time:', clock() - tic
 # Display the reconstructed phase map and holography image:
 phase_map_rec = pm(mag_data_rec)
-phase_map_rec.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter)
+phase_map_rec.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter, show=False)
+plt.savefig(dirname + "/reconstr.png")
+
 # Plot the magnetization:
-axis = (mag_data_rec*(1/mag_data_rec.magnitude.max())).quiver_plot()
+axis = (mag_data_rec*(1/mag_data_rec.magnitude.max())).quiver_plot(show=False)
 axis.set_xlim(20, 45)
 axis.set_ylim(20, 45)
+plt.savefig(dirname + "/quiver.png")
+
 # Display the Phase:
 phase_diff = phase_map_rec-phase_map
-phase_diff.display_phase('Difference')
+phase_diff.display_phase('Difference', show=False)
+plt.savefig(dirname + "/difference.png")
+
 # Get the average difference from the experimental results:
 print 'Average difference:', np.average(phase_diff.phase)
 # Plot holographic contour maps with overlayed magnetic distributions:
-axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
-mag_data_rec.quiver_plot(axis=axis)
+axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter, show=False)
+mag_data_rec.quiver_plot(axis=axis, show=False)
 axis = plt.gca()
 axis.set_xlim(20, 45)
 axis.set_ylim(20, 45)
-axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
-mag_data_rec.quiver_plot(axis=axis, log=log)
+plt.savefig(dirname + "/overlay_normal.png")
+
+axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter, show=False)
+mag_data_rec.quiver_plot(axis=axis, log=log, show=False)
 axis = plt.gca()
 axis.set_xlim(20, 45)
 axis.set_ylim(20, 45)
+plt.savefig(dirname + "/overlay_log.png")