Source Code for Module test_suite.unit_tests._lib._dispersion.test_matrix_exponential

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2014 Troels E. Linnet                                         # 
  4  # Copyright (C) 2014 Edward d'Auvergne                                        # 
  5  #                                                                             # 
  6  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  7  #                                                                             # 
  8  # This program is free software: you can redistribute it and/or modify        # 
  9  # it under the terms of the GNU General Public License as published by        # 
 10  # the Free Software Foundation, either version 3 of the License, or           # 
 11  # (at your option) any later version.                                         # 
 12  #                                                                             # 
 13  # This program is distributed in the hope that it will be useful,             # 
 14  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 15  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 16  # GNU General Public License for more details.                                # 
 17  #                                                                             # 
 18  # You should have received a copy of the GNU General Public License           # 
 19  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23  # Python module imports. 
 24  from os import sep 
 25  from numpy import complex64, load, sum 
 26  from unittest import TestCase 
 27   
 28  # relax module imports. 
 29  from lib.dispersion.ns_cpmg_2site_3d import rcpmg_3d_rankN 
 30  from lib.dispersion.ns_mmq_2site import rmmq_2site_rankN 
 31  from lib.linear_algebra.matrix_exponential import matrix_exponential as np_matrix_exponential 
 32  from lib.dispersion.matrix_exponential import matrix_exponential 
 33  from status import Status; status = Status() 
 34   
 35   
 36 -class Test_matrix_exponential(TestCase): 
 37      """Unit tests for the lib.dispersion.matrix_exponential relax module.""" 
 38   
 39 -    def setUp(self): 
 40          """Set up for all unit tests.""" 
 41   
 42          # Default parameter values. 
 43          self.data = status.install_path+sep+'test_suite'+sep+'shared_data'+sep+'dispersion'+sep+'unit_tests'+sep+'lib'+sep+'dispersion'+sep+'matrix_exponential' 
 44   
 45 -    def return_data_ns_cpmg_2site_3d(self, fname): 
 46          """Return the data structures from the data path.""" 
 47   
 48          r180x = load(fname+"_r180x"+".npy") 
 49          M0 = load(fname+"_M0"+".npy") 
 50          r10a = load(fname+"_r10a"+".npy") 
 51          r10b = load(fname+"_r10b"+".npy") 
 52          r20a = load(fname+"_r20a"+".npy") 
 53          r20b = load(fname+"_r20b"+".npy") 
 54          pA = load(fname+"_pA"+".npy") 
 55          dw = load(fname+"_dw"+".npy") 
 56          dw_orig = load(fname+"_dw_orig"+".npy") 
 57          kex = load(fname+"_kex"+".npy") 
 58          inv_tcpmg = load(fname+"_inv_tcpmg"+".npy") 
 59          tcp = load(fname+"_tcp"+".npy") 
 60          num_points = load(fname+"_num_points"+".npy") 
 61          power = load(fname+"_power"+".npy") 
 62          back_calc = load(fname+"_back_calc"+".npy") 
 63   
 64          # Once off parameter conversions. 
 65          pB = 1.0 - pA 
 66          k_BA = pA * kex 
 67          k_AB = pB * kex 
 68   
 69          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 
 70   
 71   
 72 -    def return_data_mmq_2site(self, fname): 
 73          """Return the data structures from the data path.""" 
 74   
 75          M0 = load(fname+"_M0"+".npy") 
 76          r20a = load(fname+"_r20a"+".npy") 
 77          r20b = load(fname+"_r20b"+".npy") 
 78          pA = load(fname+"_pA"+".npy") 
 79          dw = load(fname+"_dw"+".npy") 
 80          dwH = load(fname+"_dwH"+".npy") 
 81          kex = load(fname+"_kex"+".npy") 
 82          inv_tcpmg = load(fname+"_inv_tcpmg"+".npy") 
 83          tcp = load(fname+"_tcp"+".npy") 
 84          num_points = load(fname+"_num_points"+".npy") 
 85          power = load(fname+"_power"+".npy") 
 86          back_calc = load(fname+"_back_calc"+".npy") 
 87   
 88          # Once off parameter conversions. 
 89          pB = 1.0 - pA 
 90          k_BA = pA * kex 
 91          k_AB = pB * kex 
 92   
 93          return M0, r20a, r20b, pA, dw, dwH, kex, inv_tcpmg, tcp, num_points, power, back_calc, pB, k_BA, k_AB 
 94   
 95   
 96 -    def test_ns_cpmg_2site_3d_hansen_cpmg_data(self): 
 97          """Test the matrix_exponential() 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.""" 
 98   
 99          fname = self.data + sep+ "test_hansen_cpmg_data_to_ns_cpmg_2site_3D" 
