test_suite.system_tests.relax

40 """Class for testing various aspects specific to relaxation curve-fitting.""" 41

42 - def __init__(self, methodName='runTest'):

43 """Skip the tests if the C modules are non-functional. 44 45 @keyword methodName: The name of the test. 46 @type methodName: str 47 """ 48 49 # Execute the base class method. 50 super(Relax_fit, self).__init__(methodName) 51 52 # Missing module. 53 if not dep_check.C_module_exp_fn: 54 # Store in the status object. 55 status.skipped_tests.append([methodName, 'Relax curve-fitting C module', self._skip_type])

56 57

58 - def setUp(self):

59 """Set up for all the functional tests.""" 60 61 # Create the data pipe. 62 self.interpreter.pipe.create('mf', 'mf') 63 64 # Create a temporary directory for dumping files. 65 ds.tmpdir = mkdtemp()

66 67

68 - def check_curve_fitting(self):

69 """Check the results of the curve-fitting.""" 70 71 # Data. 72 relax_times = [0.0176, 0.0176, 0.0352, 0.0704, 0.0704, 0.1056, 0.1584, 0.1584, 0.1936, 0.1936] 73 chi2 = [None, None, None, 2.916952651567855, 5.4916923952919632, 16.21182245065274, 4.3591263759462926, 9.8925377583244316, None, None, None, 6.0238341559877782] 74 rx = [None, None, None, 8.0814894819820662, 8.6478971039559642, 9.5710638183013845, 10.716551838066295, 11.143793935455122, None, None, None, 12.82875370075309] 75 i0 = [None, None, None, 1996050.9679875025, 2068490.9458927638, 1611556.5194095275, 1362887.2331948928, 1877670.5623875158, None, None, None, 897044.17382064369] 76 77 # Some checks. 78 self.assertEqual(cdp.curve_type, 'exp') 79 self.assertEqual(cdp.int_method, ds.int_type) 80 self.assertEqual(len(cdp.relax_times), 10) 81 cdp_relax_times = sorted(cdp.relax_times.values()) 82 for i in range(10): 83 self.assertEqual(cdp_relax_times[i], relax_times[i]) 84 85 # Check the errors. 86 for key in cdp.sigma_I: 87 self.assertAlmostEqual(cdp.sigma_I[key], 10578.039482421433, 6) 88 self.assertAlmostEqual(cdp.var_I[key], 111894919.29166669) 89 90 # Spin data check. 91 i = 0 92 for spin in spin_loop(): 93 # No data present. 94 if chi2[i] == None: 95 self.assertTrue(not hasattr(spin, 'chi2')) 96 97 # Data present. 98 else: 99 self.assertAlmostEqual(spin.chi2, chi2[i]) 100 self.assertAlmostEqual(spin.rx, rx[i]) 101 self.assertAlmostEqual(spin.i0/1e6, i0[i]/1e6) 102 103 # Increment the spin index. 104 i = i + 1 105 if i >= 12: 106 break

107 108

109 - def check_curve_fitting_manual(self):

110 """Check the results of the curve-fitting.""" 111 112 # Data. 113 relax_times = [0.0176, 0.0176, 0.0352, 0.0704, 0.0704, 0.1056, 0.1584, 0.1584, 0.1936, 0.1936] 114 chi2 = [2.916952651567855, 5.4916923952919632, 16.21182245065274, 4.3591263759462926, 9.8925377583244316, 6.0238341559877782] 115 rx = [8.0814894819820662, 8.6478971039559642, 9.5710638183013845, 10.716551838066295, 11.143793935455122, 12.82875370075309] 116 i0 = [1996050.9679875025, 2068490.9458927638, 1611556.5194095275, 1362887.2331948928, 1877670.5623875158, 897044.17382064369] 117 118 # Some checks. 119 self.assertEqual(cdp.curve_type, 'exp') 120 self.assertEqual(cdp.int_method, ds.int_type) 121 self.assertEqual(len(cdp.relax_times), 10) 122 cdp_relax_times = sorted(cdp.relax_times.values()) 123 for i in range(10): 124 self.assertEqual(cdp_relax_times[i], relax_times[i]) 125 126 # Check the errors. 127 for key in cdp.sigma_I: 128 self.assertAlmostEqual(cdp.sigma_I[key], 10578.039482421433, 6) 129 self.assertAlmostEqual(cdp.var_I[key], 111894919.29166669) 130 131 # Spin data check. 132 i = 0 133 for spin in spin_loop(skip_desel=True): 134 self.assertAlmostEqual(spin.chi2, chi2[i]) 135 self.assertAlmostEqual(spin.rx, rx[i]) 136 self.assertAlmostEqual(spin.i0/1e6, i0[i]/1e6) 137 138 # Increment the spin index. 139 i = i + 1

