Package test_suite :: Package unit_tests :: Package _lib :: Package _dispersion :: Module test_cr72_full_cluster_three_fields
[hide private]
[frames] | no frames]

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

  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 copy import deepcopy 
 24  from numpy import array, arange, int32, float64, max, ones, pi 
 25  from unittest import TestCase 
 26   
 27  # relax module imports. 
 28  from lib.dispersion.cr72 import r2eff_CR72 
 29   
 30  """The 1H gyromagnetic ratio.""" 
 31  g1H = 26.7522212 * 1e7 
 32   
 33  """The 15N gyromagnetic ratio.""" 
 34  g15N = -2.7126 * 1e7 
 35   
 36   
37 -class Test_cr72_full_cluster_three_fields(TestCase):
38 """Unit tests for the lib.dispersion.cr72 relax module.""" 39
40 - def setUp(self):
41 """Set up for all unit tests.""" 42 43 # Default parameter values. 44 self.r2 = None 45 self.r2a = 5.0 46 self.r2b = 10.0 47 self.dw = 3.0 48 self.pA = 0.9 49 self.kex = 1000.0 50 self.spins_params = ['r2a', 'r2b', 'dw', 'pA', 'kex'] 51 52 # Define parameters 53 self.model = "CR72 full" 54 self.num_spins = 2 55 #self.fields = array([800. * 1E6]) 56 #self.fields = array([600. * 1E6, 800. * 1E6]) 57 self.fields = array([600. * 1E6, 800. * 1E6, 900. * 1E6]) 58 self.exp_type = ['SQ CPMG'] 59 self.offset = [0] 60 61 # Required data structures. 62 self.relax_times = self.fields / (100 * 100. *1E6 ) 63 self.ncycs = [] 64 self.points = [] 65 self.value = [] 66 self.error = [] 67 self.num_disp_points = [] 68 69 for i in range(len(self.fields)): 70 ncyc = arange(2, 1000. * self.relax_times[i], 4) 71 #ncyc = arange(2, 42, 2) 72 self.ncycs.append(ncyc) 73 #print("sfrq: ", self.fields[i], "number of cpmg frq", len(ncyc), ncyc) 74 75 cpmg_point = ncyc / self.relax_times[i] 76 77 nr_cpmg_points = len(cpmg_point) 78 self.num_disp_points.append(nr_cpmg_points) 79 80 self.points.append(list(cpmg_point)) 81 self.value.append([2.0]*len(cpmg_point)) 82 self.error.append([1.0]*len(cpmg_point))
83 84
85 - def calc_r2eff(self):
86 """Calculate and check the R2eff values.""" 87 88 # Assemble param vector. 89 self.params = self.assemble_param_vector(r2=self.r2, r2a=self.r2a, r2b=self.r2b, dw=self.dw, pA=self.pA, kex=self.kex, spins_params=self.spins_params) 90 91 # Make nested list arrays of data. And return them. 92 values, errors, cpmg_frqs, missing, frqs, exp_types, relax_times, offsets = self.return_r2eff_arrays() 93 94 # Unpack the parameter values. 95 # Initialise the post spin parameter indices. 96 end_index = [] 97 # The spin and frequency dependent R2 parameters. 98 end_index.append(len(self.exp_type) * self.num_spins * len(self.fields)) 99 if self.model in ["CR72 full"]: 100 end_index.append(2 * len(self.exp_type) * self.num_spins * len(self.fields)) 101 # The spin and dependent parameters (phi_ex, dw, padw2). 102 end_index.append(end_index[-1] + self.num_spins) 103 104 # Unpack the parameter values. 105 R20 = self.params[:end_index[1]].reshape(self.num_spins*2, len(self.fields)) 106 R20A = R20[::2].flatten() 107 R20B = R20[1::2].flatten() 108 dw = self.params[end_index[1]:end_index[2]] 109 pA = self.params[end_index[2]] 110 kex = self.params[end_index[2]+1] 111 112 # Copy value structure 113 self.back_calc = deepcopy(values) 114 115 # Setup special numpy array structures, for higher dimensional computation. 116 # Get the shape of back_calc structure. 117 back_calc_shape = [1, self.num_spins, len(self.fields), 1] 118 119 # Find which frequency has the maximum number of disp points. 120 # To let the numpy array operate well together, the broadcast size has to be equal for all shapes. 121 self.max_num_disp_points = max(self.num_disp_points) 122 123 # Create numpy arrays to pass to the lib function. 124 # All numpy arrays have to have same shape to allow to multiply together. 125 # The dimensions should be [ei][si][mi][oi][di]. [Experiment][spins][spec. frq][offset][disp points]. 126 # The number of disp point can change per spectrometer, so we make the maximum size. 127 self.R20A_a = ones(back_calc_shape + [self.max_num_disp_points]) 128 self.R20B_a = ones(back_calc_shape + [self.max_num_disp_points]) 129 self.dw_frq_a = ones(back_calc_shape + [self.max_num_disp_points]) 130 self.cpmg_frqs_a = ones(back_calc_shape + [self.max_num_disp_points]) 131 self.num_disp_points_a = ones(back_calc_shape + [self.max_num_disp_points]) 132 self.back_calc_a = ones(back_calc_shape + [self.max_num_disp_points]) 133 134 # Loop over the spins. 135 for si in range(self.num_spins): 136 # Loop over the spectrometer frequencies. 137 for mi in range(len(self.fields)): 138 # Extract number of dispersion points. 139 num_disp_points = self.num_disp_points[mi] 140 141 # Extract cpmg_frqs and num_disp_points from lists. 142 self.cpmg_frqs_a[0][si][mi][0][:num_disp_points] = cpmg_frqs[0][mi][0] 143 self.num_disp_points_a[0][si][mi][0][:num_disp_points] = self.num_disp_points[mi] 144 145 # Now calculate. 146 147 # Loop over the spins. 148 for si in range(self.num_spins): 149 # Loop over the spectrometer frequencies. 150 for mi in range(len(self.fields)): 151 # Extract number of dispersion points. 152 num_disp_points = len(cpmg_frqs[0][mi][0]) 153 154 # The R20 index. 155 r20_index = mi + si*len(self.fields) 156 157 # Store r20a and r20b values per disp point. 158 self.R20A_a[0][si][mi][0] = array( [R20A[r20_index]] * self.max_num_disp_points, float64) 159 self.R20B_a[0][si][mi][0] = array( [R20B[r20_index]] * self.max_num_disp_points, float64) 160 161 # Convert dw from ppm to rad/s. 162 dw_frq = dw[si] * frqs[0][si][mi] 163 164 # Store dw_frq per disp point. 165 self.dw_frq_a[0][si][mi][0] = array( [dw_frq] * self.max_num_disp_points, float64) 166 167 ## Back calculate the R2eff values. 168 r2eff_CR72(r20a_orig=self.R20A_a, r20b_orig=self.R20B_a, dw_orig=self.dw_frq_a, r20a=self.R20A_a, r20b=self.R20B_a, pA=pA, dw=self.dw_frq_a, kex=kex, cpmg_frqs=self.cpmg_frqs_a, back_calc=self.back_calc_a) 169 170 # Now return the values back to the structure of self.back_calc object. 171 ## For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 172 # Loop over the spins. 173 for si in range(self.num_spins): 174 # Loop over the spectrometer frequencies. 175 for mi in range(len(self.fields)): 176 # Extract number of dispersion points. 177 num_disp_points = self.num_disp_points[mi] 178 179 # Extract the value 180 self.back_calc[0][si][mi][0][:] = self.back_calc_a[0][si][mi][0][:num_disp_points] 181 182 # Check values. 183 for di in range(num_disp_points): 184 self.assertAlmostEqual(self.back_calc[0][si][mi][0][di], self.R20A_a[0][si][mi][0][di])
185 186
187 - def return_r2eff_arrays(self):
188 """Return numpy arrays of the R2eff/R1rho values and errors. 189 190 @return: The numpy array structures of the R2eff/R1rho values, errors, missing data, and corresponding Larmor frequencies. For each structure, the first dimension corresponds to the experiment types, the second the spins of a spin block, the third to the spectrometer field strength, and the fourth is the dispersion points. For the Larmor frequency structure, the fourth dimension is omitted. For R1rho-type data, an offset dimension is inserted between the spectrometer field strength and the dispersion points. 191 @rtype: lists of numpy float arrays, lists of numpy float arrays, lists of numpy float arrays, numpy rank-2 int array 192 """ 193 194 # Initialise the data structures for the target function. 195 exp_types = [] 196 values = [] 197 errors = [] 198 missing = [] 199 frqs = [] 200 frqs_H = [] 201 relax_times = [] 202 offsets = [] 203 for ei in range(len(self.exp_type)): 204 values.append([]) 205 errors.append([]) 206 missing.append([]) 207 frqs.append([]) 208 frqs_H.append([]) 209 relax_times.append([]) 210 offsets.append([]) 211 for si in range(self.num_spins): 212 values[ei].append([]) 213 errors[ei].append([]) 214 missing[ei].append([]) 215 frqs[ei].append([]) 216 frqs_H[ei].append([]) 217 offsets[ei].append([]) 218 for mi in range(len(self.fields)): 219 values[ei][si].append([]) 220 errors[ei][si].append([]) 221 missing[ei][si].append([]) 222 frqs[ei][si].append(0.0) 223 frqs_H[ei][si].append(0.0) 224 offsets[ei][si].append([]) 225 for oi in range(len(self.offset)): 226 values[ei][si][mi].append([]) 227 errors[ei][si][mi].append([]) 228 missing[ei][si][mi].append([]) 229 offsets[ei][si][mi].append([]) 230 for mi in range(len(self.fields)): 231 relax_times[ei].append(None) 232 233 cpmg_frqs = [] 234 for ei in range(len(self.exp_type)): 235 cpmg_frqs.append([]) 236 for mi in range(len(self.fields)): 237 cpmg_frqs[ei].append([]) 238 for oi in range(len(self.offset)): 239 #cpmg_frqs[ei][mi].append(self.points) 240 cpmg_frqs[ei][mi].append([]) 241 242 243 # Pack the R2eff/R1rho data. 244 si = 0 245 for spin_index in range(self.num_spins): 246 data_flag = True 247 248 for ei in range(len(self.exp_type)): 249 exp_type = self.exp_type[ei] 250 # Add the experiment type. 251 if exp_type not in exp_types: 252 exp_types.append(exp_type) 253 254 for mi in range(len(self.fields)): 255 # Get the frq. 256 frq = self.fields[mi] 257 258 # The Larmor frequency for this spin (and that of an attached proton for the MMQ models) and field strength (in MHz*2pi to speed up the ppm to rad/s conversion). 259 frqs[ei][si][mi] = 2.0 * pi * frq / g1H * g15N * 1e-6 260 261 # Get the cpmg frq. 262 cpmg_frqs[ei][mi][oi] = self.points[mi] 263 264 for oi in range(len(self.offset)): 265 for di in range(len(self.points[mi])): 266 267 missing[ei][si][mi][oi].append(0) 268 269 # Values 270 values[ei][si][mi][oi].append(self.value[mi][di]) 271 # The errors. 272 errors[ei][si][mi][oi].append(self.error[mi][di]) 273 274 # The relaxation times. 275 # Found. 276 relax_time = self.relax_times[mi] 277 278 # Store the time. 279 relax_times[ei][mi] = relax_time 280 281 # Increment the spin index. 282 si += 1 283 284 # Convert to numpy arrays. 285 relax_times = array(relax_times, float64) 286 for ei in range(len(self.exp_type)): 287 for si in range(self.num_spins): 288 for mi in range(len(self.fields)): 289 for oi in range(len(self.offset)): 290 values[ei][si][mi][oi] = array(values[ei][si][mi][oi], float64) 291 errors[ei][si][mi][oi] = array(errors[ei][si][mi][oi], float64) 292 missing[ei][si][mi][oi] = array(missing[ei][si][mi][oi], int32) 293 294 # Return the structures. 295 return values, errors, cpmg_frqs, missing, frqs, exp_types, relax_times, offsets
296 297
298 - def assemble_param_vector(self, r2=None, r2a=None, r2b=None, dw=None, pA=None, kex=None, spins_params=None):
299 """Assemble the dispersion relaxation dispersion curve fitting parameter vector. 300 301 @keyword r2: The transversal relaxation rate. 302 @type r2: float 303 @keyword r2a: The transversal relaxation rate for state A in the absence of exchange. 304 @type r2a: float 305 @keyword r2b: The transversal relaxation rate for state B in the absence of exchange. 306 @type r2b: float 307 @keyword dw: The chemical exchange difference between states A and B in ppm. 308 @type dw: float 309 @keyword pA: The population of state A. 310 @type pA: float 311 @keyword kex: The rate of exchange. 312 @type kex: float 313 @keyword spins_params: List of parameter strings used in dispersion model. 314 @type spins_params: array of strings 315 @return: An array of the parameter values of the dispersion relaxation model. 316 @rtype: numpy float array 317 """ 318 319 # Initialise. 320 param_vector = [] 321 322 # Loop over the parameters of the cluster. 323 for param_name, spin_index, mi in self.loop_parameters(spins_params=spins_params): 324 if param_name == 'r2': 325 value = r2 326 value = value + mi + spin_index*0.1 327 elif param_name == 'r2a': 328 value = r2a 329 value = value + mi+ spin_index*0.1 330 elif param_name == 'r2b': 331 value = r2b 332 value = value + mi + spin_index*0.1 333 elif param_name == 'dw': 334 value = dw 335 elif param_name == 'pA': 336 value = pA 337 elif param_name == 'kex': 338 value = kex 339 340 # Add to the vector. 341 param_vector.append(value) 342 343 # Return a numpy array. 344 return array(param_vector, float64)
345 346
347 - def loop_parameters(self, spins_params=None):
348 """Generator function for looping of the parameters of the cluster. 349 350 @keyword spins_params: List of parameter strings used in dispersion model. 351 @type spins_params: array of strings 352 @return: The parameter name. 353 @rtype: str 354 """ 355 356 # Loop over the parameters of the cluster. 357 # First the R2 parameters (one per spin per field strength). 358 for spin_index in range(self.num_spins): 359 360 # The R2 parameter. 361 if 'r2' in spins_params: 362 for ei in range(len(self.exp_type)): 363 for mi in range(len(self.fields)): 364 yield 'r2', spin_index, mi 365 366 # The R2A parameter. 367 if 'r2a' in spins_params: 368 for ei in range(len(self.exp_type)): 369 for mi in range(len(self.fields)): 370 yield 'r2a', spin_index, mi 371 372 373 # The R2B parameter. 374 if 'r2b' in spins_params: 375 for ei in range(len(self.exp_type)): 376 for mi in range(len(self.fields)): 377 yield 'r2b', spin_index, mi 378 379 380 # Then the chemical shift difference parameters 'phi_ex', 'phi_ex_B', 'phi_ex_C', 'padw2', 'dw', 'dw_AB', 'dw_BC', 'dw_AB' (one per spin). 381 for spin_index in range(self.num_spins): 382 383 if 'dw' in spins_params: 384 yield 'dw', spin_index, 0 385 386 # All other parameters (one per spin cluster). 387 for param in spins_params: 388 if not param in ['r2', 'r2a', 'r2b', 'phi_ex', 'phi_ex_B', 'phi_ex_C', 'padw2', 'dw', 'dw_AB', 'dw_BC', 'dw_AB', 'dwH', 'dwH_AB', 'dwH_BC', 'dwH_AB']: 389 if param == 'pA': 390 yield 'pA', 0, 0 391 elif param == 'kex': 392 yield 'kex', 0, 0
393 394
396 """Test the r2eff_cr72() function for no exchange when dw = 0.0.""" 397 398 # Parameter reset. 399 self.dw = 0.0 400 401 # Calculate and check the R2eff values. 402 self.calc_r2eff()
403 404
406 """Test the r2eff_cr72() function for no exchange when pA = 1.0.""" 407 408 # Parameter reset. 409 self.pA = 1.0 410 411 # Calculate and check the R2eff values. 412 self.calc_r2eff()
413 414
416 """Test the r2eff_cr72() function for no exchange when kex = 0.0.""" 417 418 # Parameter reset. 419 self.kex = 0.0 420 421 # Calculate and check the R2eff values. 422 self.calc_r2eff()
423 424
426 """Test the r2eff_cr72() function for no exchange when dw = 0.0 and pA = 1.0.""" 427 428 # Parameter reset. 429 self.pA = 1.0 430 self.dw = 0.0 431 432 # Calculate and check the R2eff values. 433 self.calc_r2eff()
434 435
437 """Test the r2eff_cr72() function for no exchange when dw = 0.0 and kex = 0.0.""" 438 439 # Parameter reset. 440 self.dw = 0.0 441 self.kex = 0.0 442 443 # Calculate and check the R2eff values. 444 self.calc_r2eff()
445 446
448 """Test the r2eff_cr72() function for no exchange when pA = 1.0 and kex = 0.0.""" 449 450 # Parameter reset. 451 self.pA = 1.0 452 self.kex = 0.0 453 454 # Calculate and check the R2eff values. 455 self.calc_r2eff()
456 457
459 """Test the r2eff_cr72() function for no exchange when dw = 0.0, pA = 1.0, and kex = 0.0.""" 460 461 # Parameter reset. 462 self.dw = 0.0 463 self.pA = 1.0 464 self.kex = 0.0 465 466 # Calculate and check the R2eff values. 467 self.calc_r2eff()
468 469
471 """Test the r2eff_cr72() function for no exchange when kex = 1e7.""" 472 473 # Parameter reset. 474 self.kex = 1e7 475 476 # Calculate and check the R2eff values. 477 self.calc_r2eff()
478