100          r180x, M0, r10a, r10b, r20a, r20b, pA, dw, dw_orig, kex, inv_tcpmg, tcp, num_points, power, back_calc, pB, k_BA, k_AB = self.return_data_ns_cpmg_2site_3d(fname) 
101   
102          # Extract the total numbers of experiments, number of spins, number of magnetic field strength, number of offsets, maximum number of dispersion point. 
103          NE, NS, NM, NO, ND = back_calc.shape 
104   
105          # The matrix R that contains all the contributions to the evolution, i.e. relaxation, exchange and chemical shift evolution. 
106          R_mat = rcpmg_3d_rankN(R1A=r10a, R1B=r10b, R2A=r20a, R2B=r20b, pA=pA, pB=pB, dw=dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 
107       
108          # This matrix is a propagator that will evolve the magnetization with the matrix R for a delay tcp. 
109          Rexpo_mat = matrix_exponential(R_mat) 
110       
111          # Loop over the spins 
112          for si in range(NS): 
113              # Loop over the spectrometer frequencies. 
114              for mi in range(NM): 
115                  # Extract number of points. 
116                  num_points_si_mi = int(num_points[0, si, mi, 0]) 
117       
118                  # Loop over the time points, back calculating the R2eff values. 
119                  for di in range(num_points_si_mi): 
120                      # Test the two different methods. 
121                      R_mat_i = R_mat[0, si, mi, 0, di] 
122    
123                      # The lower dimensional matrix exponential. 
124                      Rexpo = np_matrix_exponential(R_mat_i) 
125       
126                      # The higher dimensional matrix exponential. 
127                      Rexpo_mat_i = Rexpo_mat[0, si, mi, 0, di] 
128   
129                      diff = Rexpo - Rexpo_mat_i 
130                      diff_sum = sum(diff) 
131   
132                      # Test that the sum difference is zero.                                         
133                      self.assertAlmostEqual(diff_sum, 0.0) 
134   
135   
136 -    def test_ns_mmq_2site_korzhnev_2005_15n_dq_data_complex64(self): 
137          """Test the matrix_exponential() function for higher dimensional data, and compare to matrix_exponential.  This uses the data from systemtest Relax_disp.test_korzhnev_2005_15n_dq_data. 
138          This test does the matrix exponential in complex64.""" 
139   
140          fname = self.data + sep+ "test_korzhnev_2005_15n_dq_data" 
141          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) 
142   
143          # Extract the total numbers of experiments, number of spins, number of magnetic field strength, number of offsets, maximum number of dispersion point. 
144          NS, NM, NO = num_points.shape 
145   
146          # Populate the m1 and m2 matrices (only once per function call for speed). 
147          m1_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 
148          m2_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=-dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 
149       
150          # The A+/- matrices. 
151          A_pos_mat = matrix_exponential(m1_mat, dtype=complex64) 
152          A_neg_mat = matrix_exponential(m2_mat, dtype=complex64) 
153       
154          # Loop over spins. 
155          for si in range(NS): 
156              # Loop over the spectrometer frequencies. 
157              for mi in range(NM): 
158                  # Loop over offsets: 
159                  for oi in range(NO): 
160                      # Extract number of points. 
161                      num_points_i = num_points[si, mi, oi] 
162       
163                      # Loop over the time points, back calculating the R2eff values. 
164                      for i in range(num_points_i): 
165                          # Test the two different methods. 
166                          # The A+/- matrices. 
167                          A_pos_i = A_pos_mat[si, mi, oi, i] 
168                          A_neg_i = A_neg_mat[si, mi, oi, i] 
169       
170                          # The lower dimensional matrix exponential. 
171                          A_pos = np_matrix_exponential(m1_mat[si, mi, oi, i]) 
172                          A_neg = np_matrix_exponential(m2_mat[si, mi, oi, i]) 
173       
174                          # Calculate differences 
175                          diff_A_pos_real = A_pos_i.real - A_pos.real 
176                          diff_A_pos_real_sum = sum(diff_A_pos_real) 
177                          diff_A_pos_imag = A_pos_i.imag - A_pos.imag 
178                          diff_A_pos_imag_sum = sum(diff_A_pos_imag) 
179   
180                          diff_A_neg_real = A_neg_i.real - A_neg.real 
181                          diff_A_neg_real_sum = sum(diff_A_neg_real) 
182                          diff_A_neg_imag = A_neg_i.imag - A_neg.imag 
183                          diff_A_neg_imag_sum = sum(diff_A_neg_imag) 
184   
185                          # Test that the sum difference is zero.                                         
186                          self.assertAlmostEqual(diff_A_pos_real_sum, 0.0, 6) 
187                          self.assertAlmostEqual(diff_A_pos_imag_sum, 0.0, 6) 
188                          self.assertAlmostEqual(diff_A_neg_real_sum, 0.0, 6) 
189                          self.assertAlmostEqual(diff_A_neg_imag_sum, 0.0, 6) 
190   
191   
192 -    def test_ns_mmq_2site_korzhnev_2005_15n_dq_data_complex128(self): 
193          """Test the matrix_exponential() function for higher dimensional data, and compare to matrix_exponential.  This uses the data from systemtest Relax_disp.test_korzhnev_2005_15n_dq_data. 
194          This test does the matrix exponential in complex128.""" 
195   
196          fname = self.data + sep+ "test_korzhnev_2005_15n_dq_data" 
197          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) 
198   
199          # Extract the total numbers of experiments, number of spins, number of magnetic field strength, number of offsets, maximum number of dispersion point. 
200          NS, NM, NO = num_points.shape 
201   
202          # Populate the m1 and m2 matrices (only once per function call for speed). 
203          m1_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 
204          m2_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=-dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 
205       
206          # The A+/- matrices. 
207          A_pos_mat = matrix_exponential(m1_mat) 
208          A_neg_mat = matrix_exponential(m2_mat) 
209       
210          # Loop over spins. 
211          for si in range(NS): 
212              # Loop over the spectrometer frequencies. 
213              for mi in range(NM): 
214                  # Loop over offsets: 
215                  for oi in range(NO): 
216                      # Extract number of points. 
217                      num_points_i = num_points[si, mi, oi] 
218       
219                      # Loop over the time points, back calculating the R2eff values. 
220                      for i in range(num_points_i): 
221                          # Test the two different methods. 
222                          # The A+/- matrices. 
223                          A_pos_i = A_pos_mat[si, mi, oi, i] 
224                          A_neg_i = A_neg_mat[si, mi, oi, i] 
225       
226                          # The lower dimensional matrix exponential. 
227                          A_pos = np_matrix_exponential(m1_mat[si, mi, oi, i]) 
228                          A_neg = np_matrix_exponential(m2_mat[si, mi, oi, i]) 
229       
230                          # Calculate differences 
231                          diff_A_pos_real = A_pos_i.real - A_pos.real 
232                          diff_A_pos_real_sum = sum(diff_A_pos_real) 
233                          diff_A_pos_imag = A_pos_i.imag - A_pos.imag 
234                          diff_A_pos_imag_sum = sum(diff_A_pos_imag) 
235   
236                          diff_A_neg_real = A_neg_i.real - A_neg.real 
237                          diff_A_neg_real_sum = sum(diff_A_neg_real) 
238                          diff_A_neg_imag = A_neg_i.imag - A_neg.imag 
239                          diff_A_neg_imag_sum = sum(diff_A_neg_imag) 
240   
241                          # Test that the sum difference is zero.                                         
242                          self.assertAlmostEqual(diff_A_pos_real_sum, 0.0) 
243                          self.assertAlmostEqual(diff_A_pos_imag_sum, 0.0) 
244                          self.assertAlmostEqual(diff_A_neg_real_sum, 0.0) 
245                          self.assertAlmostEqual(diff_A_neg_imag_sum, 0.0) 
246