test_suite.unit_tests._lib._dispersion.test_cr72_full_cluster_one

39 """Unit tests for the lib.dispersion.cr72 relax module.""" 40

41 - def setUp(self):

42 """Set up for all unit tests.""" 43 44 # Default parameter values. 45 self.r2 = None 46 self.r2a = 5.0 47 self.r2b = 10.0 48 self.dw = 3.0 49 self.pA = 0.9 50 self.kex = 1000.0 51 self.spins_params = ['r2a', 'r2b', 'dw', 'pA', 'kex'] 52 53 # Define parameters 54 self.model = "CR72 full" 55 self.num_spins = 2 56 self.fields = array([800. * 1E6]) 57 #self.fields = array([600. * 1E6, 800. * 1E6]) 58 #self.fields = array([600. * 1E6, 800. * 1E6, 900. * 1E6]) 59 self.exp_type = ['SQ CPMG'] 60 self.offset = [0] 61 62 # Required data structures. 63 self.relax_times = self.fields / (100 * 100. *1E6 ) 64 self.ncycs = [] 65 self.points = [] 66 self.value = [] 67 self.error = [] 68 self.num_disp_points = [] 69 70 for i in range(len(self.fields)): 71 ncyc = arange(2, 1000. * self.relax_times[i], 4) 72 #ncyc = arange(2, 42, 2) 73 self.ncycs.append(ncyc) 74 #print("sfrq: ", self.fields[i], "number of cpmg frq", len(ncyc), ncyc) 75 76 cpmg_point = ncyc / self.relax_times[i] 77 78 nr_cpmg_points = len(cpmg_point) 79 self.num_disp_points.append(nr_cpmg_points) 80 81 self.points.append(list(cpmg_point)) 82 self.value.append([2.0]*len(cpmg_point)) 83 self.error.append([1.0]*len(cpmg_point))

84 85

86 - def calc_r2eff(self):