140 141

142 - def test_auto_analysis_pipe_name(self):

143 """Test the auto-analysis when the data pipe name is incorrectly supplied.""" 144 145 # Reset 146 self.interpreter.reset() 147 148 # Create the 'rx' data pipe. 149 self.interpreter.pipe.create('rx', 'relax_fit') 150 151 # Execute the auto-analysis 152 self.assertRaises(RelaxNoPipeError, relax_fit.Relax_fit, pipe_name='test', pipe_bundle='R1 analysis', file_root='rx', results_dir='.', grid_inc=11, mc_sim_num=50, view_plots=True)

153 154

155 - def test_bug_12670_12679(self):

156 """Test the relaxation curve fitting, replicating U{bug #12670<https://web.archive.org/web/https://gna.org/bugs/?12670>} and U{bug #12679<https://web.archive.org/web/https://gna.org/bugs/?12679>}.""" 157 158 # Execute the script. 159 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'1UBQ_relax_fit.py') 160 161 # Open the intensities.agr file. 162 file = open(ds.tmpdir + sep + 'intensities.agr') 163 lines = file.readlines() 164 file.close() 165 166 # Loop over all lines. 167 for i in range(len(lines)): 168 # Find the "@target G0.S0" line. 169 if search('@target', lines[i]): 170 index = i + 2 171 172 # Split up the lines. 173 lines[i] = lines[i].split() 174 175 # Check some of the Grace data. 176 self.assertEqual(len(lines[index]), 3) 177 self.assertEqual(lines[index][0], '0.004000000000000') 178 self.assertEqual(lines[index][1], '487178.000000000000000') 179 self.assertEqual(lines[index][2], '20570.000000000000000')

180 181

182 - def test_bug_18789(self):

183 """Test for zero errors in Grace plots, replicating U{bug #18789<https://web.archive.org/web/https://gna.org/bugs/?18789>}.""" 184 185 # Execute the script. 186 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'curve_fitting'+sep+'bug_18789_no_grace_errors.py') 187 188 # Open the Grace file. 189 file = open(ds.tmpdir + sep + 'rx.agr') 190 lines = file.readlines() 191 file.close() 192 193 # Loop over all lines. 194 for i in range(len(lines)): 195 # Find the "@target G0.S0" line. 196 if search('@target', lines[i]): 197 index = i + 2 198 199 # Split up the lines. 200 lines[i] = lines[i].split() 201 202 # Check for zero errors. 203 self.assertEqual(len(lines[index]), 3) 204 self.assertNotEqual(float(lines[index][2]), 0.0) 205 self.assertNotEqual(float(lines[index+1][2]), 0.0)

206 207

208 - def test_bug_19887_curvefit_fail(self):

209 """Test for the failure of relaxation curve-fitting, replicating U{bug #19887<https://web.archive.org/web/https://gna.org/bugs/?19887>}.""" 210 211 # Execute the script. 212 try: 213 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'curve_fitting'+sep+'bug_19887_curvefit_fail.py') 214 215 # A RelaxError is expected (but assertRaises() does not work with the script_exec method). 216 except RelaxError: 217 pass

218 219

220 - def test_bug_23244_Iinf_graph(self):

221 """Test that the auto-analysis creates the Iinf grace graphs, replicating U{bug #23244<https://web.archive.org/web/https://gna.org/bugs/?23244>}.""" 222 223 # Execute the script. 224 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'curve_fitting'+sep+'relax_fit_inversion_recovery.py') 225 226 # Check that the Iinf parameter files have been created. 227 self.assertTrue(access(ds.tmpdir+sep+'i0.out', F_OK)) 228 self.assertTrue(access(ds.tmpdir+sep+'iinf.out', F_OK)) 229 self.assertTrue(access(ds.tmpdir+sep+'rx.out', F_OK)) 230 231 # Check that the Iinf grace graphs have been created. 232 self.assertTrue(access(ds.tmpdir+sep+'grace'+sep+'i0.agr', F_OK)) 233 self.assertTrue(access(ds.tmpdir+sep+'grace'+sep+'iinf.agr', F_OK)) 234 self.assertTrue(access(ds.tmpdir+sep+'grace'+sep+'rx.agr', F_OK))

