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.assert_(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_bug_12670_12679(self):

143 """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>}.""" 144 145 # Execute the script. 146 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'1UBQ_relax_fit.py') 147 148 # Open the intensities.agr file. 149 file = open(ds.tmpdir + sep + 'intensities.agr') 150 lines = file.readlines() 151 file.close() 152 153 # Loop over all lines. 154 for i in range(len(lines)): 155 # Find the "@target G0.S0" line. 156 if search('@target', lines[i]): 157 index = i + 2 158 159 # Split up the lines. 160 lines[i] = lines[i].split() 161 162 # Check some of the Grace data. 163 self.assertEqual(len(lines[index]), 3) 164 self.assertEqual(lines[index][0], '0.004000000000000') 165 self.assertEqual(lines[index][1], '487178.000000000000000') 166 self.assertEqual(lines[index][2], '20570.000000000000000')

167 168

169 - def test_bug_18789(self):

170 """Test for zero errors in Grace plots, replicating U{bug #18789<https://web.archive.org/web/https://gna.org/bugs/?18789>}.""" 171 172 # Execute the script. 173 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'curve_fitting'+sep+'bug_18789_no_grace_errors.py') 174 175 # Open the Grace file. 176 file = open(ds.tmpdir + sep + 'rx.agr') 177 lines = file.readlines() 178 file.close() 179 180 # Loop over all lines. 181 for i in range(len(lines)): 182 # Find the "@target G0.S0" line. 183 if search('@target', lines[i]): 184 index = i + 2 185 186 # Split up the lines. 187 lines[i] = lines[i].split() 188 189 # Check for zero errors. 190 self.assertEqual(len(lines[index]), 3) 191 self.assertNotEqual(float(lines[index][2]), 0.0) 192 self.assertNotEqual(float(lines[index+1][2]), 0.0)

193 194

195 - def test_bug_19887_curvefit_fail(self):

196 """Test for the failure of relaxation curve-fitting, replicating U{bug #19887<https://web.archive.org/web/https://gna.org/bugs/?19887>}.""" 197 198 # Execute the script. 199 try: 200 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'curve_fitting'+sep+'bug_19887_curvefit_fail.py') 201 202 # A RelaxError is expected (but assertRaises() does not work with the script_exec method). 203 except RelaxError: 204 pass

205 206

207 - def test_bug_23244_Iinf_graph(self):

208 """Test that the auto-analysis creates the Iinf grace graphs, replicating U{bug #23244<https://web.archive.org/web/https://gna.org/bugs/?23244>}.""" 209 210 # Execute the script. 211 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'curve_fitting'+sep+'relax_fit_inversion_recovery.py') 212 213 # Check that the Iinf parameter files have been created. 214 self.assert_(access(ds.tmpdir+sep+'i0.out', F_OK)) 215 self.assert_(access(ds.tmpdir+sep+'iinf.out', F_OK)) 216 self.assert_(access(ds.tmpdir+sep+'rx.out', F_OK)) 217 218 # Check that the Iinf grace graphs have been created. 219 self.assert_(access(ds.tmpdir+sep+'grace'+sep+'i0.agr', F_OK)) 220 self.assert_(access(ds.tmpdir+sep+'grace'+sep+'iinf.agr', F_OK)) 221 self.assert_(access(ds.tmpdir+sep+'grace'+sep+'rx.agr', F_OK))

222 223

224 - def test_curve_fitting_height(self):

225 """Test the relaxation curve fitting C modules.""" 226 227 # The intensity type. 228 ds.int_type = 'height' 229 230 # Execute the script. 231 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit.py') 232 233 # Check the curve-fitting results. 234 self.check_curve_fitting()

235 236

237 - def test_curve_fitting_height_estimate_error(self):

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

459 460

461 - def test_curve_fitting_volume(self):

462 """Test the relaxation curve fitting C modules.""" 463 464 # The intensity type. 465 ds.int_type = 'point sum' 466 467 # Execute the script. 468 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit.py') 469 470 # Check the curve-fitting results. 471 self.check_curve_fitting()

472 473

474 - def test_inversion_recovery(self):

475 """Test the fitting for the inversion recovery R1 experiment.""" 476 477 # Execute the script. 478 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit_exp_3param_inv_neg.py') 479 480 # Check the curve-fitting results. 481 spin = return_spin(":4@N") 482 self.assertAlmostEqual(spin.rx, 1.2) 483 self.assertAlmostEqual(spin.i0, -30.0) 484 self.assertAlmostEqual(spin.iinf, 22.0)

485 486

487 - def test_read_sparky(self):

488 """The Sparky peak height loading test.""" 489 490 # Load the original state. 491 self.interpreter.state.load(state='basic_heights_T2_ncyc1', dir=status.install_path + sep+'test_suite'+sep+'shared_data'+sep+'saved_states', force=True) 492 493 # Create a new data pipe for the new data. 494 self.interpreter.pipe.create('new', 'relax_fit') 495 496 # Load the Lupin Ap4Aase sequence. 497 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) 498 499 # Name the spins so they can be matched to the assignments. 500 self.interpreter.spin.name(name='N') 501 502 # Read the peak heights. 503 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') 504 505 506 # Test the integrity of the data. 507 ################################# 508 509 # Get the data pipes. 510 dp_new = pipes.get_pipe('new') 511 dp_rx = pipes.get_pipe('rx') 512 513 # Loop over the spins of the original data. 514 for mol_index, res_index, spin_index in spin_index_loop(): 515 # Alias the spin containers. 516 new_spin = dp_new.mol[mol_index].res[res_index].spin[spin_index] 517 orig_spin = dp_rx.mol[mol_index].res[res_index].spin[spin_index] 518 519 # Check the sequence info. 520 self.assertEqual(dp_new.mol[mol_index].name, dp_rx.mol[mol_index].name) 521 self.assertEqual(dp_new.mol[mol_index].res[res_index].num, dp_rx.mol[mol_index].res[res_index].num) 522 self.assertEqual(dp_new.mol[mol_index].res[res_index].name, dp_rx.mol[mol_index].res[res_index].name) 523 self.assertEqual(new_spin.num, orig_spin.num) 524 self.assertEqual(new_spin.name, orig_spin.name) 525 526 # Skip deselected spins. 527 if not orig_spin.select: 528 continue 529 530 # Check intensities (if they exist). 531 if hasattr(orig_spin, 'peak_intensity'): 532 for id in dp_new.spectrum_ids: 533 self.assertEqual(orig_spin.peak_intensity[id], new_spin.peak_intensity[id])

534 535

536 - def test_saturation_recovery(self):

537 """Test the fitting for the saturation recovery R1 experiment.""" 538 539 # Execute the script. 540 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit_saturation_recovery.py') 541 542 # Check the curve-fitting results. 543 spin = return_spin(":17@H") 544 self.assertAlmostEqual(spin.chi2, 0.0) 545 self.assertAlmostEqual(spin.rx, 0.5) 546 self.assertAlmostEqual(spin.iinf/1e15, 1.0)

547 548

549 - def test_zooming_grid_search(self):

550 """Test the relaxation curve fitting C modules.""" 551 552 # Execute the script. 553 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit_zooming_grid.py') 554 555 # Check the curve-fitting results (the values are from the optimisation of test_curve_fitting_height()). 556 spin = return_spin(":4@N") 557 self.assertAlmostEqual(spin.chi2, 2.9169526515678883) 558 self.assertAlmostEqual(spin.rx, 8.0814894974893967) 559 self.assertAlmostEqual(spin.i0/1e6, 1996050.9699629977/1e6)

Source Code for Module test_suite.system_tests.relax_fit