87 """Calculate and check the R2eff values.""" 88 89 # Assemble param vector. 90 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) 91 92 # Make nested list arrays of data. And return them. 93 values, errors, cpmg_frqs, missing, frqs, exp_types, relax_times, offsets = self.return_r2eff_arrays() 94 95 # Unpack the parameter values. 96 # Initialise the post spin parameter indices. 97 end_index = [] 98 # The spin and frequency dependent R2 parameters. 99 end_index.append(len(self.exp_type) * self.num_spins * len(self.fields)) 100 if self.model in ["CR72 full"]: 101 end_index.append(2 * len(self.exp_type) * self.num_spins * len(self.fields)) 102 # The spin and dependent parameters (phi_ex, dw, padw2). 103 end_index.append(end_index[-1] + self.num_spins) 104 105 # Unpack the parameter values. 106 R20 = self.params[:end_index[1]].reshape(self.num_spins*2, len(self.fields)) 107 R20A = R20[::2].flatten() 108 R20B = R20[1::2].flatten() 109 dw = self.params[end_index[1]:end_index[2]] 110 pA = self.params[end_index[2]] 111 kex = self.params[end_index[2]+1] 112 113 # Copy value structure 114 self.back_calc = deepcopy(values) 115 116 # Setup special numpy array structures, for higher dimensional computation. 117 # Get the shape of back_calc structure. 118 back_calc_shape = list( asarray(self.back_calc).shape )[:4] 119 120 # Find which frequency has the maximum number of disp points. 121 # To let the numpy array operate well together, the broadcast size has to be equal for all shapes. 122 self.max_num_disp_points = max(self.num_disp_points) 123 124 # Create numpy arrays to pass to the lib function. 125 # All numpy arrays have to have same shape to allow to multiply together. 126 # The dimensions should be [ei][si][mi][oi][di]. [Experiment][spins][spec. frq][offset][disp points]. 127 # The number of disp point can change per spectrometer, so we make the maximum size. 128 self.R20A_a = ones(back_calc_shape + [self.max_num_disp_points]) 129 self.R20B_a = ones(back_calc_shape + [self.max_num_disp_points]) 130 self.dw_frq_a = ones(back_calc_shape + [self.max_num_disp_points]) 131 self.cpmg_frqs_a = ones(back_calc_shape + [self.max_num_disp_points]) 132 self.num_disp_points_a = ones(back_calc_shape + [self.max_num_disp_points]) 133 self.back_calc_a = ones(back_calc_shape + [self.max_num_disp_points]) 134 135 # Loop over the spins. 136 for si in range(self.num_spins): 137 # Loop over the spectrometer frequencies. 138 for mi in range(len(self.fields)): 139 # Extract number of dispersion points. 140 num_disp_points = self.num_disp_points[mi] 141 142 # Extract cpmg_frqs and num_disp_points from lists. 143 self.cpmg_frqs_a[0][si][mi][0][:num_disp_points] = cpmg_frqs[0][mi][0] 144 self.num_disp_points_a[0][si][mi][0][:num_disp_points] = self.num_disp_points[mi] 145 146 # Now calculate. 147 148 # Loop over the spins. 149 for si in range(self.num_spins): 150 # Loop over the spectrometer frequencies. 151 for mi in range(len(self.fields)): 152 # Extract number of dispersion points. 153 num_disp_points = len(cpmg_frqs[0][mi][0]) 154 155 # The R20 index. 156 r20_index = mi + si*len(self.fields) 157 158 # Store r20a and r20b values per disp point. 159 self.R20A_a[0][si][mi][0] = array( [R20A[r20_index]] * self.max_num_disp_points, float64) 160 self.R20B_a[0][si][mi][0] = array( [R20B[r20_index]] * self.max_num_disp_points, float64) 161 162 # Convert dw from ppm to rad/s. 163 dw_frq = dw[si] * frqs[0][si][mi] 164 165 # Store dw_frq per disp point. 166 self.dw_frq_a[0][si][mi][0] = array( [dw_frq] * self.max_num_disp_points, float64) 167 168 ## Back calculate the R2eff values. 169 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) 170 171 # Now return the values back to the structure of self.back_calc object. 172 ## 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. 173 # Loop over the spins. 174 for si in range(self.num_spins): 175 # Loop over the spectrometer frequencies. 176 for mi in range(len(self.fields)): 177 # Extract number of dispersion points. 178 num_disp_points = self.num_disp_points[mi] 179 180 # Extract the value 181 self.back_calc[0][si][mi][0][:] = self.back_calc_a[0][si][mi][0][:num_disp_points] 182 183 # Check values. 184 for di in range(num_disp_points): 185 self.assertAlmostEqual(self.back_calc[0][si][mi][0][di], self.R20A_a[0][si][mi][0][di])

186 187

188 - def return_r2eff_arrays(self):

