From 74619d301bda35b312ffa959a3e353b01ebeffe9 Mon Sep 17 00:00:00 2001
From: "Joern Ungermann (IEK-7 FZ-Juelich)" <j.ungermann@fz-juelich.de>
Date: Tue, 9 Jul 2013 13:10:46 +0200
Subject: [PATCH] Added draft for a numcore routine performing a faster real
space forward calculation.
---
pyramid/magdata.py | 3 ++-
pyramid/numcore/__init__.py | 1 +
pyramid/numcore/phase_mag_real.pyx | 32 ++++++++++++++++++++++++++++++
pyramid/phasemapper.py | 10 +++++++++-
setup.py | 6 +++---
5 files changed, 47 insertions(+), 5 deletions(-)
create mode 100644 pyramid/numcore/phase_mag_real.pyx
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index 797a14e..61221c7 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -5,7 +5,6 @@
import numpy as np
import tables.netcdf3 as nc
import matplotlib.pyplot as plt
-from mayavi import mlab
class MagData:
@@ -164,6 +163,8 @@ class MagData:
None
'''
+ from mayavi import mlab
+
res = self.res
dim = self.dim
# Create points and vector components as lists:
diff --git a/pyramid/numcore/__init__.py b/pyramid/numcore/__init__.py
index e69de29..c4a3052 100644
--- a/pyramid/numcore/__init__.py
+++ b/pyramid/numcore/__init__.py
@@ -0,0 +1 @@
+from phase_mag_real import *
diff --git a/pyramid/numcore/phase_mag_real.pyx b/pyramid/numcore/phase_mag_real.pyx
new file mode 100644
index 0000000..efb838f
--- /dev/null
+++ b/pyramid/numcore/phase_mag_real.pyx
@@ -0,0 +1,32 @@
+import numpy as np
+import math
+cimport cython
+cimport numpy as np
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def phase_mag_real_helper_1(
+ unsigned int v_dim, unsigned int u_dim,
+ double[:, :] phi_u,
+ double[:, :] phi_v,
+ double[:, :] u_mag,
+ double[:, :] v_mag,
+ double[:, :] phase, float threshold):
+ cdef unsigned int i, j, ii, jj, iii, jjj
+ cdef double u, v
+ for j in range(v_dim):
+ for i in range(u_dim):
+ u = u_mag[j, i]
+ v = v_mag[j, i]
+ iii = u_dim - 1 - i
+ jjj = v_dim - 1 - j
+ if abs(u) > threshold:
+ for jj in range(phase.shape[0]):
+ for ii in range(phase.shape[1]):
+ phase[jj, ii] += u * phi_u[jjj + jj, iii + ii]
+ if abs(v) > threshold:
+ for jj in range(phase.shape[0]):
+ for ii in range(phase.shape[1]):
+ phase[jj, ii] -= v * phi_v[jjj + jj, iii + ii]
+
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 2939b35..7d4a3f0 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -3,7 +3,8 @@
import numpy as np
-
+import numcore
+import time
# Physical constants
PHI_0 = -2067.83 # magnetic flux in T*nm²
@@ -102,12 +103,19 @@ def phase_mag_real(res, projection, method, b_0=1, jacobi=None):
phase -= v_mag[j, i] * phase_v
############################### TODO: NUMERICAL CORE #####################################
else: # Without Jacobi matrix (faster)
+## phasecopy = phase.copy()
+## start_time = time.time()
+## numcore.phase_mag_real_helper_1(v_dim, u_dim, phi_u, phi_v, u_mag, v_mag, phasecopy, threshold)
+## print time.time() - start_time
+## start_time = time.time()
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]
+## print time.time() - start_time
+## print ((phase - phasecopy) ** 2).sum()
# Return the phase:
return phase
diff --git a/setup.py b/setup.py
index 48a121a..b4a8f41 100644
--- a/setup.py
+++ b/setup.py
@@ -15,10 +15,10 @@ from Cython.Build import cythonize
setup(
name = 'Pyramid',
- version = '0.1',
+ version = '0.1',
description = 'PYthon based Reconstruction Algorithm for MagnetIc Distributions',
author = 'Jan Caron',
author_email = 'j.caron@fz-juelich.de',
- packages = ['pyramid'],
- ext_modules = cythonize(glob.glob(os.path.join('pyramid','numcore','*.pyx')))
+ packages = ['pyramid', 'pyramid.numcore'],
+ ext_modules = cythonize(glob.glob(os.path.join('pyramid', 'numcore', '*.pyx')))
)
--
GitLab