1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 from os import F_OK, access, sep
25 from re import search
26 from tempfile import mkdtemp
27
28
29 from auto_analyses import relax_fit
30 from data_store import Relax_data_store; ds = Relax_data_store()
31 import dep_check
32 from pipe_control.mol_res_spin import count_spins, return_spin, spin_index_loop, spin_loop
33 from pipe_control import pipes
34 from lib.errors import RelaxError
35 from status import Status; status = Status()
36 from test_suite.system_tests.base_classes import SystemTestCase
37
38
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
50 super(Relax_fit, self).__init__(methodName)
51
52
53 if not dep_check.C_module_exp_fn:
54
55 status.skipped_tests.append([methodName, 'Relax curve-fitting C module', self._skip_type])
56
57
59 """Set up for all the functional tests."""
60
61
62 self.interpreter.pipe.create('mf', 'mf')
63
64
65 ds.tmpdir = mkdtemp()
66
67
69 """Check the results of the curve-fitting."""
70
71
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
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
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
91 i = 0
92 for spin in spin_loop():
93
94 if chi2[i] == None:
95 self.assert_(not hasattr(spin, 'chi2'))
96
97
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
104 i = i + 1
105 if i >= 12:
106 break
107
108
110 """Check the results of the curve-fitting."""
111
112
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
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
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
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
139 i = i + 1
140
141
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
146 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'1UBQ_relax_fit.py')
147
148
149 file = open(ds.tmpdir + sep + 'intensities.agr')
150 lines = file.readlines()
151 file.close()
152
153
154 for i in range(len(lines)):
155
156 if search('@target', lines[i]):
157 index = i + 2
158
159
160 lines[i] = lines[i].split()
161
162
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
170 """Test for zero errors in Grace plots, replicating U{bug #18789<https://web.archive.org/web/https://gna.org/bugs/?18789>}."""
171
172
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
176 file = open(ds.tmpdir + sep + 'rx.agr')
177 lines = file.readlines()
178 file.close()
179
180
181 for i in range(len(lines)):
182
183 if search('@target', lines[i]):
184 index = i + 2
185
186
187 lines[i] = lines[i].split()
188
189
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
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
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
203 except RelaxError:
204 pass
205
206
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
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
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
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
225 """Test the relaxation curve fitting C modules."""
226
227
228 ds.int_type = 'height'
229
230
231 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit.py')
232
233
234 self.check_curve_fitting()
235
236
238 """Test the relaxation curve fitting C modules and estimate error."""
239
240
241 self.interpreter.reset()
242
243
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
249 ds.int_type = 'height'
250
251
252 data_path = status.install_path + sep+'test_suite'+sep+'shared_data'+sep+'curve_fitting'
253
254
255 self.interpreter.spectrum.read_spins(file="T2_ncyc1_ave.list", dir=data_path)
256
257
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
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
287 for i, sname in enumerate(names):
288
289 time = times[i]
290
291
292 self.interpreter.spectrum.read_intensities(file=sname+'.list', dir=data_path, spectrum_id=sname, int_method=ds.int_type)
293
294
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
303
304 min_method = 'newton'
305
306
307 self.interpreter.deselect.spin(':512@ND2')
308
309
310 self.interpreter.relax_fit.select_model('exp')
311
312
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
318
319
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
327 dublicates = [(val, [i for i in range(len(all_times)) if all_times[i] == val]) for val in all_times]
328
329
330 list_dub_mapping = []
331 for i, dub in enumerate(dublicates):
332
333 cur_spectrum_id = all_id[i]
334
335
336 time, list_index_occur = dub
337
338
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
345 list_dub_mapping.append((cur_spectrum_id, id_list))
346
347
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
353 self.assertEqual(len(cdp.replicates), 4)
354
355
356
357 self.interpreter.spectrum.error_analysis()
358
359
360 self.interpreter.minimise.grid_search(inc=GRID_INC)
361
362
363 self.interpreter.minimise.execute(min_method, scaling=False, constraints=False)
364
365
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
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
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
391 self.assertEqual(count_spins(), 6)
392
393
394 self.check_curve_fitting_manual()
395
396
397 if True:
398
399 self.interpreter.error_analysis.covariance_matrix()
400
401
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
413 MC_SIM = 200
414
415
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
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
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
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
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
451 i += 1
452
453
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
462 """Test the relaxation curve fitting C modules."""
463
464
465 ds.int_type = 'point sum'
466
467
468 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit.py')
469
470
471 self.check_curve_fitting()
472
473
475 """Test the fitting for the inversion recovery R1 experiment."""
476
477
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
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
488 """The Sparky peak height loading test."""
489
490
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
494 self.interpreter.pipe.create('new', 'relax_fit')
495
496
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
500 self.interpreter.spin.name(name='N')
501
502
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
507
508
509
510 dp_new = pipes.get_pipe('new')
511 dp_rx = pipes.get_pipe('rx')
512
513
514 for mol_index, res_index, spin_index in spin_index_loop():
515
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
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
527 if not orig_spin.select:
528 continue
529
530
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
537 """Test the fitting for the saturation recovery R1 experiment."""
538
539
540 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit_saturation_recovery.py')
541
542
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
550 """Test the relaxation curve fitting C modules."""
551
552
553 self.script_exec(status.install_path + sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_fit_zooming_grid.py')
554
555
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)
560