From c9e6285f8b5a30ed94d035ba25e8c4253dcc65b3 Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Fri, 24 May 2013 08:56:50 +0200
Subject: [PATCH] phasemapper: Cleaned up the two approaches (mx and my, m and
beta)
---
pyramid/phasemapper.py | 135 +++++++++++++++-----------------
pyramid/test/run_tests.py | 2 +-
pyramid/test/test_compliance.py | 7 +-
scripts/get_jacobi.py | 14 ++--
scripts/invert2.py | 115 +++++++++++++++++++++++++++
5 files changed, 190 insertions(+), 83 deletions(-)
create mode 100644 scripts/invert2.py
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 69a48e8..064a568 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -62,6 +62,7 @@ def phase_mag_real(res, projection, method, b_0=1, jacobi=None):
the phasemap as a 2 dimensional array
'''
+ # Function for creating the lookup-tables:
def phi_lookup(method, n, m, res, b_0):
if method == 'slab':
def F_h(n, m):
@@ -73,13 +74,6 @@ def phase_mag_real(res, projection, method, b_0=1, jacobi=None):
elif method == 'disc':
in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
-# def phi_mag(i, j): # TODO: rename
-# return (np.cos(beta[j,i])*phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-# -np.sin(beta[j,i])*phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
-#
-# def phi_mag_deriv(i, j): # TODO: rename
-# return -(np.sin(beta[j,i])*phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-# +np.cos(beta[j,i])*phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
# Process input parameters:
v_dim, u_dim = np.shape(projection[0])
v_mag, u_mag = projection
@@ -90,22 +84,35 @@ def phase_mag_real(res, projection, method, b_0=1, jacobi=None):
uu, vv = np.meshgrid(u, v)
phi_u = phi_lookup(method, uu, vv, res, b_0)
phi_v = phi_lookup(method, vv, uu, res, b_0)
- '''CALCULATE THE PHASE'''
+ # Calculation of the phase:
phase = np.zeros((v_dim, u_dim))
- for j in range(v_dim): # TODO: only iterate over pixels that have a magn. > threshold (first >0)
- for i in range(u_dim):
- if (u_mag[j,i] != 0 and v_mag[j,i] != 0) or jacobi is not None: # TODO: same result with or without?
- phase_u = phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+ threshold = 0
+ if jacobi is not None: # With Jacobian matrix (slower)
+ jacobi[:] = 0 # Jacobi matrix --> zeros
+ ############################### TODO: NUMERICAL CORE #####################################
+ for j in range(v_dim):
+ for i in range(u_dim):
+ phase_u = phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+ jacobi[:,i+u_dim*j] = phase_u.reshape(-1)
+ if abs(u_mag[j, i]) > threshold:
+ phase += u_mag[j,i] * phase_u
phase_v = phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
- phase += u_mag[j,i]*phase_u - v_mag[j,i]*phase_v
- if jacobi is not None:
- jacobi[:,i+u_dim*j] = phase_u.reshape(-1)
- jacobi[:,u_dim*v_dim+i+u_dim*j] = -phase_v.reshape(-1)
-
+ jacobi[:,u_dim*v_dim+i+u_dim*j] = -phase_v.reshape(-1)
+ if abs(v_mag[j, i]) > threshold:
+ phase -= v_mag[j,i] * phase_v
+ ############################### TODO: NUMERICAL CORE #####################################
+ else: # Without Jacobi matrix (faster)
+ for j in range(v_dim):
+ for i in range(u_dim):
+ if abs(u_mag[j, i]) > threshold:
+ phase += u_mag[j, i] * phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+ if abs(v_mag[j, i]) > threshold:
+ phase -= v_mag[j, i] * phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+ # Return the phase:
return phase
-def phase_mag_real_ANGLE(res, projection, method, b_0=1, jacobi=None): # TODO: Modify
+def phase_mag_real_alt(res, projection, method, b_0=1, jacobi=None): # TODO: Modify
'''Calculate phasemap from magnetization data (real space approach).
Arguments:
res - the resolution of the grid (grid spacing) in nm
@@ -117,78 +124,58 @@ def phase_mag_real_ANGLE(res, projection, method, b_0=1, jacobi=None): # TODO:
the phasemap as a 2 dimensional array
'''
+ # Function for creating the lookup-tables:
def phi_lookup(method, n, m, res, b_0):
if method == 'slab':
def F_h(n, m):
a = np.log(res**2 * (n**2 + m**2))
b = np.arctan(n / m)
return n*a - 2*n + 2*m*b
- coeff = -b_0 * res**2 / ( 2 * PHI_0 )
return coeff * 0.5 * ( F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
-F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5) )
elif method == 'disc':
- coeff = - b_0 * res**2 / ( 2 * PHI_0 )
in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
-
+ # Function for the phase contribution of one pixel:
+ def phi_mag(i, j):
+ return (np.cos(beta[j,i])*phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+ -np.sin(beta[j,i])*phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
+ # Function for the derivative of the phase contribution of one pixel:
+ def phi_mag_deriv(i, j):
+ return -(np.sin(beta[j,i])*phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+ +np.cos(beta[j,i])*phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
+ # Process input parameters:
v_dim, u_dim = np.shape(projection[0])
v_mag, u_mag = projection
-
beta = np.arctan2(v_mag, u_mag)
- mag = np.hypot(u_mag, v_mag)
-
- '''CREATE COORDINATE GRIDS'''
- u = np.linspace(0,(u_dim-1),num=u_dim)
- v = np.linspace(0,(v_dim-1),num=v_dim)
- uu, vv = np.meshgrid(u,v)
-
- u_lookup = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
- v_lookup = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
- uu_lookup, vv_lookup = np.meshgrid(u_lookup, v_lookup)
-
- phi_cos = phi_lookup(method, uu_lookup, vv_lookup, res, b_0)
- phi_sin = phi_lookup(method, vv_lookup, uu_lookup, res, b_0)
-
- def phi_mag(i, j): # TODO: rename
- return (np.cos(beta[j,i])*phi_cos[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
- -np.sin(beta[j,i])*phi_sin[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
-
- def phi_mag_deriv(i, j): # TODO: rename
- return -(np.sin(beta[j,i])*phi_cos[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
- +np.cos(beta[j,i])*phi_sin[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
-
- def phi_mag_fd(i, j, h): # TODO: rename
- return ((np.cos(beta[j,i]+h) - np.cos(beta[j,i])) / h
- * phi_cos[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
- -(np.sin(beta[j,i]+h) - np.sin(beta[j,i])) / h
- * phi_sin[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
-
- '''CALCULATE THE PHASE'''
+ mag = np.hypot(u_mag, v_mag)
+ coeff = -b_0 * res**2 / ( 2 * PHI_0 )
+ # Create lookup-tables for the phase of one pixel:
+ u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
+ v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
+ uu, vv = np.meshgrid(u, v)
+ phi_u = phi_lookup(method, uu, vv, res, b_0)
+ phi_v = phi_lookup(method, vv, uu, res, b_0)
+ # Calculation of the phase:
phase = np.zeros((v_dim, u_dim))
-
- # TODO: only iterate over pixels that have a magn. > threshold (first >0)
- if jacobi is not None:
- jacobi_fd = jacobi.copy()
- h = 0.0001
-
- for j in range(v_dim):
- for i in range(u_dim):
- #if (mag[j,i] != 0 ):#or jacobi is not None): # TODO: same result with or without?
- phi_mag_cache = phi_mag(i, j)
- phase += mag[j,i] * phi_mag_cache
- if jacobi is not None:
- jacobi[:,i+u_dim*j] = phi_mag_cache.reshape(-1)
+ threshold = 0
+ if jacobi is not None: # With Jacobian matrix (slower)
+ jacobi[:] = 0 # Jacobi matrix --> zeros
+ ############################### TODO: NUMERICAL CORE #####################################
+ for j in range(v_dim):
+ for i in range(u_dim):
+ phase_cache = phi_mag(i, j)
+ jacobi[:,i+u_dim*j] = phase_cache.reshape(-1)
+ if mag[j, i] > threshold:
+ phase += mag[j,i]*phase_cache
jacobi[:,u_dim*v_dim+i+u_dim*j] = (mag[j,i]*phi_mag_deriv(i,j)).reshape(-1)
-
- jacobi_fd[:,i+u_dim*j] = phi_mag_cache.reshape(-1)
- jacobi_fd[:,u_dim*v_dim+i+u_dim*j] = (mag[j,i]*phi_mag_fd(i,j,h)).reshape(-1)
-
-
-
- if jacobi is not None:
- jacobi_diff = jacobi_fd - jacobi
- assert (np.abs(jacobi_diff) < 1.0E-8).all(), 'jacobi matrix is not the same'
-
+ ############################### TODO: NUMERICAL CORE #####################################
+ else: # Without Jacobi matrix (faster)
+ for j in range(v_dim):
+ for i in range(u_dim):
+ if abs(mag[j, i]) > threshold:
+ phase += mag[j,i] * phi_mag(i, j)
+ # Return the phase:
return phase
diff --git a/pyramid/test/run_tests.py b/pyramid/test/run_tests.py
index fbe6fec..a1e2a62 100644
--- a/pyramid/test/run_tests.py
+++ b/pyramid/test/run_tests.py
@@ -31,7 +31,7 @@ def run():
suite.addTest(loader.loadTestsFromTestCase(TestCaseAnalytic))
suite.addTest(loader.loadTestsFromTestCase(TestCaseReconstructor))
runner.run(suite)
-
+
if __name__ == '__main__':
run()
\ No newline at end of file
diff --git a/pyramid/test/test_compliance.py b/pyramid/test/test_compliance.py
index bbe3e58..f6d438f 100644
--- a/pyramid/test/test_compliance.py
+++ b/pyramid/test/test_compliance.py
@@ -51,10 +51,15 @@ class TestCaseCompliance(unittest.TestCase):
self.assertTrue(self.pep8 is not None,
msg="Install Python pep8 module to fully execute testbench!")
errors = 0
- for dir_name in ["test", os.environ["BUILDDIR"]]:
+ for dir_name in "./":
errors += self.checkDirectory(dir_name)
self.assertEqual(errors, 0)
+
+if __name__ == '__main__':
+ suite = unittest.TestLoader().loadTestsFromTestCase(TestCaseCompliance)
+ unittest.TextTestRunner(verbosity=2).run(suite)
+
# def test_variables(self):
# """
# Function for checking that all attributes are present.
diff --git a/scripts/get_jacobi.py b/scripts/get_jacobi.py
index 26732ce..8df0613 100644
--- a/scripts/get_jacobi.py
+++ b/scripts/get_jacobi.py
@@ -10,12 +10,13 @@ import pyramid.magcreator as mc
import pyramid.projector as pj
import pyramid.phasemapper as pm
from pyramid.magdata import MagData
+from pyramid.phasemap import PhaseMap
import time
import pdb, traceback, sys
from numpy import pi
-def phase_from_mag():
+def get_jacobi():
'''Calculate and display the phase map from a given magnetization.
Arguments:
None
@@ -27,22 +28,21 @@ def phase_from_mag():
b_0 = 1.0 # in T
dim = (1, 3, 3) # in px (y,x)
res = 10.0 # in nm
- beta = 0*pi/4
+ beta = pi/4
center = (0, 1, 1) # in px (y,x) index starts with 0!
width = (0, 1, 1) # in px (y,x)
- mag_shape = mc.Shapes.slab(dim, center, width)
- mag_data = MagData(res, mc.create_mag_dist(mag_shape, beta))
-
+ mag_data = MagData(res, mc.create_mag_dist(mc.Shapes.slab(dim, center, width), beta))
projection = pj.simple_axis_projection(mag_data)
'''NUMERICAL SOLUTION'''
# numerical solution Real Space (Slab):
jacobi = np.zeros((dim[2]*dim[1], 2*dim[2]*dim[1]))
tic = time.clock()
- pm.phase_mag_real(res, projection, 'slab', b_0, jacobi=jacobi)
+ phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab', b_0, jacobi=jacobi))
toc = time.clock()
+ phase_map.display()
np.savetxt('../output/jacobi.npy', jacobi)
print 'Time for Real Space Approach with Jacobi-Matrix (Slab): ' + str(toc - tic)
@@ -51,7 +51,7 @@ def phase_from_mag():
if __name__ == "__main__":
try:
- jacobi = phase_from_mag()
+ jacobi = get_jacobi()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
diff --git a/scripts/invert2.py b/scripts/invert2.py
new file mode 100644
index 0000000..5837248
--- /dev/null
+++ b/scripts/invert2.py
@@ -0,0 +1,115 @@
+# -*- coding: utf-8 -*-
+"""Create random magnetic distributions."""
+
+import random as rnd
+import pdb, traceback, sys
+import numpy as np
+from numpy import pi
+import pylab
+from copy import deepcopy
+
+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
+from scipy.optimize import leastsq
+
+
+def create_random_distribution():
+ '''Calculate, display and save a random magnetic distribution to file.
+ Arguments:
+ None
+ Returns:
+ None
+
+ '''
+ # Input parameters:
+ count = 200
+ dim = (1, 32, 32)
+ b_0 = 1 # in T
+ res = 10 # in nm
+ rnd.seed(12)
+ # Create lists for magnetic objects:
+ mag_shape_list = np.zeros((count,) + dim)
+ beta_list = np.zeros(count)
+ magnitude_list = np.zeros(count)
+ for i in range(count):
+ pixel = (rnd.randrange(dim[0]), rnd.randrange(dim[1]), rnd.randrange(dim[2]))
+ mag_shape_list[i,...] = mc.Shapes.pixel(dim, pixel)
+ beta_list[i] = 2*pi*rnd.random()
+ magnitude_list[i] = rnd.random()
+ # Create magnetic distribution
+ magnitude = mc.create_mag_dist_comb(mag_shape_list, beta_list, magnitude_list)
+ mag_data = MagData(res, magnitude)
+# mag_data.quiver_plot()
+# mag_data.save_to_llg('output/mag_dist_random_pixel.txt')
+
+
+ magx, magy, magz = mag_data.magnitude
+ maskx, masky, maskz = magx != 0, magy != 0, magz !=0
+ x_t = np.concatenate([magx[maskx], magy[masky], magz[maskz]])
+ print "x_t", x_t
+
+ def F(x):
+ mag_data_temp = MagData(deepcopy(mag_data.res), deepcopy(mag_data.magnitude))
+ magx, magy, magz = mag_data_temp.magnitude
+ maskx, masky, maskz = magx != 0, magy != 0, magz !=0
+# print maskx.sum() + masky.sum() + maskz.sum()
+ assert len(x) == maskx.sum() + masky.sum() + maskz.sum()
+ magx[maskx] = x[:maskx.sum()]
+ magy[masky] = x[maskx.sum():maskx.sum() + masky.sum()]
+ magz[maskz] = x[maskx.sum() + masky.sum():]
+ projection = pj.simple_axis_projection(mag_data_temp)
+ phase_map_slab = PhaseMap(res, pm.phase_mag_real_ANGLE(res, projection, 'slab', b_0))
+ return phase_map_slab.phase.reshape(-1)
+
+ y_m = F(x_t)
+ print "y_m", y_m
+
+ lam = 1e-6
+
+ def J(x_i):
+# print "x_i", x_i
+ y_i = F(x_i)
+# dd1 = np.zeros(mx.shape)
+# dd2 = np.zeros(mx.shape)
+ term1 = (y_i - y_m)
+ term2 = lam * x_i
+# dd1[:, :-1, :] += np.diff(mx, axis=1)
+# dd1[:, -1, :] += np.diff(mx, axis=1)[:, -1, :]
+# dd1[:, :, :-1] += np.diff(mx, axis=2)
+# dd1[:, :, -1] += np.diff(mx, axis=2)[:, :, -1]
+
+# dd2[:, :-1, :] += np.diff(my, axis=1)
+# dd2[:, -1, :] += np.diff(my, axis=1)[:, -1, :]
+# dd2[:, :, :-1] += np.diff(my, axis=2)
+# dd2[:, :, -1] += np.diff(my, axis=2)[:, :, -1]
+
+# result = np.concatenate([term1, np.sqrt(abs(dd1.reshape(-1))), np.sqrt(abs(dd2.reshape(-1)))])
+# result = np.concatenate([term1, np.sqrt(abs(dd1.reshape(-1)))])
+# print result
+ return np.concatenate([term1, term2])
+
+ x_f, _ = leastsq(J, np.zeros(x_t.shape))
+ y_f = F(x_f)
+# print "y_m", y_m
+# print "y_f", y_f
+# print "dy", y_f - y_m
+# print "x_t", x_t
+# print "x_f", x_f
+# print "dx", x_f - x_t
+# pylab.show()
+ projection = pj.simple_axis_projection(mag_data)
+ phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'))
+ hi.display(hi.holo_image(phase_map, 10))
+
+
+if __name__ == "__main__":
+ try:
+ create_random_distribution()
+ except:
+ type, value, tb = sys.exc_info()
+ traceback.print_exc()
+ pdb.post_mortem(tb)
--
GitLab