189 """Return numpy arrays of the R2eff/R1rho values and errors. 190 191 @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. 192 @rtype: lists of numpy float arrays, lists of numpy float arrays, lists of numpy float arrays, numpy rank-2 int array 193 """ 194 195 # Initialise the data structures for the target function. 196 exp_types = [] 197 values = [] 198 errors = [] 199 missing = [] 200 frqs = [] 201 frqs_H = [] 202 relax_times = [] 203 offsets = [] 204 for ei in range(len(self.exp_type)): 205 values.append([]) 206 errors.append([]) 207 missing.append([]) 208 frqs.append([]) 209 frqs_H.append([]) 210 relax_times.append([]) 211 offsets.append([]) 212 for si in range(self.num_spins): 213 values[ei].append([]) 214 errors[ei].append([]) 215 missing[ei].append([]) 216 frqs[ei].append([]) 217 frqs_H[ei].append([]) 218 offsets[ei].append([]) 219 for mi in range(len(self.fields)): 220 values[ei][si].append([]) 221 errors[ei][si].append([]) 222 missing[ei][si].append([]) 223 frqs[ei][si].append(0.0) 224 frqs_H[ei][si].append(0.0) 225 offsets[ei][si].append([]) 226 for oi in range(len(self.offset)): 227 values[ei][si][mi].append([]) 228 errors[ei][si][mi].append([]) 229 missing[ei][si][mi].append([]) 230 offsets[ei][si][mi].append([]) 231 for mi in range(len(self.fields)): 232 relax_times[ei].append(None) 233 234 cpmg_frqs = [] 235 for ei in range(len(self.exp_type)): 236 cpmg_frqs.append([]) 237 for mi in range(len(self.fields)): 238 cpmg_frqs[ei].append([]) 239 for oi in range(len(self.offset)): 240 #cpmg_frqs[ei][mi].append(self.points) 241 cpmg_frqs[ei][mi].append([]) 242 243 244 # Pack the R2eff/R1rho data. 245 si = 0 246 for spin_index in range(self.num_spins): 247 data_flag = True 248 249 for ei in range(len(self.exp_type)): 250 exp_type = self.exp_type[ei] 251 # Add the experiment type. 252 if exp_type not in exp_types: 253 exp_types.append(exp_type) 254 255 for mi in range(len(self.fields)): 256 # Get the frq. 257 frq = self.fields[mi] 258 259 # 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). 260 frqs[ei][si][mi] = 2.0 * pi * frq / g1H * g15N * 1e-6 261 262 # Get the cpmg frq. 263 cpmg_frqs[ei][mi][oi] = self.points[mi] 264 265 for oi in range(len(self.offset)): 266 for di in range(len(self.points[mi])): 267 268 missing[ei][si][mi][oi].append(0) 269 270 # Values 271 values[ei][si][mi][oi].append(self.value[mi][di]) 272 # The errors. 273 errors[ei][si][mi][oi].append(self.error[mi][di]) 274 275 # The relaxation times. 276 # Found. 277 relax_time = self.relax_times[mi] 278 279 # Store the time. 280 relax_times[ei][mi] = relax_time 281 282 # Increment the spin index. 283 si += 1 284 285 # Convert to numpy arrays. 286 relax_times = array(relax_times, float64) 287 for ei in range(len(self.exp_type)): 288 for si in range(self.num_spins): 289 for mi in range(len(self.fields)): 290 for oi in range(len(self.offset)): 291 values[ei][si][mi][oi] = array(values[ei][si][mi][oi], float64) 292 errors[ei][si][mi][oi] = array(errors[ei][si][mi][oi], float64) 293 missing[ei][si][mi][oi] = array(missing[ei][si][mi][oi], int32) 294 295 # Return the structures. 296 return values, errors, cpmg_frqs, missing, frqs, exp_types, relax_times, offsets

297 298

299 - def assemble_param_vector(self, r2=None, r2a=None, r2b=None, dw=None, pA=None, kex=None, spins_params=None):

300 """Assemble the dispersion relaxation dispersion curve fitting parameter vector. 301 302 @keyword r2: The transversal relaxation rate. 303 @type r2: float 304 @keyword r2a: The transversal relaxation rate for state A in the absence of exchange. 305 @type r2a: float 306 @keyword r2b: The transversal relaxation rate for state B in the absence of exchange. 307 @type r2b: float 308 @keyword dw: The chemical exchange difference between states A and B in ppm. 309 @type dw: float 310 @keyword pA: The population of state A. 311 @type pA: float 312 @keyword kex: The rate of exchange. 313 @type kex: float 314 @keyword spins_params: List of parameter strings used in dispersion model. 315 @type spins_params: array of strings 316 @return: An array of the parameter values of the dispersion relaxation model. 317 @rtype: numpy float array 318 """ 319 320 # Initialise. 321 param_vector = [] 322 323 # Loop over the parameters of the cluster. 324 for param_name, spin_index, mi in self.loop_parameters(spins_params=spins_params): 325 if param_name == 'r2': 326 value = r2 327 value = value + mi + spin_index*0.1 328 elif param_name == 'r2a': 329 value = r2a 330 value = value + mi+ spin_index*0.1 331 elif param_name == 'r2b': 332 value = r2b 333 value = value + mi + spin_index*0.1 334 elif param_name == 'dw': 335 value = dw 336 elif param_name == 'pA': 337 value = pA 338 elif param_name == 'kex': 339 value = kex 340 341 # Add to the vector. 342 param_vector.append(value) 343 344 # Return a numpy array. 345 return array(param_vector, float64)

346 347

348 - def loop_parameters(self, spins_params=None):

