1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """The automatic relaxation dispersion protocol for repeated data for CPMG.
24
25 U{task #7826<https://web.archive.org/web/https://gna.org/task/index.php?78266>}, Write an python class for the repeated analysis of dispersion data.
26 """
27
28
29 import dep_check
30
31
32 from copy import deepcopy
33 from datetime import datetime
34 from glob import glob
35 from os import F_OK, access, chmod, getcwd, sep
36 from numpy import asarray, arange, concatenate, max, mean, min, savetxt, square, sqrt, std, sum
37 if dep_check.scipy_module:
38 from scipy.stats import pearsonr
39 from stat import S_IRWXU, S_IRGRP, S_IROTH
40 import sys
41 from warnings import warn
42
43
44 import dep_check
45 from lib.dispersion.variables import MODEL_NOREX, MODEL_PARAMS, MODEL_R2EFF, PARAMS_R20
46 from lib.io import extract_data, get_file_path, open_write_file, sort_filenames, write_data
47 from lib.text.sectioning import section, subsection, subtitle
48 from lib.warnings import RelaxWarning
49 from pipe_control.mol_res_spin import spin_loop
50 from pipe_control import pipes
51 from prompt.interpreter import Interpreter
52 from specific_analyses.relax_disp.data import generate_r20_key, has_exponential_exp_type, has_cpmg_exp_type, is_r1_optimised, loop_exp_frq_offset, loop_exp_frq_offset_point, return_param_key_from_data
53 from status import Status; status = Status()
54
55 if dep_check.matplotlib_module:
56 import pylab as plt
57 from matplotlib.font_manager import FontProperties
58 fontP = FontProperties()
59 fontP.set_size('small')
60
61
62
63 DIC_KEY_FORMAT = "%.8f"
64
65
67 """Calculate the linear correlation 'B', the intercept 'A' for the function y=A + Bx.
68
69 @keyword x: The data for the X-axis.
70 @type x: float or numpy array.
71 @keyword y: The data for the Y-axis.
72 @type y: float or numpy array.
73 @return: The intercept A, the standard deviation for A, the slope B, the standard deviation for B, standard deviation of the residuals, the linear correlation coefficient
74 @rtype: float, float, float, float, float, float
75 """
76
77
78 x_m = mean(x)
79 y_m = mean(y)
80
81
82 N = len(y)
83
84
85 denom = N * sum(x**2) - (sum(x))**2
86
87
88 A = ( sum(x**2) * sum(y) - sum(x) * sum(x*y) ) / denom
89 B = (N * sum(x*y) - sum(x) * sum(y) ) / denom
90
91
92 std_y = sqrt( (1. / (N-2) ) * sum( (y - A - B*x)**2 ) )
93
94
95
96 std_A = std_y * sqrt( sum(x**2) / denom )
97 std_B = std_y * sqrt( N / denom )
98
99
100 r_xy = sum( (x - x_m)*(y - y_m) ) / sqrt( sum((x - x_m)**2) * sum((y - y_m)**2) )
101
102 return A, std_A, B, std_B, std_y, r_xy
103
104
106
107 """The relaxation dispersion analysis for repeated data."""
108
109
110 opt_func_tol = 1e-25
111 opt_max_iterations = int(1e7)
112
113
115 """Perform a repeated dispersion analysis for settings given."""
116
117
118 self.settings = settings
119
120
121 for setting in self.settings:
122 setattr(self, setting, self.settings[setting])
123
124 if 'pipe_type' not in self.settings:
125 self.set_self(key='pipe_type', value='relax_disp')
126
127 if 'time' not in self.settings:
128 self.set_self(key='time', value=datetime.now().strftime('%Y_%m_%d_%H_%M'))
129
130
131 if 'results_dir' not in self.settings:
132 self.set_self(key='results_dir', value=getcwd() + sep + 'results' + sep + self.time )
133
134
135 if 'mc_sim_num' not in self.settings:
136 self.set_self(key='mc_sim_num', value=40)
137
138
139 if 'exp_mc_sim_num' not in self.settings:
140 self.set_self(key='exp_mc_sim_num', value=-1)
141
142
143 if 'modsel' not in self.settings:
144 self.set_self(key='modsel', value='AIC')
145
146
147 if 'insignificance' not in self.settings:
148 self.set_self(key='insignificance', value=0.0)
149
150
151
152 if 'r1_fit' not in self.settings:
153 self.set_self(key='r1_fit', value=False)
154
155
156 if 'min_algor' not in self.settings:
157 self.set_self(key='min_algor', value='simplex')
158
159
160 if 'constraints' not in self.settings:
161 self.set_self(key='constraints', value=True)
162
163
164 if 'base_setup_pipe_name' not in self.settings:
165 base_setup_pipe_name = self.name_pipe(method='setup', model='setup', analysis='setup', glob_ini='setup')
166 self.set_self(key='base_setup_pipe_name', value=base_setup_pipe_name)
167
168
169 self.interpreter_start()
170
171
172 - def set_base_cpmg(self, method=None, glob_ini=None, force=False):
173 """ Setup base information, but do not load intensity. """
174
175
176 model = 'setup'
177 analysis = 'setup'
178
179
180 found, pipe_name, resfile, path = self.check_previous_result(method='setup', model=model, analysis=analysis, glob_ini='setup', bundle='setup')
181
182
183 if found:
184 pass
185 else:
186
187 self.interpreter.pipe.create(pipe_name=pipe_name, pipe_type=self.pipe_type, bundle=None)
188
189
190 dic_spectrum_ids = {}
191 dic_spectrum_ids_replicates = {}
192 for i, sfrq in enumerate(self.sfrqs):
193
194 key = DIC_KEY_FORMAT % (sfrq)
195
196
197 cpmg_frqs = getattr(self, key)['cpmg_frqs']
198
199
200 peaks_folder = getattr(self, key)['peaks_folder'] + sep + method
201
202
203 peaks_glob_pat = '%s_%s.ser' % (glob_ini, method)
204
205
206 peaks_file_list = glob(peaks_folder + sep + peaks_glob_pat)
207
208
209 peaks_file_list = sort_filenames(filenames=peaks_file_list)
210
211
212 for peaks_file in peaks_file_list:
213 self.interpreter.spectrum.read_spins(file=peaks_file, dir=None)
214
215
216 dic_spectrum_ids[key] = []
217 for j, cpmg_frq in enumerate(cpmg_frqs):
218
219 data_key = return_param_key_from_data(exp_type=self.exp_type, frq=sfrq, point=cpmg_frq)
220 spectrum_id = data_key + '_%i'%j
221
222
223 dic_spectrum_ids[key].append(spectrum_id)
224
225
226 self.interpreter.relax_disp.exp_type(spectrum_id=spectrum_id, exp_type=self.exp_type)
227
228
229 if cpmg_frq == 0.0:
230 cpmg_frq = None
231 self.interpreter.relax_disp.cpmg_setup(spectrum_id=spectrum_id, cpmg_frq=cpmg_frq)
232
233
234 time_T2 = getattr(self, key)['time_T2']
235 self.interpreter.relax_disp.relax_time(spectrum_id=spectrum_id, time=time_T2)
236
237
238 self.interpreter.spectrometer.frequency(id=spectrum_id, frq=sfrq, units=self.sfrq_unit)
239
240
241 list_dub = self.get_dublicates(dic_spectrum_ids[key], cpmg_frqs)
242
243
244 dic_spectrum_ids_replicates[key] = list_dub
245
246
247 cdp.dic_spectrum_ids = dic_spectrum_ids
248 cdp.dic_spectrum_ids_replicates = dic_spectrum_ids_replicates
249
250
251 self.interpreter.spin.isotope(isotope=self.isotope)
252
253
254 self.interpreter.results.write(file=resfile, dir=path, force=force)
255
256
258
259
260
261 if pipes.cdp_name() != pipe_name:
262 self.interpreter.pipe.switch(pipe_name)
263
264
265 finished = len(self.sfrqs) * [False]
266 for i, sfrq in enumerate(self.sfrqs):
267
268 key = DIC_KEY_FORMAT % (sfrq)
269
270
271 spectrum_ids = cdp.dic_spectrum_ids[key]
272
273
274 peaks_folder = getattr(self, key)['peaks_folder'] + sep + self.method
275
276
277 peaks_glob_pat = '%s_%s.ser' % (glob_ini, self.method)
278
279
280 peaks_file_list = glob(peaks_folder + sep + peaks_glob_pat)
281
282
283 peaks_file_list = sort_filenames(filenames=peaks_file_list)
284
285
286 if len(peaks_file_list) == 0:
287 finished[i] = False
288 continue
289
290
291 for peaks_file in peaks_file_list:
292 self.interpreter.spectrum.read_intensities(file=peaks_file, spectrum_id=spectrum_ids, int_method=self.int_method, int_col=list(range(len(spectrum_ids))))
293
294 if set_rmsd:
295
296 rmsd_folder = getattr(self, key)['rmsd_folder']
297
298
299 rmsd_glob_pat = '%s_*_%s.rmsd' % (glob_ini, self.method)
300
301
302 rmsd_file_list = glob(rmsd_folder + sep + rmsd_glob_pat)
303
304
305 rmsd_file_list = sort_filenames(filenames=rmsd_file_list)
306
307
308 for j, spectrum_id in enumerate(spectrum_ids):
309
310 rmsd_file = rmsd_file_list[j]
311
312
313 rmsd = float(extract_data(file=rmsd_file)[0][0])
314 self.interpreter.spectrum.baseplane_rmsd(error=rmsd, spectrum_id=spectrum_id)
315
316 finished[i] = True
317
318 return all(finished)
319
320
322 """Do spectrum error analysis, where both replicates per spectrometer frequency and subset is taken into consideration."""
323
324
325
326 if pipes.cdp_name() != pipe_name:
327 self.interpreter.pipe.switch(pipe_name)
328
329
330 for i, sfrq in enumerate(self.sfrqs):
331
332 key = DIC_KEY_FORMAT % (sfrq)
333
334
335 section(file=sys.stdout, text="Error analysis for pipe='%s' and sfr:%3.2f"%(pipe_name, sfrq), prespace=2)
336
337
338 spectrum_ids = cdp.dic_spectrum_ids[key]
339
340 if set_rep:
341
342 spectrum_ids_replicates = cdp.dic_spectrum_ids_replicates[key]
343
344
345 for replicate in spectrum_ids_replicates:
346 spectrum_id, rep_list = replicate
347
348
349 if len(rep_list) > 0:
350
351 self.interpreter.spectrum.replicated(spectrum_ids=rep_list)
352
353
354 self.interpreter.spectrum.error_analysis(subset=spectrum_ids)
355
356
357 - def set_int(self, methods=None, list_glob_ini=None, set_rmsd=True, set_rep=False, force=False):
358 """Call both the setup of data and the error analysis"""
359
360
361 model = 'setup'
362 analysis = 'int'
363
364
365 finished = len(methods) * [False]
366 for i, method in enumerate(methods):
367
368 self.set_self(key='method', value=method)
369
370
371 for glob_ini in list_glob_ini:
372
373 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method)
374
375 if not found:
376 calculate = True
377 finished[i] = False
378 elif found:
379 calculate = False
380 finished[i] = True
381
382 if calculate:
383
384 self.interpreter.pipe.copy(pipe_from=self.base_setup_pipe_name, pipe_to=pipe_name, bundle_to=self.method)
385 self.interpreter.pipe.switch(pipe_name)
386
387
388 finished_int = self.set_intensity_and_error(pipe_name=pipe_name, glob_ini=glob_ini, set_rmsd=set_rmsd)
389
390 if finished_int:
391
392 self.do_spectrum_error_analysis(pipe_name=pipe_name, set_rep=set_rep)
393
394
395 cdp.settings = self.settings
396 self.interpreter.results.write(file=resfile, dir=path, force=force)
397
398 finished[i] = True
399
400 else:
401 pipe_name = pipes.cdp_name()
402 self.interpreter.pipe.delete(pipe_name=pipe_name)
403
404 return all(finished)
405
406
407 - def calc_r2eff(self, methods=None, list_glob_ini=None, force=False):
408 """Method to calculate R2eff or read previous results."""
409
410 model = MODEL_R2EFF
411 analysis = 'int'
412
413
414 for method in methods:
415
416 self.set_self(key='method', value=method)
417
418
419 for glob_ini in list_glob_ini:
420
421 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method)
422
423 if not found:
424 calculate = True
425 elif found:
426 calculate = False
427
428 if calculate:
429
430
431 intensity_pipe_name = self.name_pipe(method=self.method, model='setup', analysis='int', glob_ini=glob_ini)
432
433 if not pipes.has_pipe(intensity_pipe_name):
434 finished = self.set_int(methods=[method], list_glob_ini=[glob_ini])
435 if not finished:
436 continue
437
438 self.interpreter.pipe.copy(pipe_from=intensity_pipe_name, pipe_to=pipe_name, bundle_to=self.method)
439 self.interpreter.pipe.switch(pipe_name)
440
441
442 self.interpreter.relax_disp.select_model(model)
443
444
445 subtitle(file=sys.stdout, text="The '%s' model for pipe='%s'" % (model, pipe_name), prespace=3)
446
447
448 if model == MODEL_R2EFF and not has_exponential_exp_type():
449 self.interpreter.minimise.calculate()
450
451
452 cdp.settings = self.settings
453 self.interpreter.results.write(file=resfile, dir=path, force=force)
454
455
456 - def deselect_all(self, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):
457 """Method to deselect all spins for a pipe."""
458
459
460 if model_from == None:
461 model_from = model
462 if analysis_from == None:
463 analysis_from = analysis
464
465
466 for method in methods:
467
468 self.set_self(key='method', value=method)
469
470
471 for glob_ini in list_glob_ini:
472
473 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method)
474
475 if not found:
476
477 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini)
478
479
480 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method)
481 self.interpreter.pipe.switch(pipe_name)
482
483
484 subtitle(file=sys.stdout, text="Deselect all spins for pipe='%s'" % (pipe_name), prespace=3)
485
486
487 self.interpreter.deselect.all()
488
489
490 cdp.settings = self.settings
491
492 if found and not force:
493 file_path = get_file_path(file_name=resfile, dir=path)
494 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path)
495 warn(RelaxWarning(text))
496 else:
497 self.interpreter.results.write(file=resfile, dir=path, force=force)
498
499
500 self.interpreter.spin.display()
501
502
503 - def select_spin(self, spin_id=None, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):
504 """Method to select spins for a pipe."""
505
506
507 if model_from == None:
508 model_from = model
509 if analysis_from == None:
510 analysis_from = analysis
511
512
513 for method in methods:
514
515 self.set_self(key='method', value=method)
516
517
518 for glob_ini in list_glob_ini:
519
520 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method)
521
522 if not found:
523
524 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini)
525
526
527 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method)
528 self.interpreter.pipe.switch(pipe_name)
529
530
531 subtitle(file=sys.stdout, text="Select spins '%s' for pipe='%s'" % (spin_id, pipe_name), prespace=3)
532
533
534 self.interpreter.select.spin(spin_id=spin_id)
535
536
537 cdp.settings = self.settings
538
539 if found and not force:
540 file_path = get_file_path(file_name=resfile, dir=path)
541 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path)
542 warn(RelaxWarning(text))
543 else:
544 self.interpreter.results.write(file=resfile, dir=path, force=force)
545
546
547 self.interpreter.spin.display()
548
549
550 - def value_set(self, spin_id=None, val=None, param=None, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):
551 """Use value.set on all pipes."""
552
553
554 if model_from == None:
555 model_from = model
556 if analysis_from == None:
557 analysis_from = analysis
558
559
560 for method in methods:
561
562 self.set_self(key='method', value=method)
563
564
565 for glob_ini in list_glob_ini:
566
567 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method)
568
569 if not found:
570
571 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini)
572
573
574 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method)
575 self.interpreter.pipe.switch(pipe_name)
576
577
578 subtitle(file=sys.stdout, text="For param '%s' set value '%3.2f' for pipe='%s'" % (param, val, pipe_name), prespace=3)
579
580
581 self.interpreter.relax_disp.select_model(model)
582
583
584 self.interpreter.value.set(val=val, param=param, spin_id=spin_id)
585
586
587 cdp.settings = self.settings
588
589 if found and not force:
590 file_path = get_file_path(file_name=resfile, dir=path)
591 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path)
592 warn(RelaxWarning(text))
593 else:
594 self.interpreter.results.write(file=resfile, dir=path, force=force)
595
596
597 self.spin_display_params(pipe_name=pipe_name)
598
599
600 - def r20_from_min_r2eff(self, spin_id=None, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):
601 """Will set the grid R20 values from the minimum R2eff values through the r20_from_min_r2eff user function.
602 This will speed up the grid search with a factor GRID_INC^(Nr_spec_freq). For a CPMG experiment with two fields and standard GRID_INC = 21, the speed-up is a factor 441.
603 """
604
605
606 if model_from == None:
607 model_from = model
608 if analysis_from == None:
609 analysis_from = analysis
610
611
612 for method in methods:
613
614 self.set_self(key='method', value=method)
615
616
617 for glob_ini in list_glob_ini:
618
619 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method)
620
621 if not found:
622
623 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini)
624
625
626 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method)
627 self.interpreter.pipe.switch(pipe_name)
628
629
630 subtitle(file=sys.stdout, text="Set grid r20 for pipe='%s'" % (pipe_name), prespace=3)
631
632
633 self.interpreter.relax_disp.select_model(model)
634
635
636 self.interpreter.relax_disp.r20_from_min_r2eff(force=True)
637
638
639 cdp.settings = self.settings
640
641 if found and not force:
642 file_path = get_file_path(file_name=resfile, dir=path)
643 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path)
644 warn(RelaxWarning(text))
645 else:
646 self.interpreter.results.write(file=resfile, dir=path, force=force)
647
648
649 self.spin_display_params(pipe_name=pipe_name)
650
651
652 - def minimise_grid_search(self, inc=11, verbosity=0, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):
653 """Use value.set on all pipes."""
654
655
656 if model_from == None:
657 model_from = model
658 if analysis_from == None:
659 analysis_from = analysis
660
661
662 for method in methods:
663
664 self.set_self(key='method', value=method)
665
666
667 for glob_ini in list_glob_ini:
668
669 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method)
670
671
672 if not found:
673
674 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis_from, glob_ini=glob_ini, bundle=self.method)
675
676 if not found:
677
678 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini)
679
680
681 if not pipes.has_pipe(model_from_pipe_name):
682 model_from_pipe_name = self.name_pipe(method=self.method, model=MODEL_R2EFF, analysis='int', glob_ini=glob_ini)
683
684
685 if not pipes.has_pipe(model_from_pipe_name):
686 self.calc_r2eff(methods=[self.method], list_glob_ini=[glob_ini])
687
688
689 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method)
690 self.interpreter.pipe.switch(pipe_name)
691
692
693 self.interpreter.relax_disp.select_model(model)
694
695
696 if model not in [MODEL_R2EFF, MODEL_NOREX]:
697 self.interpreter.relax_disp.insignificance(level=self.insignificance)
698
699
700 subtitle(file=sys.stdout, text="Grid search for pipe='%s'" % (pipe_name), prespace=3)
701
702
703 if inc:
704 self.interpreter.minimise.grid_search(inc=inc, verbosity=verbosity, constraints=self.constraints, skip_preset=True)
705
706
707 else:
708
709 for param in MODEL_PARAMS[model]:
710 self.interpreter.value.set(param=param, index=None, force=False)
711
712
713 if is_r1_optimised(model=model):
714 self.interpreter.value.set(param='r1', index=None)
715
716
717
718 cdp.settings = self.settings
719
720
721 pipe_name = self.name_pipe(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini)
722 resfile = pipe_name.replace(" ", "_")
723 model_path = model.replace(" ", "_")
724 path = self.results_dir+sep+model_path
725
726 if found and not force:
727 file_path = get_file_path(file_name=resfile, dir=path)
728 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path)
729 warn(RelaxWarning(text))
730 else:
731 self.interpreter.results.write(file=resfile, dir=path, force=force)
732
733
734 self.spin_display_params(pipe_name=pipe_name)
735
736
737 - def cluster_spins(self, spin_id=None, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):
738 """Method to select spins for a pipe."""
739
740
741 if model_from == None:
742 model_from = model
743 if analysis_from == None:
744 analysis_from = analysis
745
746
747 for method in methods:
748
749 self.set_self(key='method', value=method)
750
751
752 for glob_ini in list_glob_ini:
753
754 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method)
755
756 if not found:
757
758 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini)
759
760
761 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method)
762 self.interpreter.pipe.switch(pipe_name)
763
764
765 subtitle(file=sys.stdout, text="Cluster spins '%s' for pipe='%s'" % (spin_id, pipe_name), prespace=3)
766
767
768 self.interpreter.relax_disp.cluster(cluster_id='sel', spin_id=spin_id)
769
770
771 cdp.settings = self.settings
772
773 if found and not force:
774 file_path = get_file_path(file_name=resfile, dir=path)
775 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path)
776 warn(RelaxWarning(text))
777 else:
778 self.interpreter.results.write(file=resfile, dir=path, force=force)
779
780
781 print("Clustered spins are:", cdp.clustering)
782
783
784 - def minimise_execute(self, verbosity=1, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False, mc_err_analysis=False):
785 """Use value.set on all pipes."""
786
787
788 if model_from == None:
789 model_from = model
790 if analysis_from == None:
791 analysis_from = analysis
792
793
794 for method in methods:
795
796 self.set_self(key='method', value=method)
797
798
799 for glob_ini in list_glob_ini:
800
801 found_pipe, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method)
802
803
804 if not found_pipe:
805
806 found_analysis, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis_from, glob_ini=glob_ini, bundle=self.method)
807
808 if not (found_pipe or found_analysis):
809
810 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini)
811
812
813 if not pipes.has_pipe(model_from_pipe_name):
814 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis='grid', glob_ini=glob_ini)
815
816
817 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method)
818 self.interpreter.pipe.switch(pipe_name)
819
820
821 self.interpreter.relax_disp.select_model(model)
822
823
824 subtitle(file=sys.stdout, text="Minimise for pipe='%s'" % (pipe_name), prespace=3)
825 if hasattr(cdp, "sim_number"):
826 subsection(file=sys.stdout, text="Performing Monte-Carlo minimisations on %i simulations"%(getattr(cdp, "sim_number")), prespace=0)
827
828
829 self.interpreter.minimise.execute(min_algor=self.min_algor, func_tol=self.opt_func_tol, max_iter=self.opt_max_iterations, constraints=self.constraints, scaling=True, verbosity=verbosity)
830
831
832 if mc_err_analysis:
833 self.interpreter.monte_carlo.error_analysis()
834
835
836 cdp.settings = self.settings
837
838
839 pipe_name = self.name_pipe(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini)
840 resfile = pipe_name.replace(" ", "_")
841 model_path = model.replace(" ", "_")
842 path = self.results_dir+sep+model_path
843
844 if found_pipe and not force:
845 file_path = get_file_path(file_name=resfile, dir=path)
846 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path)
847 warn(RelaxWarning(text))
848 else:
849 self.interpreter.results.write(file=resfile, dir=path, force=force)
850
851
852 self.spin_display_params(pipe_name=pipe_name)
853
854
855 - def name_pipe(self, method, model, analysis, glob_ini, clusterid=None):
856 """Generate a unique name for the data pipe."""
857
858
859
860 if clusterid == None:
861 clusterid = 'free spins'
862
863
864 name = "%s - %s - %s - %s - %s" % (method, model, analysis, glob_ini, clusterid)
865
866
867 name = name.replace(" ", "_")
868
869
870 return name
871
872
873 - def check_previous_result(self, method=None, model=None, analysis=None, glob_ini=None, clusterid=None, bundle=None):
874
875
876 found = False
877 if bundle == None:
878 bundle = self.pipe_bundle
879
880
881 pipe_name = self.name_pipe(method=method, model=model, analysis=analysis, glob_ini=glob_ini, clusterid=clusterid)
882
883
884 model_path = model.replace(" ", "_")
885 path = self.results_dir+sep+model_path
886
887
888 resfile = pipe_name.replace(" ", "_")
889
890
891 if pipes.has_pipe(pipe_name):
892 print("Detected the presence of previous '%s' pipe - switching to it." % pipe_name)
893 self.interpreter.pipe.switch(pipe_name)
894 found = True
895
896 else:
897
898 path1 = get_file_path(file_name=resfile, dir=path)
899 path2 = path1 + '.bz2'
900 path3 = path1 + '.gz'
901
902
903 if access(path1, F_OK) or access(path2, F_OK) or access(path2, F_OK):
904
905 print("Detected the presence of results files for the '%s' pipe - loading these instead of performing optimisation for a second time." % pipe_name)
906
907
908 self.interpreter.pipe.create(pipe_name=pipe_name, pipe_type=self.pipe_type, bundle=bundle)
909 self.interpreter.pipe.switch(pipe_name)
910
911
912 self.interpreter.results.read(file=resfile, dir=path)
913
914
915 found = True
916
917 return found, pipe_name, resfile, path
918
919
921 """Display parameters for model in pipe."""
922
923
924
925 if pipes.has_pipe(pipe_name):
926
927 if pipes.cdp_name() != pipe_name:
928 print("Detected the presence of previous '%s' pipe - switching to it." % pipe_name)
929 self.interpreter.pipe.switch(pipe_name)
930
931 else:
932
933 pipe_name_split = pipe_name.split("_-_")
934 method = pipe_name_split[0]
935 model = pipe_name_split[1]
936 analysis = pipe_name_split[2]
937 bundle = method
938
939 model_path = model.replace(" ", "_")
940 path = self.results_dir+sep+model_path
941
942 resfile = pipe_name.replace(" ", "_")
943
944
945 path1 = get_file_path(file_name=resfile, dir=path)
946 path2 = path1 + '.bz2'
947 path3 = path1 + '.gz'
948
949
950 if access(path1, F_OK) or access(path2, F_OK) or access(path2, F_OK):
951
952 print("Detected the presence of results files for the '%s' pipe - loading these instead of performing optimisation for a second time." % pipe_name)
953
954
955 self.interpreter.pipe.create(pipe_name=pipe_name, pipe_type=self.pipe_type, bundle=bundle)
956 self.interpreter.pipe.switch(pipe_name)
957
958
959 self.interpreter.results.read(file=resfile, dir=path)
960
961
962 found = True
963
964
965
966 my_dic = {}
967 spin_id_list = []
968
969
970 data = []
971
972 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=spin_id, full_info=True, return_id=True, skip_desel=True):
973
974 my_dic[spin_id] = {}
975
976 spin_id_list.append(spin_id)
977
978
979 cur_data_list = [repr(mol_name), repr(resi), repr(resn), repr(cur_spin.num), repr(cur_spin.name), spin_id]
980
981
982 params = cur_spin.params
983 my_dic[spin_id]['params'] = params
984 my_dic[spin_id]['params_err'] = []
985
986
987 for i, param in enumerate(params):
988
989 param_err = param + '_err'
990 my_dic[spin_id]['params_err'].append(param_err)
991
992
993 param_key_list = []
994 if param in PARAMS_R20:
995
996 for exp_type, frq, offset, ei, mi, oi, in loop_exp_frq_offset(return_indices=True):
997
998 param_key = generate_r20_key(exp_type=exp_type, frq=frq)
999 param_key_list.append(param_key)
1000 my_dic[spin_id]['param_key_list'] = param_key_list
1001 my_dic[spin_id][param] = {}
1002
1003
1004 if len(getattr(cur_spin, param)) == 0:
1005 param_val = None
1006 else:
1007 param_val = deepcopy(getattr(cur_spin, param)[param_key])
1008 my_dic[spin_id][param][param_key] = param_val
1009
1010
1011 if param_val == None:
1012 cur_data_list.append("%s"%param_val)
1013 else:
1014 cur_data_list.append("%3.3f"%param_val)
1015
1016 else:
1017
1018 param_val = deepcopy(getattr(cur_spin, param))
1019 my_dic[spin_id][param] = param_val
1020
1021
1022 if param_val == None:
1023 cur_data_list.append("%s"%param_val)
1024 else:
1025 cur_data_list.append("%.3f"%param_val)
1026
1027
1028 data.append(cur_data_list)
1029
1030
1031 cur_spin_id = spin_id_list[0]
1032 cur_spin_params = my_dic[cur_spin_id]['params']
1033 cur_param_keys = my_dic[cur_spin_id]['param_key_list']
1034
1035
1036 param_header = ["Molecule", "Res number", "Res name", "Spin number", "Spin name", "Spin id"]
1037
1038
1039 for param in cur_spin_params:
1040 if param in PARAMS_R20:
1041 for param_key in cur_param_keys:
1042
1043 cur_key = "%3.1f" % float(param_key.split()[-2])
1044 hstring = "%s_%s" % (param, cur_key)
1045 param_header.append(hstring)
1046 else:
1047 hstring = "%s" % (param)
1048 param_header.append(hstring)
1049
1050 write_data(out=sys.stdout, headings=param_header, data=data)
1051
1052
1054 """Method which return a list of tubles, where each tuble is a spectrum id and a list of spectrum ids which are replicated"""
1055
1056
1057 dublicates = [(val, [i for i in range(len(cpmg_frqs)) if cpmg_frqs[i] == val]) for val in cpmg_frqs]
1058
1059
1060 list_dub_mapping = []
1061 for i, dub in enumerate(dublicates):
1062
1063 cur_spectrum_id = spectrum_ids[i]
1064
1065
1066 cpmg_frq, list_index_occur = dub
1067
1068
1069 id_list = []
1070 if len(list_index_occur) > 1:
1071 for list_index in list_index_occur:
1072 id_list.append(spectrum_ids[list_index])
1073
1074
1075 list_dub_mapping.append((cur_spectrum_id, id_list))
1076
1077 return list_dub_mapping
1078
1079
1080 - def col_int(self, method=None, list_glob_ini=None, selection=None):
1081
1082
1083 res_dic = {}
1084 res_dic['method'] = method
1085 res_dic['selection'] = selection
1086 for glob_ini in list_glob_ini:
1087
1088 pipe_name = self.name_pipe(method=method, model='setup', analysis='int', glob_ini=glob_ini)
1089
1090
1091 if not pipes.has_pipe(pipe_name):
1092 self.set_int(methods=[method], list_glob_ini=[glob_ini])
1093
1094 if pipes.cdp_name() != pipe_name and pipes.has_pipe(pipe_name):
1095 self.interpreter.pipe.switch(pipe_name)
1096
1097 elif pipes.has_pipe(pipe_name) == False:
1098 continue
1099
1100
1101 res_dic[str(glob_ini)] = {}
1102 res_dic[str(glob_ini)]['peak_intensity'] = {}
1103 res_dic[str(glob_ini)]['peak_intensity_err'] = {}
1104 spin_point_peak_intensity_list = []
1105 spin_point_peak_intensity_err_list = []
1106
1107
1108 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=True):
1109
1110 res_dic[str(glob_ini)]['peak_intensity'][spin_id] = {}
1111 res_dic[str(glob_ini)]['peak_intensity_err'][spin_id] = {}
1112
1113
1114 for s_id in cdp.spectrum_ids:
1115
1116 if s_id in cur_spin.peak_intensity:
1117 peak_intensity_point = cur_spin.peak_intensity[s_id]
1118 peak_intensity_err_point = cur_spin.peak_intensity_err[s_id]
1119
1120 res_dic[str(glob_ini)]['peak_intensity'][spin_id][s_id] = peak_intensity_point
1121 res_dic[str(glob_ini)]['peak_intensity_err'][spin_id][s_id] = peak_intensity_err_point
1122 spin_point_peak_intensity_list.append(peak_intensity_point)
1123 spin_point_peak_intensity_err_list.append(peak_intensity_err_point)
1124
1125 res_dic[str(glob_ini)]['peak_intensity_arr'] = asarray(spin_point_peak_intensity_list)
1126 res_dic[str(glob_ini)]['peak_intensity_err_arr'] = asarray(spin_point_peak_intensity_err_list)
1127
1128 return res_dic
1129
1130
1131 - def plot_int_corr(self, corr_data, show=False, write_stats=False):
1132
1133
1134
1135 nr_cols = len(corr_data)
1136
1137 nr_rows = 4
1138
1139
1140 fig, axises = plt.subplots(nrows=nr_rows, ncols=nr_cols)
1141 fig.suptitle('Correlation plot')
1142
1143
1144
1145
1146
1147 data_dic = {}
1148
1149
1150 for i, row_axises in enumerate(axises):
1151
1152 for j, ax in enumerate(row_axises) :
1153
1154 data, methods, glob_inis = corr_data[j]
1155 data_x, data_y = data
1156 method_x, method_y = methods
1157 glob_ini_x, glob_ini_y = glob_inis
1158
1159 x = data_x[str(glob_ini_x)]['peak_intensity_arr']
1160 x_err = data_x[str(glob_ini_x)]['peak_intensity_err_arr']
1161 np = len(x)
1162
1163 y = data_y[str(glob_ini_y)]['peak_intensity_arr']
1164 y_err = data_y[str(glob_ini_y)]['peak_intensity_err_arr']
1165
1166
1167 if i == 0:
1168
1169 method_xy_NI = "int_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y)
1170 data_dic[method_xy_NI] = []
1171
1172 ax.plot(x, x, 'o', label='%s vs. %s' % (method_x, method_x))
1173 ax.plot(x, y, '.', label='%s vs. %s' % (method_y, method_x) )
1174
1175 np = len(y)
1176 ax.set_title(r'$I$' + ' for %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10)
1177 ax.legend(loc='upper center', shadow=True, prop=fontP)
1178 ax.ticklabel_format(style='sci', axis='x', scilimits=(0, 0))
1179 ax.ticklabel_format(style='sci', axis='y', scilimits=(0, 0))
1180 ax.set_xlabel(r'$I$')
1181 ax.set_ylabel(r'$I$')
1182
1183
1184
1185 a = sum(x * y) / sum(x**2)
1186 min_x = min(x)
1187 max_x = max(x)
1188 step_x = (max_x - min_x) / np
1189 x_arange = arange(min_x, max_x, step_x)
1190 y_arange = a * x_arange
1191
1192
1193 for k, x_k in enumerate(x):
1194 y_k = y[k]
1195 x_arange_k = x_arange[k]
1196 y_arange_k = y_arange[k]
1197 data_dic[method_xy_NI].append(["%3.5f"%x_k, "%3.5f"%y_k, "%3.5f"%x_arange_k, "%3.5f"%y_arange_k])
1198
1199
1200 if i == 1:
1201
1202 x_norm = x / x.max()
1203 y_norm = y / y.max()
1204
1205 ax.plot(x_norm, x_norm, 'o', label='%s vs. %s' % (method_x, method_x))
1206 ax.plot(x_norm, y_norm, '.', label='%s vs. %s' % (method_y, method_x))
1207
1208 np = len(y_norm)
1209 ax.set_title('Normalised intensity for %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10)
1210 ax.legend(loc='upper center', shadow=True, prop = fontP)
1211 ax.set_xlabel(r'$\mathrm{Norm.} I$')
1212 ax.set_ylabel(r'$\mathrm{Norm.} I$')
1213
1214
1215
1216 if i == 2:
1217
1218 method_xy_NI = "int_err_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y)
1219 data_dic[method_xy_NI] = []
1220
1221 ax.plot(x_err, x_err, 'o', label='%s vs. %s' % (method_x, method_x))
1222 ax.plot(x_err, y_err, '.', label='%s vs. %s' % (method_y, method_x))
1223
1224 np = len(y_err)
1225 ax.set_title(r'$\sigma(I)$' + ' for %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10)
1226 ax.legend(loc='upper center', shadow=True, prop = fontP)
1227 ax.set_xlabel(r'$\sigma(I)$')
1228 ax.set_ylabel(r'$\sigma(I)$')
1229
1230
1231
1232 a = sum(x_err * y_err) / sum(x_err**2)
1233 min_x = min(x_err)
1234 max_x = max(x_err)
1235 step_x = (max_x - min_x) / np
1236 x_err_arange = arange(min_x, max_x, step_x)
1237 y_err_arange = a * x_err_arange
1238
1239
1240 for k, x_err_k in enumerate(x_err):
1241 y_err_k = y_err[k]
1242 x_err_arange_k = x_err_arange[k]
1243 y_err_arange_k = y_err_arange[k]
1244 data_dic[method_xy_NI].append(["%3.5f"%x_err_k, "%3.5f"%y_err_k, "%3.5f"%x_err_arange_k, "%3.5f"%y_err_arange_k])
1245
1246
1247
1248 if i == 3:
1249
1250 x_to_x_err = x / x_err
1251 y_to_y_err = y / y_err
1252
1253 ax.plot(x_to_x_err, x_to_x_err, 'o', label='%s vs. %s' % (method_x, method_x))
1254 ax.plot(x_to_x_err, y_to_y_err, '.', label='%s vs. %s' % (method_y, method_x))
1255
1256 np = len(y_to_y_err)
1257 ax.set_title(r'$I/\sigma(I)$' + ' for %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10)
1258 ax.legend(loc='upper center', shadow=True, prop = fontP)
1259 ax.set_xlabel(r'$I/\sigma(I)$')
1260 ax.set_ylabel(r'$I/\sigma(I)$')
1261
1262 plt.tight_layout()
1263
1264
1265
1266 if write_stats:
1267
1268 headings_all = []
1269 method_xy_NI_all = []
1270
1271 for j in range(nr_cols):
1272 headings_j = []
1273 method_xy_NI_j = []
1274
1275 for i in range(nr_rows):
1276
1277 data, methods, glob_inis = corr_data[j]
1278 method_x, method_y = methods
1279 glob_ini_x, glob_ini_y = glob_inis
1280
1281
1282 if i == 0:
1283
1284 method_x_NI = "int_%s%s" % (method_x, glob_ini_x)
1285 method_y_NI = "int_%s%s" % (method_y, glob_ini_y)
1286 method_x_NI_lin = "int_lin_%s%s" % (method_x, glob_ini_x)
1287 method_y_NI_lin = "int_lin_%s%s" % (method_y, glob_ini_y)
1288 headings_j = headings_j + [method_x_NI, method_y_NI, method_x_NI_lin, method_y_NI_lin]
1289
1290 method_xy_NI = "int_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y)
1291 method_xy_NI_j.append(method_xy_NI)
1292
1293
1294 if i == 2:
1295
1296 method_x_NI = "int_err_%s%s" % (method_x, glob_ini_x)
1297 method_y_NI = "int_err_%s%s" % (method_y, glob_ini_y)
1298 method_x_NI_lin = "int_err_lin_%s%s" % (method_x, glob_ini_x)
1299 method_y_NI_lin = "int_err_lin_%s%s" % (method_y, glob_ini_y)
1300 headings_j = headings_j + [method_x_NI, method_y_NI, method_x_NI_lin, method_y_NI_lin]
1301
1302 method_xy_NI = "int_err_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y)
1303 method_xy_NI_j.append(method_xy_NI)
1304
1305 headings_all.append(headings_j)
1306 method_xy_NI_all.append(method_xy_NI_j)
1307
1308
1309 for j, headings_j in enumerate(headings_all):
1310 method_xy_NI_j = method_xy_NI_all[j]
1311
1312 data_w = []
1313 data_int = data_dic[method_xy_NI_j[0]]
1314 data_int_err = data_dic[method_xy_NI_j[1]]
1315
1316 for k, data_int_k in enumerate(data_int):
1317 data_int_err_k = data_int_err[k]
1318 data_w.append(data_int_k + data_int_err_k)
1319
1320
1321 data, methods, glob_inis = corr_data[j]
1322 data_x, data_y = data
1323 method_x, method_y = methods
1324 glob_ini_x, glob_ini_y = glob_inis
1325 np = len(data_int)
1326
1327
1328 selection = data_x['selection']
1329
1330 file_name_ini = 'int_corr_%s_%s_%s_%s_NP_%i' % (method_x, glob_ini_x, method_y, glob_ini_y, np)
1331 if selection == None:
1332 file_name_ini = file_name_ini + '_all'
1333 else:
1334 file_name_ini = file_name_ini + '_sel'
1335
1336 file_name = file_name_ini + '.txt'
1337 path = self.results_dir
1338
1339
1340
1341 png_file_name = file_name_ini + '.png'
1342 png_file_path = get_file_path(file_name=png_file_name, dir=path)
1343 plt.savefig(png_file_path, bbox_inches='tight')
1344
1345
1346 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True)
1347
1348
1349 write_data(out=file_obj, headings=headings_j, data=data_w)
1350
1351
1352 file_obj.close()
1353
1354 if show:
1355 plt.show()
1356
1357
1359
1360
1361 res_dic = {}
1362 for i, int_dic in enumerate(list_int_dics):
1363
1364 int_dic_ref = list_int_dics[0]
1365 method_ref = int_dic_ref['method']
1366 res_dic['method_ref'] = method_ref
1367 glob_ini_ref = list_glob_ini[0]
1368 res_dic['glob_ini_ref'] = str(glob_ini_ref)
1369 selection = int_dic_ref['selection']
1370 res_dic['selection'] = selection
1371
1372
1373 int_arr_ref = int_dic_ref[str(glob_ini_ref)]['peak_intensity_arr']
1374 res_dic['int_arr_ref'] = int_arr_ref
1375 int_err_arr_ref = int_dic_ref[str(glob_ini_ref)]['peak_intensity_err_arr']
1376 res_dic['int_err_arr_ref'] = int_err_arr_ref
1377
1378
1379 method_cur = int_dic['method']
1380 res_dic[method_cur] = {}
1381 res_dic[method_cur]['method'] = method_cur
1382 res_dic[method_cur]['sampling_sparseness'] = []
1383 res_dic[method_cur]['glob_ini'] = []
1384 res_dic[method_cur]['int_norm_std'] = []
1385
1386
1387 res_dic[method_cur]['r_xy_int'] = []
1388 res_dic[method_cur]['a_int'] = []
1389 res_dic[method_cur]['r_xy_int_err'] = []
1390 res_dic[method_cur]['a_int_err'] = []
1391
1392
1393 for glob_ini in list_glob_ini:
1394
1395 if str(glob_ini) not in int_dic:
1396 continue
1397
1398
1399 int_arr = int_dic[str(glob_ini)]['peak_intensity_arr']
1400 int_err_arr = int_dic[str(glob_ini)]['peak_intensity_err_arr']
1401
1402
1403
1404 if len(int_arr) != len(int_arr_ref):
1405 continue
1406
1407
1408 sampling_sparseness = float(glob_ini) / float(glob_ini_ref) * 100.
1409 res_dic[method_cur]['sampling_sparseness'].append(sampling_sparseness)
1410 res_dic[method_cur]['glob_ini'].append(glob_ini)
1411
1412
1413 res_dic[method_cur][str(glob_ini)] = {}
1414 res_dic[method_cur][str(glob_ini)]['sampling_sparseness'] = sampling_sparseness
1415 res_dic[method_cur][str(glob_ini)]['int_arr'] = int_arr
1416 res_dic[method_cur][str(glob_ini)]['int_err_arr'] = int_err_arr
1417
1418
1419
1420 x = int_arr_ref
1421 y = int_arr
1422
1423 a_int = sum(x*y) / sum(x**2)
1424 r_xy_int = sum(x*y) / sqrt(sum(x**2) * sum(y**2))
1425
1426 x = int_err_arr_ref
1427 y = int_err_arr
1428 a_int_err = sum(x*y) / sum(x**2)
1429 r_xy_int_err = sum(x*y) / sqrt(sum(x**2) * sum(y**2))
1430
1431 print(method_ref, method_cur, sampling_sparseness, glob_ini, r_xy_int**2, a_int, r_xy_int_err**2, a_int_err)
1432
1433
1434 res_dic[method_cur][str(glob_ini)]['r_xy_int'] = r_xy_int
1435 res_dic[method_cur]['r_xy_int'].append(r_xy_int)
1436 res_dic[method_cur][str(glob_ini)]['a_int'] = a_int
1437 res_dic[method_cur]['a_int'].append(a_int)
1438
1439 res_dic[method_cur][str(glob_ini)]['r_xy_int_err'] = r_xy_int_err
1440 res_dic[method_cur]['r_xy_int_err'].append(r_xy_int_err)
1441 res_dic[method_cur][str(glob_ini)]['a_int_err'] = a_int_err
1442 res_dic[method_cur]['a_int_err'].append(a_int_err)
1443
1444 res_dic[method_cur]['sampling_sparseness'] = asarray(res_dic[method_cur]['sampling_sparseness'])
1445 res_dic[method_cur]['glob_ini'] = asarray(res_dic[method_cur]['glob_ini'])
1446
1447 res_dic[method_cur]['r_xy_int'] = asarray(res_dic[method_cur]['r_xy_int'])
1448 res_dic[method_cur]['a_int'] = asarray(res_dic[method_cur]['a_int'])
1449 res_dic[method_cur]['r_xy_int_err'] = asarray(res_dic[method_cur]['r_xy_int_err'])
1450 res_dic[method_cur]['a_int_err'] = asarray(res_dic[method_cur]['a_int_err'])
1451
1452 return res_dic
1453
1454
1455 - def plot_int_stat(self, int_stat_dic=None, methods=[], list_glob_ini=[], show=False, write_stats=False):
1456
1457
1458 fig, axises = plt.subplots(nrows=2, ncols=1)
1459 fig.suptitle('Stats per NI')
1460 ax1, ax2 = axises
1461
1462
1463 min_a = 1.0
1464 max_a = 0.0
1465
1466 min_r_xy2 = 1.0
1467 max_r_xy2 = 0.0
1468
1469
1470 selection = int_stat_dic['selection']
1471
1472
1473 headings = []
1474 data_dic = {}
1475 data_dic_methods = []
1476 i_max = 0
1477
1478 for method in methods:
1479 if method not in int_stat_dic:
1480 continue
1481
1482
1483 NI = int_stat_dic[method]['glob_ini']
1484
1485 SS = int_stat_dic[method]['sampling_sparseness']
1486
1487
1488 headings = headings + ['method', 'SS', 'NI', 'slope_int', 'rxy2_int', 'slope_int_err', 'rxy2_int_err']
1489
1490
1491
1492 a_int = int_stat_dic[method]['a_int']
1493
1494 if max(a_int) > max_a:
1495 max_a = max(a_int)
1496 if min(a_int) < min_a:
1497 min_a = min(a_int)
1498
1499
1500 r_xy_int = int_stat_dic[method]['r_xy_int']
1501 r_xy_int2 = r_xy_int**2
1502
1503 if max(r_xy_int2) > max_r_xy2:
1504 max_r_xy2 = max(r_xy_int2)
1505 if min(r_xy_int2) < min_r_xy2:
1506 min_r_xy2 = min(r_xy_int2)
1507
1508
1509 a_int_err = int_stat_dic[method]['a_int_err']
1510 r_xy_int_err = int_stat_dic[method]['r_xy_int_err']
1511 r_xy_int_err2 = r_xy_int_err**2
1512
1513
1514 data_dic[method] = {}
1515 data_dic_methods.append(method)
1516 for i, NI_i in enumerate(NI):
1517 SS_i = SS[i]
1518 a_int_i = a_int[i]
1519 r_xy_int2_i = r_xy_int2[i]
1520 a_int_err_i = a_int_err[i]
1521 r_xy_int_err2_i = r_xy_int_err2[i]
1522 data_dic[method][str(i)] = ["%3.5f"%SS_i, "%i"%NI_i, "%3.5f"%a_int_i, "%3.5f"%r_xy_int2_i, "%3.5f"%a_int_err_i, "%3.5f"%r_xy_int_err2_i]
1523 if i > i_max:
1524 i_max = i
1525
1526 t = ax1.plot(SS, a_int, ".--", label='%s slope int'%method)
1527 color = t[0].get_color()
1528 ax1.plot(SS, a_int_err, ".-", label='%s slope int_err'%method, color=color)
1529
1530 t = ax2.plot(SS, r_xy_int2, "o--", label='%s r2 int'%method)
1531 color = t[0].get_color()
1532 ax2.plot(SS, r_xy_int_err2, "o-", label='%s r2 int_err'%method, color=color)
1533
1534
1535 data = []
1536
1537 for i in range(0, i_max+1):
1538 data_i = []
1539 for method in data_dic_methods:
1540 data_dic_m = data_dic[method]
1541
1542 if str(i) in data_dic_m:
1543 data_i = data_i + [method] + data_dic_m[str(i)]
1544 else:
1545 data_i = data_i + [method] + ["0", "0", "0", "0", "0", "0"]
1546
1547 data.append(data_i)
1548
1549
1550 ax1.legend(loc='lower left', shadow=True, prop = fontP)
1551
1552 ax1.set_xlabel('SS')
1553
1554 ax1.set_ylabel('Linear regression slope, without intercept')
1555
1556
1557 ax1.set_ylim(min_a*0.95, max_a*1.05)
1558 ax1.invert_xaxis()
1559
1560 ax2.legend(loc='lower right', shadow=True, prop = fontP)
1561 ax2.set_ylabel('Sample correlation ' + r'$r_{xy}^2$')
1562
1563
1564 ax2.set_ylim(min_r_xy2*0.95, max_r_xy2*1.05)
1565 ax2.invert_xaxis()
1566
1567
1568 if selection == None:
1569 file_name_ini = 'int_stat_all'
1570 else:
1571 file_name_ini = 'int_stat_sel'
1572
1573
1574 png_file_name = file_name_ini + '.png'
1575 png_file_path = get_file_path(file_name=png_file_name, dir=self.results_dir)
1576
1577
1578 if write_stats:
1579
1580 plt.savefig(png_file_path, bbox_inches='tight')
1581
1582 file_name = file_name_ini + '.txt'
1583 path = self.results_dir
1584 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True)
1585
1586
1587 write_data(out=file_obj, headings=headings, data=data)
1588
1589
1590 file_obj.close()
1591
1592
1593 if show:
1594 plt.show()
1595
1596
1597 - def col_r2eff(self, method=None, list_glob_ini=None, selection=None):
1598
1599
1600 res_dic = {}
1601 res_dic['method'] = method
1602 res_dic['selection'] = selection
1603 for glob_ini in list_glob_ini:
1604
1605 pipe_name = self.name_pipe(method=method, model=MODEL_R2EFF, analysis='int', glob_ini=glob_ini)
1606
1607
1608 if not pipes.has_pipe(pipe_name):
1609 self.calc_r2eff(methods=[method], list_glob_ini=[glob_ini])
1610
1611 if pipes.cdp_name() != pipe_name and pipes.has_pipe(pipe_name):
1612 self.interpreter.pipe.switch(pipe_name)
1613
1614 elif pipes.has_pipe(pipe_name) == False:
1615 continue
1616
1617
1618 res_dic[str(glob_ini)] = {}
1619 res_dic[str(glob_ini)]['r2eff'] = {}
1620 res_dic[str(glob_ini)]['r2eff_err'] = {}
1621 spin_point_r2eff_list = []
1622 spin_point_r2eff_err_list = []
1623
1624
1625 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=True):
1626
1627 res_dic[str(glob_ini)]['r2eff'][spin_id] = {}
1628 res_dic[str(glob_ini)]['r2eff_err'][spin_id] = {}
1629
1630
1631 for exp_type, frq, offset, point, ei, mi, oi, di in loop_exp_frq_offset_point(return_indices=True):
1632
1633 data_key = return_param_key_from_data(exp_type=self.exp_type, frq=frq, offset=offset, point=point)
1634
1635
1636 if data_key in cur_spin.r2eff:
1637 r2eff_point = cur_spin.r2eff[data_key]
1638 r2eff_err_point = cur_spin.r2eff_err[data_key]
1639
1640 res_dic[str(glob_ini)]['r2eff'][spin_id][data_key] = r2eff_point
1641 res_dic[str(glob_ini)]['r2eff_err'][spin_id][data_key] = r2eff_err_point
1642 spin_point_r2eff_list.append(r2eff_point)
1643 spin_point_r2eff_err_list.append(r2eff_err_point)
1644
1645 res_dic[str(glob_ini)]['r2eff_arr'] = asarray(spin_point_r2eff_list)
1646 res_dic[str(glob_ini)]['r2eff_err_arr'] = asarray(spin_point_r2eff_err_list)
1647
1648 return res_dic
1649
1650
1652
1653
1654
1655 nr_cols = len(corr_data)
1656
1657 nr_rows = 2
1658
1659
1660 fig, axises = plt.subplots(nrows=nr_rows, ncols=nr_cols)
1661 fig.suptitle('Correlation plot')
1662
1663
1664
1665
1666
1667 data_dic = {}
1668
1669
1670 for i, row_axises in enumerate(axises):
1671
1672 for j, ax in enumerate(row_axises):
1673
1674 data, methods, glob_inis = corr_data[j]
1675 data_x, data_y = data
1676 method_x, method_y = methods
1677 glob_ini_x, glob_ini_y = glob_inis
1678
1679 x = data_x[str(glob_ini_x)]['r2eff_arr']
1680 x_err = data_x[str(glob_ini_x)]['r2eff_err_arr']
1681 np = len(x)
1682
1683 y = data_y[str(glob_ini_y)]['r2eff_arr']
1684 y_err = data_y[str(glob_ini_y)]['r2eff_err_arr']
1685
1686
1687 if i == 0:
1688
1689 method_xy_NI = "r2eff_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y)
1690 data_dic[method_xy_NI] = []
1691
1692
1693
1694 a = sum(x * y) / sum(x**2)
1695 min_x = min(x)
1696 max_x = max(x)
1697 step_x = (max_x - min_x) / np
1698 x_arange = arange(min_x, max_x, step_x)
1699 y_arange = a * x_arange
1700
1701 ax.plot(x, x, 'o', label='%s vs. %s' % (method_x, method_x))
1702 ax.plot(x_arange, x_arange, 'b--')
1703 if len(x) != len(y):
1704 print(len(x), len(y))
1705 ax.plot(x, y, '.', label='%s vs. %s' % (method_y, method_x) )
1706 ax.plot(x_arange, y_arange, 'g--')
1707
1708 ax.set_title(r'$R_{2,\mathrm{eff}}$' + ' for %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10)
1709 ax.legend(loc='upper left', shadow=True, prop = fontP)
1710 ax.ticklabel_format(style='sci', axis='x', scilimits=(0, 0))
1711 ax.ticklabel_format(style='sci', axis='y', scilimits=(0, 0))
1712 ax.set_xlabel(r'$R_{2,\mathrm{eff}}$')
1713 ax.set_ylabel(r'$R_{2,\mathrm{eff}}$')
1714
1715
1716 for k, x_k in enumerate(x):
1717 y_k = y[k]
1718 x_arange_k = x_arange[k]
1719 y_arange_k = y_arange[k]
1720 data_dic[method_xy_NI].append(["%3.5f"%x_k, "%3.5f"%y_k, "%3.5f"%x_arange_k, "%3.5f"%y_arange_k])
1721
1722
1723 if i == 1:
1724
1725 method_xy_NI = "r2eff_to_err_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y)
1726 data_dic[method_xy_NI] = []
1727
1728 x_to_x_err = x / x_err
1729 y_to_y_err = y / y_err
1730
1731
1732
1733 a = sum(x_to_x_err * y_to_y_err) / sum(x_to_x_err**2)
1734 min_x = min(x_to_x_err)
1735 max_x = max(x_to_x_err)
1736 step_x = (max_x - min_x) / np
1737 x_to_x_err_arange = arange(min_x, max_x, step_x)
1738 y_to_x_err_arange = a * x_to_x_err_arange
1739
1740 ax.plot(x_to_x_err, x_to_x_err, 'o', label='%s vs. %s' % (method_x, method_x))
1741 ax.plot(x_to_x_err_arange, x_to_x_err_arange, 'b--')
1742 ax.plot(x_to_x_err, y_to_y_err, '.', label='%s vs. %s' % (method_y, method_x) )
1743 ax.plot(x_to_x_err_arange, y_to_x_err_arange, 'g--')
1744
1745 ax.set_title(r'$R_{2,\mathrm{eff}}/\sigma(R_{2,\mathrm{eff}})$' + ' for %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10)
1746 ax.legend(loc='upper left', shadow=True, prop = fontP)
1747 ax.set_xlabel(r'$R_{2,\mathrm{eff}}/\sigma(R_{2,\mathrm{eff}})$')
1748 ax.set_ylabel(r'$R_{2,\mathrm{eff}}/\sigma(R_{2,\mathrm{eff}})$')
1749
1750
1751 for k, x_to_x_err_k in enumerate(x_to_x_err):
1752 y_to_y_err_k = y_to_y_err[k]
1753 x_to_x_err_arange_k = x_to_x_err_arange[k]
1754 y_to_x_err_arange_k = y_to_x_err_arange[k]
1755 data_dic[method_xy_NI].append(["%3.5f"%x_to_x_err_k, "%3.5f"%y_to_y_err_k, "%3.5f"%x_to_x_err_arange_k, "%3.5f"%y_to_x_err_arange_k])
1756
1757
1758 plt.tight_layout()
1759
1760
1761
1762 if write_stats:
1763
1764 headings_all = []
1765 method_xy_NI_all = []
1766
1767 for j in range(nr_cols):
1768 headings_j = []
1769 method_xy_NI_j = []
1770
1771 for i in range(nr_rows):
1772
1773 data, methods, glob_inis = corr_data[j]
1774 method_x, method_y = methods
1775 glob_ini_x, glob_ini_y = glob_inis
1776
1777
1778 if i == 0:
1779
1780 method_x_NI = "r2eff_%s%s" % (method_x, glob_ini_x)
1781 method_y_NI = "r2eff_%s%s" % (method_y, glob_ini_y)
1782 method_x_NI_lin = "r2eff_lin_%s%s" % (method_x, glob_ini_x)
1783 method_y_NI_lin = "r2eff_lin_%s%s" % (method_y, glob_ini_y)
1784 headings_j = headings_j + [method_x_NI, method_y_NI, method_x_NI_lin, method_y_NI_lin]
1785
1786 method_xy_NI = "r2eff_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y)
1787 method_xy_NI_j.append(method_xy_NI)
1788
1789
1790 if i == 1:
1791
1792 method_x_NI = "r2eff_to_err_%s%s" % (method_x, glob_ini_x)
1793 method_y_NI = "r2eff_to_err_%s%s" % (method_y, glob_ini_y)
1794 method_x_NI_lin = "r2eff_to_err_lin_%s%s" % (method_x, glob_ini_x)
1795 method_y_NI_lin = "r2eff_to_err_lin_%s%s" % (method_y, glob_ini_y)
1796 headings_j = headings_j + [method_x_NI, method_y_NI, method_x_NI_lin, method_y_NI_lin]
1797
1798 method_xy_NI = "r2eff_to_err_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y)
1799 method_xy_NI_j.append(method_xy_NI)
1800
1801 headings_all.append(headings_j)
1802 method_xy_NI_all.append(method_xy_NI_j)
1803
1804
1805 for j, headings_j in enumerate(headings_all):
1806 method_xy_NI_j = method_xy_NI_all[j]
1807
1808 data_w = []
1809 data_r2eff = data_dic[method_xy_NI_j[0]]
1810 data_r2eff_to_err = data_dic[method_xy_NI_j[1]]
1811
1812 for k, data_r2eff_k in enumerate(data_r2eff):
1813 data_r2eff_to_err_k = data_r2eff_to_err[k]
1814 data_w.append(data_r2eff_k + data_r2eff_to_err_k)
1815
1816
1817 data, methods, glob_inis = corr_data[j]
1818 data_x, data_y = data
1819 method_x, method_y = methods
1820 glob_ini_x, glob_ini_y = glob_inis
1821 np = len(data_r2eff)
1822
1823
1824 selection = data_x['selection']
1825
1826 file_name_ini = 'r2eff_corr_%s_%s_%s_%s_NP_%i' % (method_x, glob_ini_x, method_y, glob_ini_y, np)
1827 if selection == None:
1828 file_name_ini = file_name_ini + '_all'
1829 else:
1830 file_name_ini = file_name_ini + '_sel'
1831
1832 file_name = file_name_ini + '.txt'
1833 path = self.results_dir
1834
1835
1836
1837 png_file_name = file_name_ini + '.png'
1838 png_file_path = get_file_path(file_name=png_file_name, dir=path)
1839 plt.savefig(png_file_path, bbox_inches='tight')
1840
1841
1842 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True)
1843
1844
1845 write_data(out=file_obj, headings=headings_j, data=data_w)
1846
1847
1848 file_obj.close()
1849
1850 if show:
1851 plt.show()
1852
1853
1855
1856
1857 res_dic = {}
1858 for i, r2eff_dic in enumerate(list_r2eff_dics):
1859
1860 r2eff_dic_ref = list_r2eff_dics[0]
1861 method_ref = r2eff_dic_ref['method']
1862 res_dic['method_ref'] = method_ref
1863 glob_ini_ref = list_glob_ini[0]
1864 res_dic['glob_ini_ref'] = str(glob_ini_ref)
1865 selection = r2eff_dic_ref['selection']
1866 res_dic['selection'] = selection
1867
1868
1869 r2eff_arr_ref = r2eff_dic_ref[str(glob_ini_ref)]['r2eff_arr']
1870 res_dic['r2eff_arr_ref'] = r2eff_arr_ref
1871 r2eff_err_arr_ref = r2eff_dic_ref[str(glob_ini_ref)]['r2eff_err_arr']
1872 res_dic['r2eff_err_arr_ref'] = r2eff_err_arr_ref
1873
1874
1875 method_cur = r2eff_dic['method']
1876 res_dic[method_cur] = {}
1877 res_dic[method_cur]['method'] = method_cur
1878 res_dic[method_cur]['sampling_sparseness'] = []
1879 res_dic[method_cur]['glob_ini'] = []
1880 res_dic[method_cur]['r2eff_norm_std'] = []
1881
1882
1883 res_dic[method_cur]['pearsons_correlation_coefficient'] = []
1884 res_dic[method_cur]['two_tailed_p_value'] = []
1885 res_dic[method_cur]['r_xy'] = []
1886 res_dic[method_cur]['a'] = []
1887 res_dic[method_cur]['r_xy_2'] = []
1888 res_dic[method_cur]['a_2'] = []
1889 res_dic[method_cur]['r_xy_int'] = []
1890 res_dic[method_cur]['a_int'] = []
1891 res_dic[method_cur]['b_int'] = []
1892
1893
1894 for glob_ini in list_glob_ini:
1895
1896 if str(glob_ini) not in r2eff_dic:
1897 continue
1898
1899
1900 r2eff_arr = r2eff_dic[str(glob_ini)]['r2eff_arr']
1901 r2eff_err_arr = r2eff_dic[str(glob_ini)]['r2eff_err_arr']
1902
1903
1904
1905
1906 if len(r2eff_arr) != len(r2eff_arr_ref):
1907 continue
1908
1909
1910 r2eff_norm_arr = r2eff_arr/r2eff_arr_ref
1911
1912
1913 r2eff_norm_std = std(r2eff_norm_arr, ddof=1)
1914
1915
1916 r2eff_diff_norm_arr = (r2eff_arr - r2eff_arr_ref) / r2eff_arr_ref
1917 r2eff_diff_norm_std = std(r2eff_diff_norm_arr, ddof=1)
1918
1919
1920 sampling_sparseness = float(glob_ini) / float(glob_ini_ref) * 100.
1921 res_dic[method_cur]['sampling_sparseness'].append(sampling_sparseness)
1922 res_dic[method_cur]['glob_ini'].append(glob_ini)
1923
1924
1925 res_dic[method_cur][str(glob_ini)] = {}
1926 res_dic[method_cur][str(glob_ini)]['sampling_sparseness'] = sampling_sparseness
1927 res_dic[method_cur][str(glob_ini)]['r2eff_arr'] = r2eff_arr
1928 res_dic[method_cur][str(glob_ini)]['r2eff_norm_arr'] = r2eff_norm_arr
1929 res_dic[method_cur][str(glob_ini)]['r2eff_norm_std'] = r2eff_norm_std
1930 res_dic[method_cur]['r2eff_norm_std'].append(r2eff_norm_std)
1931
1932 res_dic[method_cur][str(glob_ini)]['r2eff_diff_norm_arr'] = r2eff_diff_norm_arr
1933 res_dic[method_cur][str(glob_ini)]['r2eff_diff_norm_std'] = r2eff_diff_norm_std
1934
1935
1936
1937
1938 r2eff_vs_err_ref = r2eff_arr_ref / r2eff_err_arr_ref
1939 r2eff_vs_err = r2eff_arr / r2eff_err_arr
1940
1941
1942 pearsons_correlation_coefficient, two_tailed_p_value = pearsonr(r2eff_vs_err_ref, r2eff_vs_err)
1943
1944
1945
1946 x = r2eff_vs_err_ref
1947 x_m = mean(x)
1948 y = r2eff_vs_err
1949 y_m = mean(r2eff_vs_err)
1950
1951
1952 n = len(y)
1953 a_int = (sum(x*y) - 1./n * sum(x) * sum(y) ) / ( sum(x**2) - 1./n * (sum(x))**2 )
1954 b_int = 1./n * sum(y) - a_int * 1./n * sum(x)
1955
1956 r_xy_int = sum( (x - x_m)*(y - y_m) ) / sqrt( sum((x - x_m)**2) * sum((y - y_m)**2) )
1957
1958
1959 a = sum(x*y) / sum(x**2)
1960 r_xy = sum(x*y) / sqrt(sum(x**2) * sum(y**2))
1961
1962
1963 x_2 = r2eff_arr_ref
1964 y_2 = r2eff_arr
1965 a_2 = sum(x_2*y_2) / sum(x_2**2)
1966 r_xy_2 = sum(x_2*y_2) / sqrt(sum(x_2**2) * sum(y_2**2))
1967
1968
1969 print(method_ref, method_cur, sampling_sparseness, glob_ini, pearsons_correlation_coefficient, r_xy**2, a, r_xy_2**2, a_2)
1970
1971
1972 res_dic[method_cur][str(glob_ini)]['pearsons_correlation_coefficient'] = pearsons_correlation_coefficient
1973 res_dic[method_cur]['pearsons_correlation_coefficient'].append(pearsons_correlation_coefficient)
1974 res_dic[method_cur][str(glob_ini)]['two_tailed_p_value'] = two_tailed_p_value
1975 res_dic[method_cur]['two_tailed_p_value'].append(two_tailed_p_value)
1976 res_dic[method_cur][str(glob_ini)]['r_xy'] = r_xy
1977 res_dic[method_cur]['r_xy'].append(r_xy)
1978 res_dic[method_cur][str(glob_ini)]['a'] = a
1979 res_dic[method_cur]['a'].append(a)
1980
1981 res_dic[method_cur][str(glob_ini)]['r_xy_2'] = r_xy_2
1982 res_dic[method_cur]['r_xy_2'].append(r_xy_2)
1983 res_dic[method_cur][str(glob_ini)]['a_2'] = a_2
1984 res_dic[method_cur]['a_2'].append(a_2)
1985
1986 res_dic[method_cur][str(glob_ini)]['r_xy_int'] = r_xy_int
1987 res_dic[method_cur]['r_xy_int'].append(r_xy_int)
1988 res_dic[method_cur][str(glob_ini)]['a_int'] = a_int
1989 res_dic[method_cur]['a_int'].append(a_int)
1990 res_dic[method_cur][str(glob_ini)]['b_int'] = b_int
1991 res_dic[method_cur]['b_int'].append(b_int)
1992
1993 res_dic[method_cur]['sampling_sparseness'] = asarray(res_dic[method_cur]['sampling_sparseness'])
1994 res_dic[method_cur]['glob_ini'] = asarray(res_dic[method_cur]['glob_ini'])
1995 res_dic[method_cur]['r2eff_norm_std'] = asarray(res_dic[method_cur]['r2eff_norm_std'])
1996
1997 res_dic[method_cur]['pearsons_correlation_coefficient'] = asarray(res_dic[method_cur]['pearsons_correlation_coefficient'])
1998 res_dic[method_cur]['two_tailed_p_value'] = asarray(res_dic[method_cur]['two_tailed_p_value'])
1999 res_dic[method_cur]['r_xy'] = asarray(res_dic[method_cur]['r_xy'])
2000 res_dic[method_cur]['a'] = asarray(res_dic[method_cur]['a'])
2001 res_dic[method_cur]['r_xy_2'] = asarray(res_dic[method_cur]['r_xy_2'])
2002 res_dic[method_cur]['a_2'] = asarray(res_dic[method_cur]['a_2'])
2003 res_dic[method_cur]['r_xy_int'] = asarray(res_dic[method_cur]['r_xy_int'])
2004 res_dic[method_cur]['a_int'] = asarray(res_dic[method_cur]['a_int'])
2005 res_dic[method_cur]['b_int'] = asarray(res_dic[method_cur]['b_int'])
2006
2007 return res_dic
2008
2009
2010 - def plot_r2eff_stat(self, r2eff_stat_dic=None, methods=[], list_glob_ini=[], show=False, write_stats=False):
2011
2012
2013 fig, axises = plt.subplots(nrows=2, ncols=1)
2014 fig.suptitle('Stats per NI')
2015 ax1, ax2 = axises
2016
2017
2018 min_a = 1.0
2019 max_a = 0.0
2020
2021 min_r_xy2 = 1.0
2022 max_r_xy2 = 0.0
2023
2024
2025 selection = r2eff_stat_dic['selection']
2026
2027
2028 headings = []
2029 data_dic = {}
2030 data_dic_methods = []
2031 i_max = 0
2032
2033 for method in methods:
2034 if method not in r2eff_stat_dic:
2035 continue
2036
2037
2038 NI = r2eff_stat_dic[method]['glob_ini']
2039
2040 SS = r2eff_stat_dic[method]['sampling_sparseness']
2041
2042
2043 headings = headings + ['method', 'SS', 'NI', 'slope_r2eff', 'rxy2_r2eff', 'slope_r2eff_vs_err', 'rxy2_r2eff_vs_err']
2044
2045
2046
2047 a = r2eff_stat_dic[method]['a']
2048
2049 if max(a) > max_a:
2050 max_a = max(a)
2051 if min(a) < min_a:
2052 min_a = min(a)
2053
2054
2055 r_xy = r2eff_stat_dic[method]['r_xy']
2056 r_xy2 = r_xy**2
2057
2058 if max(r_xy2) > max_r_xy2:
2059 max_r_xy2 = max(r_xy2)
2060 if min(r_xy2) < min_r_xy2:
2061 min_r_xy2 = min(r_xy2)
2062
2063
2064 a_r2eff = r2eff_stat_dic[method]['a_2']
2065 r_xy_r2eff = r2eff_stat_dic[method]['r_xy_2']
2066 r_xy_r2eff2 = r_xy_r2eff**2
2067
2068
2069 data_dic[method] = {}
2070 data_dic_methods.append(method)
2071 for i, NI_i in enumerate(NI):
2072 SS_i = SS[i]
2073 a_i = a[i]
2074 r_xy2_i = r_xy2[i]
2075 a_r2eff_i = a_r2eff[i]
2076 r_xy_r2eff2_i = r_xy_r2eff2[i]
2077 data_dic[method][str(i)] = ["%3.5f"%SS_i, "%i"%NI_i, "%3.5f"%a_r2eff_i, "%3.5f"%r_xy_r2eff2_i, "%3.5f"%a_i, "%3.5f"%r_xy2_i]
2078 if i > i_max:
2079 i_max = i
2080
2081
2082
2083 t = ax1.plot(SS, a_r2eff, ".--", label='%s slope R2eff'%method)
2084 color = t[0].get_color()
2085 ax1.plot(SS, a, ".-", label='%s slope'%method, color=color)
2086
2087 t = ax2.plot(SS, r_xy_r2eff2, "o--", label='%s r2 R2eff'%method)
2088 color = t[0].get_color()
2089 ax2.plot(SS, r_xy2, "o-", label='%s r2'%method, color=color)
2090
2091
2092 data = []
2093
2094 for i in range(0, i_max+1):
2095 data_i = []
2096 for method in data_dic_methods:
2097 data_dic_m = data_dic[method]
2098
2099 if str(i) in data_dic_m:
2100 data_i = data_i + [method] + data_dic_m[str(i)]
2101 else:
2102 data_i = data_i + [method] + ["0", "0", "0", "0", "0", "0"]
2103
2104 data.append(data_i)
2105
2106
2107 ax1.legend(loc='lower left', shadow=True, prop = fontP)
2108
2109 ax1.set_xlabel('SS')
2110
2111 ax1.set_ylabel('Linear regression slope, without intercept')
2112
2113
2114 ax1.set_ylim(min_a*0.95, max_a*1.05)
2115 ax1.invert_xaxis()
2116
2117 ax2.legend(loc='lower right', shadow=True, prop = fontP)
2118 ax2.set_ylabel('Sample correlation ' + r'$r_{xy}^2$')
2119
2120
2121 ax2.set_ylim(min_r_xy2*0.95, max_r_xy2*1.05)
2122 ax2.invert_xaxis()
2123
2124
2125 if selection == None:
2126 file_name_ini = 'r2eff_stat_all'
2127 else:
2128 file_name_ini = 'r2eff_stat_sel'
2129
2130
2131 png_file_name = file_name_ini + '.png'
2132 png_file_path = get_file_path(file_name=png_file_name, dir=self.results_dir)
2133
2134
2135 if write_stats:
2136
2137 plt.savefig(png_file_path, bbox_inches='tight')
2138
2139 file_name = file_name_ini + '.txt'
2140 path = self.results_dir
2141 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True)
2142
2143
2144 write_data(out=file_obj, headings=headings, data=data)
2145
2146
2147 file_obj.close()
2148
2149
2150 if show:
2151 plt.show()
2152
2153
2154 - def col_min(self, method=None, model=None, analysis=None, list_glob_ini=None, selection=None):
2155
2156
2157 res_dic = {}
2158 res_dic['method'] = method
2159 res_dic['selection'] = selection
2160 res_dic['analysis'] = analysis
2161
2162 for glob_ini in list_glob_ini:
2163
2164 found, pipe_name, resfile, path = self.check_previous_result(method=method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=method)
2165
2166 if pipes.cdp_name() != pipe_name:
2167 self.interpreter.pipe.switch(pipe_name)
2168
2169
2170 res_dic[str(glob_ini)] = {}
2171 res_dic[str(glob_ini)]['params'] = {}
2172
2173
2174 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=True):
2175 params_list = cur_spin.params
2176
2177
2178 res_dic[str(glob_ini)]['params']['params_list'] = params_list
2179
2180
2181 for param in params_list:
2182
2183 res_dic[str(glob_ini)]['params'][param] = []
2184 res_dic[str(glob_ini)]['params'][param+'_resi'] = []
2185
2186
2187 res_dic[str(glob_ini)][param] = {}
2188
2189
2190 break
2191
2192
2193 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=True):
2194
2195 for param in params_list:
2196
2197
2198 res_dic[str(glob_ini)][param][spin_id] = {}
2199
2200
2201 param_key_list = []
2202 if param in PARAMS_R20:
2203
2204 for exp_type, frq, offset, ei, mi, oi, in loop_exp_frq_offset(return_indices=True):
2205
2206 param_key = generate_r20_key(exp_type=exp_type, frq=frq)
2207 param_key_list.append(param_key)
2208 res_dic[str(glob_ini)][param][spin_id]['param_key_list'] = param_key_list
2209
2210
2211 param_val = deepcopy(getattr(cur_spin, param)[param_key])
2212
2213 res_dic[str(glob_ini)]['params'][param].append(param_val)
2214 res_dic[str(glob_ini)][param][spin_id][param_key] = param_val
2215
2216 res_dic[str(glob_ini)]['params'][param+'_resi'].append("%i_%3.0f"%(resi, frq/1E6))
2217
2218 else:
2219
2220 param_val = deepcopy(getattr(cur_spin, param))
2221
2222 res_dic[str(glob_ini)]['params'][param].append(param_val)
2223 res_dic[str(glob_ini)][param][spin_id][param_key] = param_val
2224
2225 res_dic[str(glob_ini)]['params'][param+'_resi'].append("%i"%(resi))
2226
2227
2228 subtitle(file=sys.stdout, text="The minimised valus for '%s' model for pipe='%s at %s'" % (model, pipe_name, glob_ini), prespace=3)
2229
2230
2231 for param in params_list:
2232 res_dic[str(glob_ini)]['params'][param] = asarray(res_dic[str(glob_ini)]['params'][param])
2233 param_vals_str = " ".join(format(x, "2.3f") for x in res_dic[str(glob_ini)]['params'][param])
2234 print(param, param_vals_str)
2235
2236 return res_dic
2237
2238
2239 - def plot_min_corr(self, corr_data, show=False, write_stats=False):
2240
2241
2242 nr_cols = len(corr_data)
2243
2244 data_xy_0, methods_0, glob_inis_0 = corr_data[0]
2245 glob_ini_0 = glob_inis_0[0]
2246 params_list = data_xy_0[0][str(glob_ini_0)]['params']['params_list']
2247 nr_rows = len(params_list)
2248 analysis = data_xy_0[0]['analysis']
2249
2250
2251 fig, axises = plt.subplots(nrows=nr_rows, ncols=nr_cols)
2252 fig.suptitle('Correlation plot')
2253
2254
2255
2256
2257
2258 data_dic = {}
2259
2260
2261 for i, row_axises in enumerate(axises):
2262 param = params_list[i]
2263
2264
2265 for j, ax in enumerate(row_axises):
2266
2267 data, methods, glob_inis = corr_data[j]
2268 data_x, data_y = data
2269 method_x, method_y = methods
2270 glob_ini_x, glob_ini_y = glob_inis
2271
2272 x = data_x[str(glob_ini_x)]['params'][param]
2273 y = data_y[str(glob_ini_y)]['params'][param]
2274
2275 precision = abs(y-x) / ((x+y)/2)
2276
2277 precision_outlier = precision > 1.00
2278
2279 precision_outlier_nr = sum(precision_outlier)
2280 resis = data_x[str(glob_ini_x)]['params'][param+'_resi']
2281 resi_outlier_arr = asarray(resis)[precision_outlier]
2282 resi_outlier_arr_str = ":"+','.join(str(e) for e in resi_outlier_arr)
2283
2284 np = len(y)
2285
2286
2287 a = sum(x * y) / sum(x**2)
2288 min_xy = min(concatenate((x, y)))
2289 max_xy = max(concatenate((x, y)))
2290
2291 dx = (max_xy - min_xy) / np
2292 x_arange = arange(min_xy, max_xy + dx, dx)
2293 y_arange = a * x_arange
2294
2295 ax.plot(x, x, 'o', label='%s vs. %s' % (method_x, method_x))
2296 ax.plot(x, y, '.', label='%s vs. %s' % (method_y, method_x) )
2297
2298
2299 y_label = '%s'%param
2300
2301
2302 ax.set_ylabel(y_label)
2303
2304 ax.set_title('For %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10)
2305 ax.legend(loc='upper left', shadow=True, prop=fontP)
2306
2307
2308 if param == 'kex':
2309 ax.ticklabel_format(style='sci', axis='x', scilimits=(0, 0))
2310 ax.ticklabel_format(style='sci', axis='y', scilimits=(0, 0))
2311
2312
2313
2314
2315
2316
2317
2318
2319 ax.plot(x_arange, x_arange, 'b--')
2320 ax.plot(x_arange, y_arange, 'g--')
2321
2322
2323 method_xy_NI = "%s_%s_%s%s_%s%s" % (analysis, param, method_x, glob_ini_x, method_y, glob_ini_y)
2324 data_dic[method_xy_NI] = []
2325
2326
2327 for k, x_k in enumerate(x):
2328 y_k = y[k]
2329 x_arange_k = x_arange[k]
2330 y_arange_k = y_arange[k]
2331 precision_k = precision[k]
2332 resi_k = resis[k]
2333
2334 data_dic[method_xy_NI].append(["%3.5f"%x_k, "%3.5f"%y_k, "%3.5f"%x_arange_k, "%3.5f"%y_arange_k, "%3.5f"%precision_k, "%s"%resi_k, "%i"%precision_outlier_nr, "%s"%resi_outlier_arr_str])
2335
2336 plt.tight_layout()
2337
2338
2339
2340 if write_stats:
2341
2342 headings_all = []
2343 method_xy_NI_all = []
2344
2345 for j in range(nr_cols):
2346 headings_j = []
2347 method_xy_NI_j = []
2348
2349 for i in range(nr_rows):
2350
2351 data, methods, glob_inis = corr_data[j]
2352 method_x, method_y = methods
2353 glob_ini_x, glob_ini_y = glob_inis
2354 param = params_list[i]
2355
2356
2357 method_x_NI = "%s_%s_%s%s" % (analysis, param, method_x, glob_ini_x)
2358 method_y_NI = "%s_%s_%s%s" % (analysis, param, method_y, glob_ini_y)
2359 method_x_NI_lin = "%s_%s_lin_%s%s" % (analysis, param, method_x, glob_ini_x)
2360 method_y_NI_lin = "%s_%s_lin_%s%s" % (analysis, param, method_y, glob_ini_y)
2361 headings_j = headings_j + [method_x_NI, method_y_NI, method_x_NI_lin, method_y_NI_lin, "abs_diff_frac", "resi", "outlier_nr", "outlier_resi"]
2362
2363 method_xy_NI = "%s_%s_%s%s_%s%s" % (analysis, param, method_x, glob_ini_x, method_y, glob_ini_y)
2364 method_xy_NI_j.append(method_xy_NI)
2365
2366 headings_all.append(headings_j)
2367 method_xy_NI_all.append(method_xy_NI_j)
2368
2369
2370 for j, headings_j in enumerate(headings_all):
2371 method_xy_NI_j = method_xy_NI_all[j]
2372
2373 data_w = []
2374 method_xy_NI_r2 = method_xy_NI_j[0]
2375 data_r2 = data_dic[method_xy_NI_r2]
2376
2377
2378 for k, data_r2_k in enumerate(data_r2):
2379 data_row = data_r2_k
2380
2381 for method_xy_NI in method_xy_NI_j[1:]:
2382 data_param = data_dic[method_xy_NI]
2383
2384 try:
2385 data_param_row = data_param[k]
2386 except IndexError:
2387 data_param_row = len(data_param[0]) * ['0.0']
2388 data_param_row[-1] = ':'
2389
2390 data_row = data_row + data_param_row
2391
2392 data_w.append(data_row)
2393
2394
2395 data, methods, glob_inis = corr_data[j]
2396 data_x, data_y = data
2397 method_x, method_y = methods
2398 glob_ini_x, glob_ini_y = glob_inis
2399
2400
2401 selection = data_x['selection']
2402
2403 file_name_ini = '%s_%s_%s_%s_%s' % (analysis, method_x, glob_ini_x, method_y, glob_ini_y)
2404 if selection == None:
2405 file_name_ini = file_name_ini + '_all'
2406 else:
2407 file_name_ini = file_name_ini + '_sel'
2408
2409 file_name = file_name_ini + '.txt'
2410 path = self.results_dir
2411
2412
2413
2414 png_file_name = file_name_ini + '.png'
2415 png_file_path = get_file_path(file_name=png_file_name, dir=path)
2416 plt.savefig(png_file_path, bbox_inches='tight')
2417
2418
2419 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True)
2420
2421
2422 write_data(out=file_obj, headings=headings_j, data=data_w)
2423
2424
2425 file_obj.close()
2426
2427 if show:
2428 plt.show()
2429
2430
2432
2433
2434 res_dic = {}
2435 for i, min_dic in enumerate(list_r2eff_dics):
2436
2437 min_dic_ref = list_r2eff_dics[0]
2438 method_ref = min_dic_ref['method']
2439 res_dic['method_ref'] = method_ref
2440 glob_ini_ref = list_glob_ini[0]
2441 res_dic['glob_ini_ref'] = str(glob_ini_ref)
2442 selection = min_dic_ref['selection']
2443 res_dic['selection'] = selection
2444 params_list = min_dic_ref[str(glob_ini_ref)]['params']['params_list']
2445 res_dic['params_list'] = params_list
2446 analysis = min_dic_ref['analysis']
2447 res_dic['analysis'] = analysis
2448
2449
2450 method_cur = min_dic['method']
2451 res_dic[method_cur] = {}
2452
2453
2454 for j, param in enumerate(params_list):
2455 res_dic[param] = {}
2456
2457
2458 param_arr_ref = min_dic_ref[str(glob_ini_ref)]['params'][param]
2459 res_dic[param]['param_arr_ref'] = param_arr_ref
2460
2461 res_dic[method_cur][param] = {}
2462 res_dic[method_cur][param]['method'] = method_cur
2463 res_dic[method_cur][param]['sampling_sparseness'] = []
2464 res_dic[method_cur][param]['glob_ini'] = []
2465
2466
2467 res_dic[method_cur][param]['r_xy'] = []
2468 res_dic[method_cur][param]['a'] = []
2469 res_dic[method_cur][param]['precision_outlier_nr'] = []
2470 res_dic[method_cur][param]['resi_outlier'] = []
2471
2472
2473 for glob_ini in list_glob_ini:
2474
2475 if str(glob_ini) not in min_dic:
2476 continue
2477
2478
2479 param_arr = min_dic[str(glob_ini)]['params'][param]
2480 resis = min_dic[str(glob_ini)]['params'][param+'_resi']
2481
2482
2483
2484 if len(param_arr) != len(param_arr_ref):
2485 continue
2486
2487
2488 sampling_sparseness = float(glob_ini) / float(glob_ini_ref) * 100.
2489 res_dic[method_cur][param]['sampling_sparseness'].append(sampling_sparseness)
2490 res_dic[method_cur][param]['glob_ini'].append(glob_ini)
2491
2492
2493 res_dic[method_cur][param][str(glob_ini)] = {}
2494 res_dic[method_cur][param][str(glob_ini)]['sampling_sparseness'] = sampling_sparseness
2495 res_dic[method_cur][param][str(glob_ini)]['param_arr'] = param_arr
2496
2497
2498
2499 x = param_arr_ref
2500 x_m = mean(x)
2501 y = param_arr
2502 y_m = mean(y)
2503
2504
2505 a = sum(x*y) / sum(x**2)
2506 r_xy = sum(x*y) / sqrt(sum(x**2) * sum(y**2))
2507
2508
2509 precision = abs(y-x) / ((x+y)/2)
2510
2511 precision_outlier = precision > 1.00
2512
2513 precision_outlier_nr = sum(precision_outlier)
2514
2515 resi_outlier_arr = asarray(resis)[precision_outlier]
2516 resi_outlier_arr_str = ":"+','.join(str(e) for e in resi_outlier_arr)
2517
2518 print(param, method_ref, method_cur, sampling_sparseness, glob_ini, r_xy**2, a, precision_outlier_nr, resi_outlier_arr_str)
2519
2520
2521 res_dic[method_cur][param][str(glob_ini)]['r_xy'] = r_xy
2522 res_dic[method_cur][param]['r_xy'].append(r_xy)
2523 res_dic[method_cur][param][str(glob_ini)]['a'] = a
2524 res_dic[method_cur][param]['a'].append(a)
2525 res_dic[method_cur][param][str(glob_ini)]['precision_outlier_nr'] = precision_outlier_nr
2526 res_dic[method_cur][param]['precision_outlier_nr'].append(precision_outlier_nr)
2527 res_dic[method_cur][param]['resi_outlier'].append(resi_outlier_arr_str)
2528
2529 res_dic[method_cur][param]['sampling_sparseness'] = asarray(res_dic[method_cur][param]['sampling_sparseness'])
2530 res_dic[method_cur][param]['glob_ini'] = asarray(res_dic[method_cur][param]['glob_ini'])
2531
2532 res_dic[method_cur][param]['r_xy'] = asarray(res_dic[method_cur][param]['r_xy'])
2533 res_dic[method_cur][param]['a'] = asarray(res_dic[method_cur][param]['a'])
2534 res_dic[method_cur][param]['precision_outlier_nr'] = asarray(res_dic[method_cur][param]['precision_outlier_nr'])
2535 res_dic[method_cur][param]['resi_outlier'] = asarray(res_dic[method_cur][param]['resi_outlier'])
2536
2537 return res_dic
2538
2539
2540 - def plot_min_stat(self, min_stat_dic=None, methods=[], list_glob_ini=[], show=False, write_stats=False):
2541
2542
2543 min_a = 1.0
2544 max_a = 0.0
2545
2546 min_r_xy2 = 1.0
2547 max_r_xy2 = 0.0
2548
2549
2550 selection = min_stat_dic['selection']
2551 params_list = min_stat_dic['params_list']
2552 analysis = min_stat_dic['analysis']
2553
2554
2555 headings = []
2556 data_dic = {}
2557 data_dic_methods = []
2558
2559 i_max = 0
2560
2561 for method in methods:
2562 if method not in min_stat_dic:
2563 continue
2564
2565
2566 fig, axises = plt.subplots(nrows=len(params_list), ncols=1)
2567 fig.suptitle('Stats per NI %s' % method)
2568
2569
2570 data_dic[method] = {}
2571 data_dic_methods.append(method)
2572
2573 for j, param in enumerate(params_list):
2574 data_dic[method][param] = {}
2575
2576
2577 NI = min_stat_dic[method][param]['glob_ini']
2578
2579
2580 SS = min_stat_dic[method][param]['sampling_sparseness']
2581
2582
2583 headings = headings + ['method_%s'%param, 'SS', 'NI', 'slope', 'rxy2', 'outlier_nr', 'resi_outlier']
2584
2585
2586
2587 a = min_stat_dic[method][param]['a']
2588
2589 if max(a) > max_a:
2590 max_a = max(a)
2591 if min(a) < min_a:
2592 min_a = min(a)
2593
2594
2595 r_xy = min_stat_dic[method][param]['r_xy']
2596 r_xy2 = r_xy**2
2597
2598 if max(r_xy2) > max_r_xy2:
2599 max_r_xy2 = max(r_xy2)
2600 if min(r_xy2) < min_r_xy2:
2601 min_r_xy2 = min(r_xy2)
2602
2603
2604 precision_outlier_nr = min_stat_dic[method][param]['precision_outlier_nr']
2605 resi_outlier = min_stat_dic[method][param]['resi_outlier']
2606
2607
2608 for i, NI_i in enumerate(NI):
2609 SS_i = SS[i]
2610 a_i = a[i]
2611 r_xy2_i = r_xy2[i]
2612 precision_outlier_nr_i = precision_outlier_nr[i]
2613 resi_outlier_i = resi_outlier[i]
2614 data_dic[method][param][str(i)] = ["%3.5f"%SS_i, "%i"%NI_i, "%3.5f"%a_i, "%3.5f"%r_xy2_i, "%i"%precision_outlier_nr_i, "%s"%resi_outlier_i]
2615 if i > i_max:
2616 i_max = i
2617
2618 ax = axises[j]
2619 ax.plot(SS, a, ".-", label='%s_%s_a' % (method, param) )
2620 ax.plot(SS, r_xy2, "o--", label='%s_%s_r_xy2' % (method, param) )
2621 ax.legend(loc='lower left', shadow=True, prop = fontP)
2622 ax.set_xlabel('SS')
2623 ax.invert_xaxis()
2624
2625
2626
2627
2628 data = []
2629
2630
2631 for i in range(0, i_max+1):
2632 data_i = []
2633 for method in data_dic_methods:
2634 data_dic_m = data_dic[method]
2635
2636 for j, param in enumerate(params_list):
2637
2638 if str(i) in data_dic_m[param]:
2639 data_i = data_i + ["%s_%s" % (method, param)] + data_dic_m[param][str(i)]
2640 else:
2641 data_i = data_i + ["%s_%s" % (method, param)] + ["0", "0", "0", "0", "0", ":"]
2642
2643 data.append(data_i)
2644
2645
2646 if selection == None:
2647 file_name_ini = '%s_stat_all' % analysis
2648 else:
2649 file_name_ini = '%s_stat_sel' % analysis
2650
2651
2652 png_file_name = file_name_ini + '.png'
2653 png_file_path = get_file_path(file_name=png_file_name, dir=self.results_dir)
2654
2655
2656 if write_stats:
2657
2658 plt.savefig(png_file_path, bbox_inches='tight')
2659
2660 file_name = file_name_ini + '.txt'
2661 path = self.results_dir
2662 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True)
2663
2664
2665 write_data(out=file_obj, headings=headings, data=data)
2666
2667
2668 file_obj.close()
2669
2670
2671 if show:
2672 plt.show()
2673
2674
2675 - def write_results(self, method=None, model=None, analysis=None, list_glob_ini=None, selection=None, write_disp=True):
2676
2677 for glob_ini in list_glob_ini:
2678
2679 found, pipe_name, resfile, path = self.check_previous_result(method=method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=method)
2680
2681 if pipes.cdp_name() != pipe_name:
2682 self.interpreter.pipe.switch(pipe_name)
2683
2684
2685 section(file=sys.stdout, text="Results writing for pipe='%s"%(pipe_name), prespace=2, postspace=0)
2686 model_params = deepcopy(MODEL_PARAMS[model])
2687 subsection(file=sys.stdout, text="Model %s, with params='%s"%(model, model_params), prespace=0)
2688
2689
2690 model_path = model.replace(" ", "_")
2691 analysis_path = analysis.replace(" ", "_")
2692 path = self.results_dir+sep+model_path+sep+analysis_path
2693
2694
2695 if write_disp:
2696 path_disp = path+sep+"disp_curves"+sep+method+sep+str(glob_ini)
2697 self.interpreter.relax_disp.plot_disp_curves(dir=path_disp, force=True)
2698 self.interpreter.relax_disp.write_disp_curves(dir=path_disp, force=True)
2699
2700
2701 self.interpreter.value.write(param='model', file='model.out', dir=path, force=True)
2702
2703 models_tested = None
2704
2705
2706 filep = str(glob_ini)+"_"+method+"_"
2707 path_par = path+sep+"r2"
2708 if has_cpmg_exp_type():
2709
2710 self.write_results_test(path=path_par, model=model, models_tested=models_tested, param='r2', file_name_ini=filep+'r20')
2711
2712
2713 self.write_results_test(path=path_par, model=model, models_tested=models_tested, param='r2a', file_name_ini=filep+'r20a')
2714 self.write_results_test(path=path_par, model=model, models_tested=models_tested, param='r2b', file_name_ini=filep+'r20b')
2715
2716
2717 path_par = path+sep+"pop"
2718 search = method+"_"+"pA"
2719 self.write_results_test(path=path_par, model=model, models_tested=models_tested, search=search, param='pA', file_name_ini=filep+'pA')
2720 self.write_results_test(path=path_par, model=model, models_tested=models_tested, param='pB', file_name_ini=filep+'pB')
2721
2722
2723 path_par = path+sep+"dw"
2724 search = method+"_"+"dw"
2725 self.write_results_test(path=path_par, model=model, models_tested=models_tested, search=search, param='dw', file_name_ini=filep+'dw')
2726
2727
2728 path_par = path+sep+"rate"
2729 params = ['k_AB', 'kex', 'tex']
2730 for param in params:
2731 search = method+"_"+param
2732 self.write_results_test(path=path_par, model=model, models_tested=models_tested, search=search, param=param, file_name_ini=filep+param)
2733
2734
2735 if not (model == MODEL_R2EFF and has_fixed_time_exp_type()):
2736 path_par = path+sep+"chi2"
2737 self.interpreter.value.write(param='chi2', file=filep+'chi2.out', dir=path_par, force=True)
2738 search = method+"_"+"chi2"
2739 col_file_name="collect_%s.sh"%search
2740 self.write_convert_file(file_name=col_file_name, path=path_par, search=search)
2741
2742 self.interpreter.grace.write(y_data_type='chi2', file='chi2.agr', dir=path_par+sep+"grace", force=True)
2743
2744
2745 - def write_results_test(self, path=None, model=None, models_tested=None, search=None, param=None, file_name_ini=None):
2746 """Create a set of results, text and Grace files for the current data pipe.
2747
2748 @keyword path: The directory to place the files into.
2749 @type path: str
2750 @keyword model: The model tested.
2751 @type model: None or str
2752 @keyword model_tested: List of models tested, if the pipe is final.
2753 @type model_tested: None or list of str.
2754 @keyword param: The param to write out.
2755 @type param: None or list of str.
2756 @keyword file_name_ini: The initial part of the file name for the grace and text files.
2757 @type file_name_ini: None or str.
2758 """
2759
2760
2761 if file_name_ini == None:
2762 file_name_ini = param
2763
2764
2765 write_result = False
2766 if model != None:
2767
2768 model_params = deepcopy(MODEL_PARAMS[model])
2769
2770 if param in model_params:
2771 write_result = True
2772
2773
2774 elif model == None:
2775
2776 for model_tested in models_tested:
2777
2778 model_params = deepcopy(MODEL_PARAMS[model_tested])
2779
2780 if param in model_params:
2781 write_result = True
2782 break
2783
2784
2785 if write_result:
2786 self.interpreter.value.write(param=param, file='%s.out'%file_name_ini, dir=path, force=True)
2787
2788 if search != None:
2789 col_file_name="collect_%s.sh"%search
2790 self.write_convert_file(file_name=col_file_name, path=path, search=search)
2791
2792
2793 self.interpreter.grace.write(x_data_type='res_num', y_data_type=param, file='%s.agr'%file_name_ini, dir=path+sep+"grace", force=True)
2794
2796 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True)
2797
2798
2799 file_obj.write('#! /bin/bash' + '\n')
2800 file_obj.write('SEARCH=%s'%(search) + '\n')
2801 file_obj.write('FILES=(*_${SEARCH}.out)' + '\n')
2802 file_obj.write('readarray -t FILESSORT < <(for a in "${FILES[@]}"; do echo "$a"; done | sort -Vr)' + '\n')
2803 file_obj.write('# Skip the first two lines of header' + '\n')
2804 file_obj.write("tail -n+3 ${FILESSORT[0]} | sed 's,^# ,,' | grep -v 'None None' | awk '{print $2,$3,$5}' | column -t > collect_${SEARCH}.tmp" + '\n')
2805 file_obj.write('# Make array' + '\n')
2806 file_obj.write('ACUT=(collect_${SEARCH}.tmp)' + '\n')
2807 file_obj.write('for f in "${FILESSORT[@]}"; do' + '\n')
2808 file_obj.write(' FNAME="${f%.*}"' + '\n')
2809 file_obj.write(' NI=`echo $f | cut -d"_" -f1`' + '\n')
2810 file_obj.write(' echo "Processing $f with NI=$NI"' + '\n')
2811 file_obj.write(' tail -n+3 $f | sed "s,^# ,," | grep -v "None None" | sed "s,value,${NI}," | sed "s,error,${NI}," | awk %s{print $6,$7}%s | column -t > ${FNAME}.tmp'%("'","'") + '\n')
2812 file_obj.write(' ACUT+=(${FNAME}.tmp)' + '\n')
2813 file_obj.write('done' + '\n')
2814 file_obj.write('#paste "${ACUT[@]}" | column -t > collect_${SEARCH}.txt' + '\n')
2815 file_obj.write('paste "${ACUT[@]}" > collect_${SEARCH}.txt' + '\n')
2816 file_obj.write('rm ${ACUT[@]}' + '\n')
2817
2818
2819
2820 file_obj.close()
2821
2822 chmod(file_path, S_IRWXU|S_IRGRP|S_IROTH)
2823
2824
2825 - def create_mc_data(self, number=500, distribution="measured", fixed_error=None, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):
2826 """Create MC data."""
2827
2828
2829 if model_from == None:
2830 model_from = model
2831 if analysis_from == None:
2832 analysis_from = analysis
2833
2834
2835 for method in methods:
2836
2837 self.set_self(key='method', value=method)
2838
2839
2840 for glob_ini in list_glob_ini:
2841
2842 found_pipe, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method)
2843
2844
2845 if not found_pipe:
2846
2847 found_analysis, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis_from, glob_ini=glob_ini, bundle=self.method)
2848
2849
2850 subtitle(file=sys.stdout, text="MC data for pipe='%s'" % (pipe_name), prespace=3)
2851
2852
2853 self.interpreter.relax_disp.select_model(model)
2854
2855
2856 self.interpreter.monte_carlo.setup(number=number)
2857 self.interpreter.monte_carlo.create_data(distribution=distribution, fixed_error=fixed_error)
2858 self.interpreter.monte_carlo.initial_values()
2859
2860
2861 cdp.settings = self.settings
2862
2863
2864 pipe_name = self.name_pipe(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini)
2865 resfile = pipe_name.replace(" ", "_")
2866 model_path = model.replace(" ", "_")
2867 path = self.results_dir+sep+model_path
2868
2869 if found_pipe and not force:
2870 file_path = get_file_path(file_name=resfile, dir=path)
2871 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path)
2872 warn(RelaxWarning(text))
2873 else:
2874 self.interpreter.results.write(file=resfile, dir=path, force=force)
2875
2876
2877 - def summarize_results(self, methods=None, model=None, analysis=None, list_glob_ini=None, selections=None):
2878 """ Collect all results and make statistics"""
2879
2880 ref_method, ni_method = methods
2881 ref_ni, ni_list = list_glob_ini
2882
2883 def collect_data(pipe_name=None, selections=None, ni=None):
2884
2885 self.interpreter.pipe.switch(pipe_name)
2886
2887
2888 all_intensities = []
2889 all_intensities_err = []
2890
2891 for selection in selections:
2892 sel_intensities = []
2893 sel_intensities_err = []
2894 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=False):
2895
2896
2897 for s_id in cdp.spectrum_ids:
2898
2899 if s_id in cur_spin.peak_intensity:
2900 peak_intensity_point = cur_spin.peak_intensity[s_id]
2901 peak_intensity_err_point = cur_spin.peak_intensity_err[s_id]
2902
2903 sel_intensities.append(peak_intensity_point)
2904 sel_intensities_err.append(peak_intensity_err_point)
2905
2906
2907 all_intensities.append(sel_intensities)
2908 all_intensities_err.append(sel_intensities_err)
2909
2910
2911 all_r2eff = []
2912 all_r2eff_err = []
2913
2914 for selection in selections:
2915 sel_r2eff = []
2916 sel_r2eff_err = []
2917
2918
2919 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=False):
2920
2921 for exp_type, frq, offset, point, ei, mi, oi, di in loop_exp_frq_offset_point(return_indices=True):
2922
2923 data_key = return_param_key_from_data(exp_type=self.exp_type, frq=frq, offset=offset, point=point)
2924
2925
2926 if data_key in cur_spin.r2eff:
2927 r2eff_point = cur_spin.r2eff[data_key]
2928 r2eff_err_point = cur_spin.r2eff_err[data_key]
2929
2930 sel_r2eff.append(r2eff_point)
2931 sel_r2eff_err.append(r2eff_err_point)
2932
2933
2934 all_r2eff.append(sel_r2eff)
2935 all_r2eff_err.append(sel_r2eff_err)
2936
2937
2938 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=True):
2939 params_list = cur_spin.params
2940
2941
2942 break
2943
2944
2945 all_pars = []
2946 all_pars_err = []
2947
2948
2949 for selection in selections:
2950 sel_pars = []
2951 sel_pars_err = []
2952
2953
2954 for param in params_list:
2955
2956 sel_pars_list = []
2957 sel_pars_list_err = []
2958
2959
2960 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=True):
2961 if param in PARAMS_R20:
2962
2963 for exp_type, frq, offset, ei, mi, oi, in loop_exp_frq_offset(return_indices=True):
2964
2965 param_key = generate_r20_key(exp_type=exp_type, frq=frq)
2966
2967
2968 param_val = deepcopy(getattr(cur_spin, param)[param_key])
2969 param_err = deepcopy(getattr(cur_spin, param+"_err")[param_key])
2970
2971 else:
2972
2973 param_val = deepcopy(getattr(cur_spin, param))
2974 param_err = deepcopy(getattr(cur_spin, param+"_err"))
2975
2976
2977 sel_pars_list.append(param_val)
2978 sel_pars_list_err.append(param_err)
2979
2980
2981 sel_pars.append(sel_pars_list)
2982 sel_pars_err.append(sel_pars_list_err)
2983
2984
2985 all_pars.append(sel_pars)
2986 all_pars_err.append(sel_pars_err)
2987
2988
2989 data_dic[pipe_name] = {}
2990 data_dic[pipe_name]['ni'] = ni
2991 data_dic[pipe_name]['all_intensities'] = all_intensities
2992 data_dic[pipe_name]['all_intensities_err'] = all_intensities_err
2993 data_dic[pipe_name]['all_r2eff'] = all_r2eff
2994 data_dic[pipe_name]['all_r2eff_err'] = all_r2eff_err
2995 data_dic['params_list'] = params_list
2996 data_dic[pipe_name]['all_pars'] = all_pars
2997 data_dic[pipe_name]['all_pars_err'] = all_pars_err
2998
2999
3000 found, ref_pipe_name, resfile, path = self.check_previous_result(method=ref_method, model=model, analysis=analysis, glob_ini=ref_ni, bundle=ref_method)
3001
3002
3003 data_dic = {}
3004
3005 collect_data(pipe_name=ref_pipe_name, selections=selections, ni=ref_ni)
3006
3007
3008
3009 ni_pipe_names = []
3010 for i, ni in enumerate(ni_list):
3011
3012 found, ni_pipe_name, resfile, path = self.check_previous_result(method=ni_method, model=model, analysis=analysis, glob_ini=ni, bundle=ref_method)
3013 ni_pipe_names.append(ni_pipe_name)
3014
3015
3016 collect_data(pipe_name=ni_pipe_name, selections=selections, ni=ni)
3017
3018
3019 self.interpreter.pipe.switch(ref_pipe_name)
3020
3021
3022 if ni_pipe_name != ref_pipe_name:
3023 self.interpreter.pipe.delete(pipe_name=ni_pipe_name)
3024
3025
3026 def get_stats(x=None, y=None, x_err=None, y_err=None):
3027
3028 x = asarray(x)
3029 y = asarray(y)
3030
3031
3032 g = x/x
3033 h = y/x
3034
3035
3036 d_xy = x - y
3037 d_gh = g - h
3038
3039
3040 mean_d_xy = mean(d_xy)
3041 mean_d_gh = mean(d_gh)
3042
3043
3044 std_d_xy = std(d_xy, ddof=1)
3045 std_d_gh = std(d_gh, ddof=1)
3046
3047
3048 rmsd_xy = sqrt(mean(square(d_xy)))
3049 rmsd_gh = sqrt(mean(square(d_gh)))
3050
3051 if x_err != None and y_err != None:
3052 pool_std_xy = sqrt( mean(square(x_err) + square(y_err)) )
3053
3054 g_err = x_err / x
3055 h_err = y_err / x
3056
3057 pool_std_gh= sqrt( mean(square(g_err) + square(h_err)) )
3058
3059 return mean_d_xy, mean_d_gh, std_d_xy, std_d_gh, pool_std_xy, pool_std_gh, rmsd_xy, rmsd_gh
3060
3061 else:
3062 return mean_d_xy, mean_d_gh, std_d_xy, std_d_gh, rmsd_xy, rmsd_gh
3063
3064
3065 headers = ["ni", "pct"]
3066 types = ["int", "r2eff"] + data_dic['params_list']
3067
3068
3069 all_data = []
3070
3071 for i, ni in enumerate(ni_list):
3072
3073 pipe_name = ni_pipe_names[i]
3074
3075
3076 d_row = []
3077 pct = float(ni)/float(ref_ni)*100
3078
3079
3080 d_row.append(ni)
3081 d_row.append(pct)
3082
3083
3084 for j, selection in enumerate(selections):
3085
3086 k = 0
3087 for type_i in types:
3088 if type_i == 'int':
3089
3090 if i == 0:
3091 headers += ["%s_%s_mean_d_xy"%(j, type_i), "%s_%s_mean_d_gh"%(j, type_i), "%s_%s_std_d_xy"%(j, type_i), "%s_%s_std_d_gh"%(j, type_i), "%s_%s_pool_d_xy"%(j, type_i), "%s_%s_pool_d_gh"%(j, type_i), "%s_%s_rmsd_d_xy"%(j, type_i), "%s_%s_rmsd_d_gh"%(j, type_i)]
3092 x = data_dic[ref_pipe_name]['all_intensities'][j]
3093 y = data_dic[pipe_name]['all_intensities'][j]
3094 x_err = data_dic[ref_pipe_name]['all_intensities_err'][j]
3095 y_err = data_dic[pipe_name]['all_intensities_err'][j]
3096
3097 mean_d_xy, mean_d_gh, std_d_xy, std_d_gh, pool_std_xy, pool_std_gh, rmsd_xy, rmsd_gh = get_stats(x=x, y=y, x_err=x_err, y_err=y_err)
3098 d_row += [mean_d_xy, mean_d_gh, std_d_xy, std_d_gh, pool_std_xy, pool_std_gh, rmsd_xy, rmsd_gh]
3099
3100 elif type_i == 'r2eff':
3101
3102 if i == 0:
3103 headers += ["%s_%s_mean_d_xy"%(j, type_i), "%s_%s_mean_d_gh"%(j, type_i), "%s_%s_std_d_xy"%(j, type_i), "%s_%s_std_d_gh"%(j, type_i), "%s_%s_pool_d_xy"%(j, type_i), "%s_%s_pool_d_gh"%(j, type_i), "%s_%s_rmsd_d_xy"%(j, type_i), "%s_%s_rmsd_d_gh"%(j, type_i)]
3104 x = data_dic[ref_pipe_name]['all_r2eff'][j]
3105 y = data_dic[pipe_name]['all_r2eff'][j]
3106 x_err = data_dic[ref_pipe_name]['all_r2eff_err'][j]
3107 y_err = data_dic[pipe_name]['all_r2eff_err'][j]
3108
3109 if len(x) != len(y):
3110 d_row += [9999.0, 9999.0, 9999.0, 9999.0, 9999.0, 9999.0, 9999.0, 9999.0]
3111 else:
3112 mean_d_xy, mean_d_gh, std_d_xy, std_d_gh, pool_std_xy, pool_std_gh, rmsd_xy, rmsd_gh = get_stats(x=x, y=y, x_err=x_err, y_err=y_err)
3113 d_row += [mean_d_xy, mean_d_gh, std_d_xy, std_d_gh, pool_std_xy, pool_std_gh, rmsd_xy, rmsd_gh]
3114
3115 elif type_i in data_dic['params_list']:
3116
3117
3118 if i == 0:
3119 headers += ["%s_%s_mean_d_xy"%(j, type_i), "%s_%s_mean_d_gh"%(j, type_i), "%s_%s_std_d_xy"%(j, type_i), "%s_%s_std_d_gh"%(j, type_i), "%s_%s_pool_d_xy"%(j, type_i), "%s_%s_pool_d_gh"%(j, type_i), "%s_%s_rmsd_d_xy"%(j, type_i), "%s_%s_rmsd_d_gh"%(j, type_i)]
3120 x = data_dic[ref_pipe_name]['all_pars'][j][k]
3121 y = data_dic[pipe_name]['all_pars'][j][k]
3122 x_err = data_dic[ref_pipe_name]['all_pars_err'][j][k]
3123 y_err = data_dic[pipe_name]['all_pars_err'][j][k]
3124
3125 mean_d_xy, mean_d_gh, std_d_xy, std_d_gh, pool_std_xy, pool_std_gh, rmsd_xy, rmsd_gh = get_stats(x=x, y=y, x_err=x_err, y_err=y_err)
3126 d_row += [mean_d_xy, mean_d_gh, std_d_xy, std_d_gh, pool_std_xy, pool_std_gh, rmsd_xy, rmsd_gh]
3127
3128 if type_i in ['kex', 'k_AB', 'pA']:
3129 if i == 0:
3130 headers += ["%s_%s_ref"%(j, type_i), "%s_%s_ref_std"%(j, type_i), "%s_%s_ni"%(j, type_i), "%s_%s_ni_std"%(j, type_i)]
3131
3132 ref_val = data_dic[ref_pipe_name]['all_pars'][j][k][0]
3133 ni_val = data_dic[pipe_name]['all_pars'][j][k][0]
3134 ref_err = data_dic[ref_pipe_name]['all_pars_err'][j][k][0]
3135 ni_err = data_dic[pipe_name]['all_pars_err'][j][k][0]
3136
3137 d_row += [ref_val, ref_err, ni_val, ni_err]
3138
3139 k += 1
3140
3141
3142
3143 all_data.append(d_row)
3144
3145 file_name = ni_pipe_names[0] + '.txt'
3146 path = self.results_dir
3147 DAT = asarray(all_data)
3148 fmt = " ".join(["%6d"] + ["%8.2f"] + ["%21.6e"] * (DAT.shape[1]-2))
3149 headers_fmt = " ".join(["%4s"] + ["%8s"] + ["%21s"] * (DAT.shape[1]-2))
3150 headers_num = []
3151 for k, header_k in enumerate(headers):
3152 headers_num.append("%i_%s"%(k, header_k))
3153 headers_out = headers_fmt % tuple(headers_num)
3154 savetxt(path+sep+file_name, DAT, fmt=fmt, header=headers_out, delimiter=' ')
3155
3156
3162
3163
3165 """Store to self and settings dictionary"""
3166
3167
3168 self.settings[key] = value
3169
3170
3171 setattr(self, key, value)
3172
3173
3177
3178
3182
3183
3188
3189
3194