235 236

237 - def test_curve_fitting_height(self):

238 """Test the relaxation curve fitting C modules.""" 239 240 # The intensity type. 241 ds.int_type = 'height' 242 243 # Execute the script. 244 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit.py') 245 246 # Check the curve-fitting results. 247 self.check_curve_fitting()

248 249

250 - def test_curve_fitting_height_estimate_error(self):

251 """Test the relaxation curve fitting C modules and estimate error.""" 252 253 # Reset 254 self.interpreter.reset() 255 256 # Create pipe. 257 pipe_name = 'base pipe' 258 pipe_bundle = 'relax_fit' 259 self.interpreter.pipe.create(pipe_name=pipe_name, bundle=pipe_bundle, pipe_type='relax_fit') 260 261 # The intensity type. 262 ds.int_type = 'height' 263 264 # Create the data pipe and load the base data. 265 data_path = status.install_path + sep+'test_suite'+sep+'shared_data'+sep+'curve_fitting' 266 267 # Create the spins 268 self.interpreter.spectrum.read_spins(file="T2_ncyc1_ave.list", dir=data_path) 269 270 # Relaxation times (in seconds). 271 times = [ 272 0.0176, 273 0.0176, 274 0.0352, 275 0.0704, 276 0.0704, 277 0.1056, 278 0.1584, 279 0.1584, 280 0.1936, 281 0.1936 282 ] 283 284 # Spectrum names. 285 names = [ 286 'T2_ncyc1_ave', 287 'T2_ncyc1b_ave', 288 'T2_ncyc2_ave', 289 'T2_ncyc4_ave', 290 'T2_ncyc4b_ave', 291 'T2_ncyc6_ave', 292 'T2_ncyc9_ave', 293 'T2_ncyc9b_ave', 294 'T2_ncyc11_ave', 295 'T2_ncyc11b_ave' 296 ] 297 298 299 # Loop over Spectrum names. 300 for i, sname in enumerate(names): 301 # Get the time. 302 time = times[i] 303 304 # Load the peak intensities. 305 self.interpreter.spectrum.read_intensities(file=sname+'.list', dir=data_path, spectrum_id=sname, int_method=ds.int_type) 306 307 # Set the relaxation times. 308 self.interpreter.relax_fit.relax_time(time=time, spectrum_id=sname) 309 310 self.interpreter.deselect.spin(':3,11,18,19,23,31,42,44,54,66,82,92,94,99,101,113,124,126,136,141,145,147,332,345,346,358,361') 311 312 GRID_INC = 11 313 MC_SIM = 3 314 results_dir = mkdtemp() 315 #min_method = 'simplex' 316 #min_method = 'BFGS' 317 min_method = 'newton' 318 319 # De select one more. 320 self.interpreter.deselect.spin(':512@ND2') 321 322 # Set the relaxation curve type. 323 self.interpreter.relax_fit.select_model('exp') 324 325 # Do automatic 326 if True: 327 relax_fit.Relax_fit(pipe_name=pipe_name, pipe_bundle=pipe_bundle, file_root='R2', results_dir=results_dir, grid_inc=GRID_INC, mc_sim_num=MC_SIM, view_plots=False) 328 329 else: 330 # Prepare for finding dublictes. 331 332 # Collect all times, and matching spectrum ID. 333 all_times = [] 334 all_id = [] 335 for spectrum_id in cdp.relax_times: 336 all_times.append(cdp.relax_times[spectrum_id]) 337 all_id.append(spectrum_id) 338 339 # Get the dublicates. 340 dublicates = [(val, [i for i in range(len(all_times)) if all_times[i] == val]) for val in all_times] 341 342 # Loop over the list of the mapping of times and duplications. 343 list_dub_mapping = [] 344 for i, dub in enumerate(dublicates): 345 # Get current spectum id. 346 cur_spectrum_id = all_id[i] 347 348 # Get the tuple of time and indexes of duplications. 349 time, list_index_occur = dub 350 351 # Collect mapping of index to id. 352 id_list = [] 353 if len(list_index_occur) > 1: 354 for list_index in list_index_occur: 355 id_list.append(all_id[list_index]) 356 357 # Store to list 358 list_dub_mapping.append((cur_spectrum_id, id_list)) 359 360 # Assign dublicates. 361 for spectrum_id, dub_pair in list_dub_mapping: 362 if len(dub_pair) > 0: 363 self.interpreter.spectrum.replicated(spectrum_ids=dub_pair) 364 365 # Test the number of replicates stored in cdp, is 4. 366 self.assertEqual(len(cdp.replicates), 4) 367 368 369 # Peak intensity error analysis. 370 self.interpreter.spectrum.error_analysis() 371 372 # Grid search. 373 self.interpreter.minimise.grid_search(inc=GRID_INC) 374 375 # Minimise. 376 self.interpreter.minimise.execute(min_method, scaling=False, constraints=False) 377 378 # Monte Carlo simulations. 379 self.interpreter.monte_carlo.setup(number=MC_SIM) 380 self.interpreter.monte_carlo.create_data() 381 self.interpreter.monte_carlo.initial_values() 382 self.interpreter.minimise.execute(min_method, scaling=False, constraints=False) 383 self.interpreter.monte_carlo.error_analysis() 384 385 # Test seq 386 tseq = [ [4, 'GLY', ':4@N'], 387 [5, 'SER', ':5@N'], 388 [6, 'MET', ':6@N'], 389 [7, 'ASP', ':7@N'], 390 [8, 'SER', ':8@N'], 391 [12, 'GLY', ':12@N']] 392 393 # Print spins 394 i = 0 395 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True): 396 print(resi, resn, spin_id) 397 self.assertEqual(resi, tseq[i][0]) 398 self.assertEqual(resn, tseq[i][1]) 399 self.assertEqual(spin_id, tseq[i][2]) 400 401 i += 1 402 403 # Test the number of spins. 404 self.assertEqual(count_spins(), 6) 405 406 # Check the curve-fitting results. 407 self.check_curve_fitting_manual() 408 409 # Compare rx errors. 410 if True: 411 # Estimate rx and i0 errors. 412 self.interpreter.error_analysis.covariance_matrix() 413 414 # Collect: 415 i0_est = [] 416 i0_err_est = [] 417 rx_est = [] 418 rx_err_est = [] 419 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True): 420 i0_est.append(cur_spin.i0) 421 i0_err_est.append(cur_spin.i0_err) 422 rx_est.append(cur_spin.rx) 423 rx_err_est.append(cur_spin.rx_err) 424 425 # Set number of MC simulati0ns 426 MC_SIM = 200 427 428 # Monte Carlo simulations. 429 self.interpreter.monte_carlo.setup(number=MC_SIM) 430 self.interpreter.monte_carlo.create_data() 431 self.interpreter.monte_carlo.initial_values() 432 self.interpreter.minimise.execute(min_method, scaling=False, constraints=False) 433 self.interpreter.monte_carlo.error_analysis() 434 435 # Collect: 436 i0_mc = [] 437 i0_err_mc = [] 438 rx_mc = [] 439 rx_err_mc = [] 440 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True): 441 i0_mc.append(cur_spin.i0) 442 i0_err_mc.append(cur_spin.i0_err) 443 rx_mc.append(cur_spin.rx) 444 rx_err_mc.append(cur_spin.rx_err) 445 446 # Now print and compare 447 i = 0 448 print("Comparison between error estimation from Jacobian co-variance matrix and Monte-Carlo simulations.") 449 print("Spin ID: rx_err_diff=est-MC, i0_err_diff=est-MC, rx_err=est/MC, i0_err=est/MC, i0=est/MC, rx=est/MC.") 450 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True): 451 # Extract for estimation. 452 i0_est_i = i0_est[i] 453 i0_err_est_i = i0_err_est[i] 454 rx_est_i = rx_est[i] 455 rx_err_est_i = rx_err_est[i] 456 457 # Extract from monte carlo. 458 i0_mc_i = i0_mc[i] 459 i0_err_mc_i = i0_err_mc[i] 460 rx_mc_i = rx_mc[i] 461 rx_err_mc_i = rx_err_mc[i] 462 463 # Add to counter. 464 i += 1 465 466 # Prepare text. 467 rx_err_diff = rx_err_est_i - rx_err_mc_i 468 i0_err_diff = i0_err_est_i - i0_err_mc_i 469 470 text = "Spin '%s': rx_err_diff=%3.4f, i0_err_diff=%3.3f, rx_err=%3.4f/%3.4f, i0_err=%3.3f/%3.3f, rx=%3.3f/%3.3f, i0=%3.3f/%3.3f" % (spin_id, rx_err_diff, i0_err_diff, rx_err_est_i, rx_err_mc_i, i0_err_est_i, i0_err_mc_i, rx_est_i, rx_mc_i, i0_est_i, i0_mc_i) 471 print(text)

