Trees       Indices       Help   
relax
Package test_suite :: Package unit_tests :: Package _lib :: Package _dispersion :: Module test_matrix_exponential
[hide private]
[frames] | no frames]

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

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2014 Troels E. Linnet                                         # 
  4  #                                                                             # 
  5  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  6  #                                                                             # 
  7  # This program is free software: you can redistribute it and/or modify        # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation, either version 3 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # This program is distributed in the hope that it will be useful,             # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 19  #                                                                             # 
 20  ############################################################################### 
 21   
 22  # Python module imports. 
 23  from os import sep 
 24  from numpy import complex64, load, sum 
 25  from unittest import TestCase 
 26   
 27  # relax module imports. 
 28  from lib.dispersion.ns_cpmg_2site_3d import rcpmg_3d_rankN 
 29  from lib.dispersion.ns_mmq_2site import rmmq_2site_rankN 
 30  from lib.linear_algebra.matrix_exponential import matrix_exponential as np_matrix_exponential 
 31  from lib.dispersion.matrix_exponential import matrix_exponential 
 32  from status import Status; status = Status() 
 33   
 34   

 35 -class Test_matrix_exponential(TestCase): 

 36      """Unit tests for the lib.dispersion.matrix_exponential relax module.""" 
 37   

 38 -    def setUp(self): 

 39          """Set up for all unit tests.""" 
 40   
 41          # Default parameter values. 
 42          self.data = status.install_path+sep+'test_suite'+sep+'shared_data'+sep+'dispersion'+sep+'unit_tests'+sep+'lib'+sep+'dispersion'+sep+'matrix_exponential' 

 43   

 44 -    def return_data_ns_cpmg_2site_3d(self, fname): 

 45          """Return the data structures from the data path.""" 
 46   
 47          r180x = load(fname+"_r180x"+".npy") 
 48          M0 = load(fname+"_M0"+".npy") 
 49          r10a = load(fname+"_r10a"+".npy") 
 50          r10b = load(fname+"_r10b"+".npy") 
 51          r20a = load(fname+"_r20a"+".npy") 
 52          r20b = load(fname+"_r20b"+".npy") 
 53          pA = load(fname+"_pA"+".npy") 
 54          dw = load(fname+"_dw"+".npy") 
 55          dw_orig = load(fname+"_dw_orig"+".npy") 
 56          kex = load(fname+"_kex"+".npy") 
 57          inv_tcpmg = load(fname+"_inv_tcpmg"+".npy") 
 58          tcp = load(fname+"_tcp"+".npy") 
 59          num_points = load(fname+"_num_points"+".npy") 
 60          power = load(fname+"_power"+".npy") 
 61          back_calc = load(fname+"_back_calc"+".npy") 
 62   
 63          # Once off parameter conversions. 
 64          pB = 1.0 - pA 
 65          k_BA = pA * kex 
 66          k_AB = pB * kex 
 67   
 68          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 

 69   
 70   

 71 -    def return_data_mmq_2site(self, fname): 

 72          """Return the data structures from the data path.""" 
 73   
 74          M0 = load(fname+"_M0"+".npy") 
 75          r20a = load(fname+"_r20a"+".npy") 
 76          r20b = load(fname+"_r20b"+".npy") 
 77          pA = load(fname+"_pA"+".npy") 
 78          dw = load(fname+"_dw"+".npy") 
 79          dwH = load(fname+"_dwH"+".npy") 
 80          kex = load(fname+"_kex"+".npy") 
 81          inv_tcpmg = load(fname+"_inv_tcpmg"+".npy") 
 82          tcp = load(fname+"_tcp"+".npy") 
 83          num_points = load(fname+"_num_points"+".npy") 
 84          power = load(fname+"_power"+".npy") 
 85          back_calc = load(fname+"_back_calc"+".npy") 
 86   
 87          # Once off parameter conversions. 
 88          pB = 1.0 - pA 
 89          k_BA = pA * kex 
 90          k_AB = pB * kex 
 91   
 92          return M0, r20a, r20b, pA, dw, dwH, kex, inv_tcpmg, tcp, num_points, power, back_calc, pB, k_BA, k_AB 

 93   
 94   

 96          """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.""" 
 97   
 98          fname = self.data + sep+ "test_hansen_cpmg_data_to_ns_cpmg_2site_3D" 
 99          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) 
100   
101          # Extract the total numbers of experiments, number of spins, number of magnetic field strength, number of offsets, maximum number of dispersion point. 
102          NE, NS, NM, NO, ND = back_calc.shape 
103   
104          # The matrix R that contains all the contributions to the evolution, i.e. relaxation, exchange and chemical shift evolution. 
105          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) 
106       
107          # This matrix is a propagator that will evolve the magnetization with the matrix R for a delay tcp. 
108          Rexpo_mat = matrix_exponential(R_mat) 
109       
110          # Loop over the spins 
111          for si in range(NS): 
112              # Loop over the spectrometer frequencies. 
113              for mi in range(NM): 
114                  # Extract number of points. 
115                  num_points_si_mi = int(num_points[0, si, mi, 0]) 
116       
117                  # Loop over the time points, back calculating the R2eff values. 
118                  for di in range(num_points_si_mi): 
119                      # Test the two different methods. 
120                      R_mat_i = R_mat[0, si, mi, 0, di] 
121    
122                      # The lower dimensional matrix exponential. 
123                      Rexpo = np_matrix_exponential(R_mat_i) 
124       
125                      # The higher dimensional matrix exponential. 
126                      Rexpo_mat_i = Rexpo_mat[0, si, mi, 0, di] 
127   
128                      diff = Rexpo - Rexpo_mat_i 
129                      diff_sum = sum(diff) 
130   
131                      # Test that the sum difference is zero.                                         
132                      self.assertAlmostEqual(diff_sum, 0.0) 

