Author: tlinnet
Date: Mon Jun 23 14:09:27 2014
New Revision: 24243
URL: http://svn.gna.org/viewcvs/relax?rev=24243&view=rev
Log:
Added unit test for doing the matrix exponential for complex data.
This test shows, that the dtype=complex64, should be removed from
lib/dispersion/ns_mmq_2site.py.
Task #7807 (https://gna.org/task/index.php?7807): Speed-up of dispersion
models for Clustered analysis.
Modified:
branches/disp_spin_speed/test_suite/unit_tests/_lib/_dispersion/test_matrix_exponential.py
Modified:
branches/disp_spin_speed/test_suite/unit_tests/_lib/_dispersion/test_matrix_exponential.py
URL:
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/test_suite/unit_tests/_lib/_dispersion/test_matrix_exponential.py?rev=24243&r1=24242&r2=24243&view=diff
==============================================================================
---
branches/disp_spin_speed/test_suite/unit_tests/_lib/_dispersion/test_matrix_exponential.py
(original)
+++
branches/disp_spin_speed/test_suite/unit_tests/_lib/_dispersion/test_matrix_exponential.py
Mon Jun 23 14:09:27 2014
@@ -21,11 +21,11 @@
# Python module imports.
from os import path, sep
-from numpy import dot, einsum, load, sum
+from numpy import complex64, complex128, dot, einsum, load, sum
from unittest import TestCase
# relax module imports.
-from lib.dispersion.ns_matrices import rcpmg_3d_rankN
+from lib.dispersion.ns_matrices import rcpmg_3d_rankN, rmmq_2site_rankN
from lib.linear_algebra.matrix_exponential import matrix_exponential
from lib.dispersion.matrix_exponential import matrix_exponential_rankN
from status import Status; status = Status()
@@ -67,6 +67,30 @@
return r180x, M0, r10a, r10b, r20a, r20b, pA, dw, dw_orig, kex,
inv_tcpmg, tcp, num_points, power, back_calc, pB, k_BA, k_AB
+ def return_data_mmq_2site(self, fname):
+ """Return the data structures from the data path."""
+
+ M0 = load(fname+"_M0"+".npy")
+ r20a = load(fname+"_r20a"+".npy")
+ r20b = load(fname+"_r20b"+".npy")
+ pA = load(fname+"_pA"+".npy")
+ dw = load(fname+"_dw"+".npy")
+ dwH = load(fname+"_dwH"+".npy")
+ kex = load(fname+"_kex"+".npy")
+ inv_tcpmg = load(fname+"_inv_tcpmg"+".npy")
+ tcp = load(fname+"_tcp"+".npy")
+ num_points = load(fname+"_num_points"+".npy")
+ power = load(fname+"_power"+".npy")
+ back_calc = load(fname+"_back_calc"+".npy")
+
+ # Once off parameter conversions.
+ pB = 1.0 - pA
+ k_BA = pA * kex
+ k_AB = pB * kex
+
+ return M0, r20a, r20b, pA, dw, dwH, kex, inv_tcpmg, tcp, num_points,
power, back_calc, pB, k_BA, k_AB
+
+
def test_ns_cpmg_2site_3d_hansen_cpmg_data(self):
"""Test the matrix_exponential_rankN() function for higher
dimensional data, and compare to matrix_exponential. This uses the data from
systemtest Relax_disp.test_hansen_cpmg_data_to_ns_cpmg_2site_3D."""
@@ -106,3 +130,114 @@
# Test that the sum difference is zero.
self.assertAlmostEqual(diff_sum, 0.0)
+
+ def test_ns_mmq_2site_korzhnev_2005_15n_dq_data_complex64(self):
+ """Test the matrix_exponential_rankN() function for higher
dimensional data, and compare to matrix_exponential. This uses the data from
systemtest Relax_disp.test_korzhnev_2005_15n_dq_data.
+ This test does the matrix exponential in complex64."""
+
+ fname = self.data + sep+ "test_korzhnev_2005_15n_dq_data"
+ M0, R20A, R20B, pA, dw, dwH, kex, inv_tcpmg, tcp, num_points, power,
back_calc, pB, k_BA, k_AB = self.return_data_mmq_2site(fname)
+
+ # Extract the total numbers of experiments, number of spins, number
of magnetic field strength, number of offsets, maximum number of dispersion
point.
+ NS, NM, NO = num_points.shape
+
+ # Populate the m1 and m2 matrices (only once per function call for
speed).
+ m1_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=dw, k_AB=k_AB,
k_BA=k_BA, tcp=tcp)
+ m2_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=-dw, k_AB=k_AB,
k_BA=k_BA, tcp=tcp)
+
+ # The A+/- matrices.
+ A_pos_mat = matrix_exponential_rankN(m1_mat, dtype=complex64)
+ A_neg_mat = matrix_exponential_rankN(m2_mat, dtype=complex64)
+
+ # Loop over spins.
+ for si in range(NS):
+ # Loop over the spectrometer frequencies.
+ for mi in range(NM):
+ # Loop over offsets:
+ for oi in range(NO):
+ # Extract number of points.
+ num_points_i = num_points[si, mi, oi]
+
+ # Loop over the time points, back calculating the R2eff
values.
+ for i in range(num_points_i):
+ # Test the two different methods.
+ # The A+/- matrices.
+ A_pos_i = A_pos_mat[si, mi, oi, i]
+ A_neg_i = A_neg_mat[si, mi, oi, i]
+
+ # The lower dimensional matrix exponential.
+ A_pos = matrix_exponential(m1_mat[si, mi, oi, i])
+ A_neg = matrix_exponential(m2_mat[si, mi, oi, i])
+
+ # Calculate differences
+ diff_A_pos_real = A_pos_i.real - A_pos.real
+ diff_A_pos_real_sum = sum(diff_A_pos_real)
+ diff_A_pos_imag = A_pos_i.imag - A_pos.imag
+ diff_A_pos_imag_sum = sum(diff_A_pos_imag)
+
+ diff_A_neg_real = A_neg_i.real - A_neg.real
+ diff_A_neg_real_sum = sum(diff_A_neg_real)
+ diff_A_neg_imag = A_neg_i.imag - A_neg.imag
+ diff_A_neg_imag_sum = sum(diff_A_neg_imag)
+
+ # Test that the sum difference is zero.
+ self.assertAlmostEqual(diff_A_pos_real_sum , 0.0, 6)
+ self.assertAlmostEqual(diff_A_pos_imag_sum , 0.0, 6)
+ self.assertAlmostEqual(diff_A_neg_real_sum , 0.0, 6)
+ self.assertAlmostEqual(diff_A_neg_imag_sum , 0.0, 6)
+
+
+ def test_ns_mmq_2site_korzhnev_2005_15n_dq_data_complex128(self):
+ """Test the matrix_exponential_rankN() function for higher
dimensional data, and compare to matrix_exponential. This uses the data from
systemtest Relax_disp.test_korzhnev_2005_15n_dq_data.
+ This test does the matrix exponential in complex128."""
+
+ fname = self.data + sep+ "test_korzhnev_2005_15n_dq_data"
+ M0, R20A, R20B, pA, dw, dwH, kex, inv_tcpmg, tcp, num_points, power,
back_calc, pB, k_BA, k_AB = self.return_data_mmq_2site(fname)
+
+ # Extract the total numbers of experiments, number of spins, number
of magnetic field strength, number of offsets, maximum number of dispersion
point.
+ NS, NM, NO = num_points.shape
+
+ # Populate the m1 and m2 matrices (only once per function call for
speed).
+ m1_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=dw, k_AB=k_AB,
k_BA=k_BA, tcp=tcp)
+ m2_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=-dw, k_AB=k_AB,
k_BA=k_BA, tcp=tcp)
+
+ # The A+/- matrices.
+ A_pos_mat = matrix_exponential_rankN(m1_mat)
+ A_neg_mat = matrix_exponential_rankN(m2_mat)
+
+ # Loop over spins.
+ for si in range(NS):
+ # Loop over the spectrometer frequencies.
+ for mi in range(NM):
+ # Loop over offsets:
+ for oi in range(NO):
+ # Extract number of points.
+ num_points_i = num_points[si, mi, oi]
+
+ # Loop over the time points, back calculating the R2eff
values.
+ for i in range(num_points_i):
+ # Test the two different methods.
+ # The A+/- matrices.
+ A_pos_i = A_pos_mat[si, mi, oi, i]
+ A_neg_i = A_neg_mat[si, mi, oi, i]
+
+ # The lower dimensional matrix exponential.
+ A_pos = matrix_exponential(m1_mat[si, mi, oi, i])
+ A_neg = matrix_exponential(m2_mat[si, mi, oi, i])
+
+ # Calculate differences
+ diff_A_pos_real = A_pos_i.real - A_pos.real
+ diff_A_pos_real_sum = sum(diff_A_pos_real)
+ diff_A_pos_imag = A_pos_i.imag - A_pos.imag
+ diff_A_pos_imag_sum = sum(diff_A_pos_imag)
+
+ diff_A_neg_real = A_neg_i.real - A_neg.real
+ diff_A_neg_real_sum = sum(diff_A_neg_real)
+ diff_A_neg_imag = A_neg_i.imag - A_neg.imag
+ diff_A_neg_imag_sum = sum(diff_A_neg_imag)
+
+ # Test that the sum difference is zero.
+ self.assertEqual(diff_A_pos_real_sum , 0.0)
+ self.assertEqual(diff_A_pos_imag_sum , 0.0)
+ self.assertEqual(diff_A_neg_real_sum , 0.0)
+ self.assertEqual(diff_A_neg_imag_sum , 0.0)
r24243 - /branches/disp_spin_speed/test_suite/unit_tests/_lib/_dispersion/test_matrix_exponential.py