472 473

474 - def test_curve_fitting_volume(self):

475 """Test the relaxation curve fitting C modules.""" 476 477 # The intensity type. 478 ds.int_type = 'point sum' 479 480 # Execute the script. 481 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit.py') 482 483 # Check the curve-fitting results. 484 self.check_curve_fitting()

485 486

487 - def test_inversion_recovery(self):

488 """Test the fitting for the inversion recovery R1 experiment.""" 489 490 # Execute the script. 491 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit_exp_3param_inv_neg.py') 492 493 # Check the curve-fitting results. 494 spin = return_spin(spin_id=":4@N") 495 self.assertAlmostEqual(spin.rx, 1.2) 496 self.assertAlmostEqual(spin.i0, -30.0) 497 self.assertAlmostEqual(spin.iinf, 22.0)

498 499

500 - def test_read_sparky(self):

501 """The Sparky peak height loading test.""" 502 503 # Load the original state. 504 self.interpreter.state.load(state='basic_heights_T2_ncyc1', dir=status.install_path + sep+'test_suite'+sep+'shared_data'+sep+'saved_states', force=True) 505 506 # Create a new data pipe for the new data. 507 self.interpreter.pipe.create('new', 'relax_fit') 508 509 # Load the Lupin Ap4Aase sequence. 510 self.interpreter.sequence.read(file="Ap4Aase.seq", dir=status.install_path + sep+'test_suite'+sep+'shared_data', res_num_col=1, res_name_col=2) 511 512 # Name the spins so they can be matched to the assignments. 513 self.interpreter.spin.name(name='N') 514 515 # Read the peak heights. 516 self.interpreter.spectrum.read_intensities(file="T2_ncyc1_ave.list", dir=status.install_path + sep+'test_suite'+sep+'shared_data'+sep+'curve_fitting', spectrum_id='0.0176') 517 518 519 # Test the integrity of the data. 520 ################################# 521 522 # Get the data pipes. 523 dp_new = pipes.get_pipe('new') 524 dp_rx = pipes.get_pipe('rx') 525 526 # Loop over the spins of the original data. 527 for mol_index, res_index, spin_index in spin_index_loop(): 528 # Alias the spin containers. 529 new_spin = dp_new.mol[mol_index].res[res_index].spin[spin_index] 530 orig_spin = dp_rx.mol[mol_index].res[res_index].spin[spin_index] 531 532 # Check the sequence info. 533 self.assertEqual(dp_new.mol[mol_index].name, dp_rx.mol[mol_index].name) 534 self.assertEqual(dp_new.mol[mol_index].res[res_index].num, dp_rx.mol[mol_index].res[res_index].num) 535 self.assertEqual(dp_new.mol[mol_index].res[res_index].name, dp_rx.mol[mol_index].res[res_index].name) 536 self.assertEqual(new_spin.num, orig_spin.num) 537 self.assertEqual(new_spin.name, orig_spin.name) 538 539 # Skip deselected spins. 540 if not orig_spin.select: 541 continue 542 543 # Check intensities (if they exist). 544 if hasattr(orig_spin, 'peak_intensity'): 545 for id in dp_new.spectrum_ids: 546 self.assertEqual(orig_spin.peak_intensity[id], new_spin.peak_intensity[id])

547 548

549 - def test_saturation_recovery(self):

550 """Test the fitting for the saturation recovery R1 experiment.""" 551 552 # Execute the script. 553 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit_saturation_recovery.py') 554 555 # Check the curve-fitting results. 556 spin = return_spin(spin_id=":17@H") 557 self.assertAlmostEqual(spin.chi2, 0.0) 558 self.assertAlmostEqual(spin.rx, 0.5) 559 self.assertAlmostEqual(spin.iinf/1e15, 1.0)

560 561

562 - def test_zooming_grid_search(self):

563 """Test the relaxation curve fitting C modules.""" 564 565 # Execute the script. 566 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit_zooming_grid.py') 567 568 # Check the curve-fitting results (the values are from the optimisation of test_curve_fitting_height()). 569 spin = return_spin(spin_id=":4@N") 570 self.assertAlmostEqual(spin.chi2, 2.9169526515678883) 571 self.assertAlmostEqual(spin.rx, 8.0814894974893967) 572 self.assertAlmostEqual(spin.i0/1e6, 1996050.9699629977/1e6)

Source Code for Module test_suite.system_tests.relax_fit