349 """Generator function for looping of the parameters of the cluster. 350 351 @keyword spins_params: List of parameter strings used in dispersion model. 352 @type spins_params: array of strings 353 @return: The parameter name. 354 @rtype: str 355 """ 356 357 # Loop over the parameters of the cluster. 358 # First the R2 parameters (one per spin per field strength). 359 for spin_index in range(self.num_spins): 360 361 # The R2 parameter. 362 if 'r2' in spins_params: 363 for ei in range(len(self.exp_type)): 364 for mi in range(len(self.fields)): 365 yield 'r2', spin_index, mi 366 367 # The R2A parameter. 368 if 'r2a' in spins_params: 369 for ei in range(len(self.exp_type)): 370 for mi in range(len(self.fields)): 371 yield 'r2a', spin_index, mi 372 373 374 # The R2B parameter. 375 if 'r2b' in spins_params: 376 for ei in range(len(self.exp_type)): 377 for mi in range(len(self.fields)): 378 yield 'r2b', spin_index, mi 379 380 381 # 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). 382 for spin_index in range(self.num_spins): 383 384 if 'dw' in spins_params: 385 yield 'dw', spin_index, 0 386 387 # All other parameters (one per spin cluster). 388 for param in spins_params: 389 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']: 390 if param == 'pA': 391 yield 'pA', 0, 0 392 elif param == 'kex': 393 yield 'kex', 0, 0

394 395

396 - def test_cr72_full_cluster_one_field_no_rex1(self):

397 """Test the r2eff_cr72() function for no exchange when dw = 0.0.""" 398 399 # Parameter reset. 400 self.dw = 0.0 401 402 # Calculate and check the R2eff values. 403 self.calc_r2eff()

404 405

406 - def test_cr72_full_cluster_one_field_no_rex2(self):

407 """Test the r2eff_cr72() function for no exchange when pA = 1.0.""" 408 409 # Parameter reset. 410 self.pA = 1.0 411 412 # Calculate and check the R2eff values. 413 self.calc_r2eff()

414 415

416 - def test_cr72_full_cluster_one_field_no_rex3(self):

417 """Test the r2eff_cr72() function for no exchange when kex = 0.0.""" 418 419 # Parameter reset. 420 self.kex = 0.0 421 422 # Calculate and check the R2eff values. 423 self.calc_r2eff()

424 425

426 - def test_cr72_full_cluster_one_field_no_rex4(self):

427 """Test the r2eff_cr72() function for no exchange when dw = 0.0 and pA = 1.0.""" 428 429 # Parameter reset. 430 self.pA = 1.0 431 self.dw = 0.0 432 433 # Calculate and check the R2eff values. 434 self.calc_r2eff()

435 436

437 - def test_cr72_full_cluster_one_field_no_rex5(self):

438 """Test the r2eff_cr72() function for no exchange when dw = 0.0 and kex = 0.0.""" 439 440 # Parameter reset. 441 self.dw = 0.0 442 self.kex = 0.0 443 444 # Calculate and check the R2eff values. 445 self.calc_r2eff()

446 447

448 - def test_cr72_full_cluster_one_field_no_rex6(self):

449 """Test the r2eff_cr72() function for no exchange when pA = 1.0 and kex = 0.0.""" 450 451 # Parameter reset. 452 self.pA = 1.0 453 self.kex = 0.0 454 455 # Calculate and check the R2eff values. 456 self.calc_r2eff()

457 458

459 - def test_cr72_full_cluster_one_field_no_rex7(self):

460 """Test the r2eff_cr72() function for no exchange when dw = 0.0, pA = 1.0, and kex = 0.0.""" 461 462 # Parameter reset. 463 self.dw = 0.0 464 self.pA = 1.0 465 self.kex = 0.0 466 467 # Calculate and check the R2eff values. 468 self.calc_r2eff()

469 470

471 - def test_cr72_full_cluster_one_field_no_rex8(self):

472 """Test the r2eff_cr72() function for no exchange when kex = 1e7.""" 473 474 # Parameter reset. 475 self.kex = 1e7 476 477 # Calculate and check the R2eff values. 478 self.calc_r2eff()

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