133   
134   

136          """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. 
137          This test does the matrix exponential in complex64.""" 
138   
139          fname = self.data + sep+ "test_korzhnev_2005_15n_dq_data" 
140          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) 
141   
142          # Extract the total numbers of experiments, number of spins, number of magnetic field strength, number of offsets, maximum number of dispersion point. 
143          NS, NM, NO = num_points.shape 
144   
145          # Populate the m1 and m2 matrices (only once per function call for speed). 
146          m1_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 
147          m2_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=-dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 
148       
149          # The A+/- matrices. 
150          A_pos_mat = matrix_exponential(m1_mat, dtype=complex64) 
151          A_neg_mat = matrix_exponential(m2_mat, dtype=complex64) 
152       
153          # Loop over spins. 
154          for si in range(NS): 
155              # Loop over the spectrometer frequencies. 
156              for mi in range(NM): 
157                  # Loop over offsets: 
158                  for oi in range(NO): 
159                      # Extract number of points. 
160                      num_points_i = num_points[si, mi, oi] 
161       
162                      # Loop over the time points, back calculating the R2eff values. 
163                      for i in range(num_points_i): 
164                          # Test the two different methods. 
165                          # The A+/- matrices. 
166                          A_pos_i = A_pos_mat[si, mi, oi, i] 
167                          A_neg_i = A_neg_mat[si, mi, oi, i] 
168       
169                          # The lower dimensional matrix exponential. 
170                          A_pos = np_matrix_exponential(m1_mat[si, mi, oi, i]) 
171                          A_neg = np_matrix_exponential(m2_mat[si, mi, oi, i]) 
172       
173                          # Calculate differences 
174                          diff_A_pos_real = A_pos_i.real - A_pos.real 
175                          diff_A_pos_real_sum = sum(diff_A_pos_real) 
176                          diff_A_pos_imag = A_pos_i.imag - A_pos.imag 
177                          diff_A_pos_imag_sum = sum(diff_A_pos_imag) 
178   
179                          diff_A_neg_real = A_neg_i.real - A_neg.real 
180                          diff_A_neg_real_sum = sum(diff_A_neg_real) 
181                          diff_A_neg_imag = A_neg_i.imag - A_neg.imag 
182                          diff_A_neg_imag_sum = sum(diff_A_neg_imag) 
183   
184                          # Test that the sum difference is zero.                                         
185                          self.assertAlmostEqual(diff_A_pos_real_sum, 0.0, 6) 
186                          self.assertAlmostEqual(diff_A_pos_imag_sum, 0.0, 6) 
187                          self.assertAlmostEqual(diff_A_neg_real_sum, 0.0, 6) 
188                          self.assertAlmostEqual(diff_A_neg_imag_sum, 0.0, 6) 

189   
190   

192          """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. 
193          This test does the matrix exponential in complex128.""" 
194   
195          fname = self.data + sep+ "test_korzhnev_2005_15n_dq_data" 
196          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) 
197   
198          # Extract the total numbers of experiments, number of spins, number of magnetic field strength, number of offsets, maximum number of dispersion point. 
199          NS, NM, NO = num_points.shape 
200   
201          # Populate the m1 and m2 matrices (only once per function call for speed). 
202          m1_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 
203          m2_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=-dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 
204       
205          # The A+/- matrices. 
206          A_pos_mat = matrix_exponential(m1_mat) 
207          A_neg_mat = matrix_exponential(m2_mat) 
208       
209          # Loop over spins. 
210          for si in range(NS): 
211              # Loop over the spectrometer frequencies. 
212              for mi in range(NM): 
213                  # Loop over offsets: 
214                  for oi in range(NO): 
215                      # Extract number of points. 
216                      num_points_i = num_points[si, mi, oi] 
217       
218                      # Loop over the time points, back calculating the R2eff values. 
219                      for i in range(num_points_i): 
220                          # Test the two different methods. 
221                          # The A+/- matrices. 
222                          A_pos_i = A_pos_mat[si, mi, oi, i] 
223                          A_neg_i = A_neg_mat[si, mi, oi, i] 
224       
225                          # The lower dimensional matrix exponential. 
226                          A_pos = np_matrix_exponential(m1_mat[si, mi, oi, i]) 
227                          A_neg = np_matrix_exponential(m2_mat[si, mi, oi, i]) 
228       
229                          # Calculate differences 
230                          diff_A_pos_real = A_pos_i.real - A_pos.real 
231                          diff_A_pos_real_sum = sum(diff_A_pos_real) 
232                          diff_A_pos_imag = A_pos_i.imag - A_pos.imag 
233                          diff_A_pos_imag_sum = sum(diff_A_pos_imag) 
234   
235                          diff_A_neg_real = A_neg_i.real - A_neg.real 
236                          diff_A_neg_real_sum = sum(diff_A_neg_real) 
237                          diff_A_neg_imag = A_neg_i.imag - A_neg.imag 
238                          diff_A_neg_imag_sum = sum(diff_A_neg_imag) 
239   
240                          # Test that the sum difference is zero.                                         
241                          self.assertAlmostEqual(diff_A_pos_real_sum, 0.0) 
242                          self.assertAlmostEqual(diff_A_pos_imag_sum, 0.0) 
243                          self.assertAlmostEqual(diff_A_neg_real_sum, 0.0) 
244                          self.assertAlmostEqual(diff_A_neg_imag_sum, 0.0) 

245

   Trees       Indices       Help   
relax
Generated by Epydoc 3.0.1 on Fri Jun 14 11:30:14 2019 http://epydoc.sourceforge.net