1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """The model-free analysis parameter functions."""
24
25
26 from math import pi
27 from numpy import array, float64, identity, int8, zeros
28 from re import match, search
29
30
31 from lib.errors import RelaxError
32 from pipe_control import diffusion_tensor
33 from pipe_control.mol_res_spin import spin_loop
34 from specific_analyses.model_free.model import determine_model_type
35
36
93
94
96 """Function for assembling a list of all the model parameter names.
97
98 @param model_type: The model-free model type. This must be one of 'mf', 'local_tm',
99 'diff', or 'all'.
100 @type model_type: str
101 @param spin_id: The spin identification string.
102 @type spin_id: str
103 @return: A list containing all the parameters of the model-free model.
104 @rtype: list of str
105 """
106
107
108 param_names = []
109
110
111 if model_type == 'diff' or model_type == 'all':
112
113 if cdp.diff_tensor.type == 'sphere':
114 param_names.append('tm')
115
116
117 elif cdp.diff_tensor.type == 'spheroid':
118 param_names.append('tm')
119 param_names.append('Da')
120 param_names.append('theta')
121 param_names.append('phi')
122
123
124 elif cdp.diff_tensor.type == 'ellipsoid':
125 param_names.append('tm')
126 param_names.append('Da')
127 param_names.append('Dr')
128 param_names.append('alpha')
129 param_names.append('beta')
130 param_names.append('gamma')
131
132
133 if model_type != 'diff':
134
135 for spin in spin_loop(spin_id):
136
137 if not spin.select:
138 continue
139
140
141 param_names = param_names + spin.params
142
143
144 return param_names
145
146
318
319
321 """Create and return the scaling matrix.
322
323 If the spin argument is supplied, then the spin_id argument will be ignored.
324
325 @param num_params: The number of parameters in the model.
326 @type num_params: int
327 @keyword model_type: The model type, one of 'all', 'diff', 'mf', or 'local_tm'.
328 @type model_type: str
329 @keyword spin: The spin data container.
330 @type spin: SpinContainer instance
331 @keyword spin_id: The spin identification string.
332 @type spin_id: str
333 @return: The diagonal and square scaling matrix.
334 @rtype: numpy diagonal matrix
335 """
336
337
338 if num_params == 0:
339 scaling_matrix = zeros((0, 0), float64)
340 else:
341 scaling_matrix = identity(num_params, float64)
342 i = 0
343
344
345 if not scaling:
346 return scaling_matrix
347
348
349 ti_scaling = 1e-12
350
351
352 if model_type == 'diff' or model_type == 'all':
353
354 if cdp.diff_tensor.type == 'sphere':
355
356 scaling_matrix[i, i] = ti_scaling
357
358
359 i = i + 1
360
361
362 elif cdp.diff_tensor.type == 'spheroid':
363
364 scaling_matrix[i, i] = ti_scaling
365 scaling_matrix[i+1, i+1] = 1e7
366 scaling_matrix[i+2, i+2] = 1.0
367 scaling_matrix[i+3, i+3] = 1.0
368
369
370 i = i + 4
371
372
373 elif cdp.diff_tensor.type == 'ellipsoid':
374
375 scaling_matrix[i, i] = ti_scaling
376 scaling_matrix[i+1, i+1] = 1e7
377 scaling_matrix[i+2, i+2] = 1.0
378 scaling_matrix[i+3, i+3] = 1.0
379 scaling_matrix[i+4, i+4] = 1.0
380 scaling_matrix[i+5, i+5] = 1.0
381
382
383 i = i + 6
384
385
386 if model_type != 'diff':
387
388 if spin:
389 loop = [spin]
390 else:
391 loop = spin_loop(spin_id)
392
393
394 for spin in loop:
395
396 if not spin.select:
397 continue
398
399
400 for k in range(len(spin.params)):
401
402 if spin.params[k] == 'local_tm' or search('^t', spin.params[k]):
403 scaling_matrix[i, i] = ti_scaling
404
405
406 elif spin.params[k] == 'rex':
407 scaling_matrix[i, i] = 1.0 / (2.0 * pi * cdp.spectrometer_frq[cdp.ri_ids[0]]) ** 2
408
409
410 elif spin.params[k] == 'r':
411 scaling_matrix[i, i] = 1e-10
412
413
414 elif spin.params[k] == 'csa':
415 scaling_matrix[i, i] = 1e-4
416
417
418 i = i + 1
419
420
421 return scaling_matrix
422
423
425 """Calculate and return the Rex conversion factor.
426
427 @return: The Rex conversion factor.
428 @rtype: float
429 """
430
431
432 if not hasattr(cdp, 'spectrometer_frq'):
433 raise RelaxError("No spectrometer frequency information is present in the current data pipe.")
434
435
436 if hasattr(cdp, 'ri_ids'):
437 frq = cdp.spectrometer_frq[cdp.ri_ids[0]]
438
439
440 else:
441 frqs = sorted(cdp.spectrometer_frq.values())
442 frq = frqs[-1]
443
444
445 return 1.0 / (2.0 * pi * frq)**2
446
447
449 """Disassemble the model-free parameter vector.
450
451 @param model_type: The model-free model type. This must be one of 'mf', 'local_tm',
452 'diff', or 'all'.
453 @type model_type: str
454 @keyword param_vector: The model-free parameter vector.
455 @type param_vector: numpy array
456 @keyword spin: The spin data container. If this argument is supplied, then the spin_id
457 argument will be ignored.
458 @type spin: SpinContainer instance
459 @keyword spin_id: The spin identification string.
460 @type spin_id: str
461 @keyword sim_index: The optional MC simulation index.
462 @type sim_index: int
463 """
464
465
466 param_index = 0
467
468
469 if sim_index != None and (model_type == 'diff' or model_type == 'all'):
470
471 if cdp.diff_tensor.type == 'sphere':
472
473 cdp.diff_tensor.set(param='tm', value=param_vector[0], category='sim', sim_index=sim_index)
474
475
476 param_index = param_index + 1
477
478
479 elif cdp.diff_tensor.type == 'spheroid':
480
481 cdp.diff_tensor.set(param='tm', value=param_vector[0], category='sim', sim_index=sim_index)
482 cdp.diff_tensor.set(param='Da', value=param_vector[1], category='sim', sim_index=sim_index)
483 cdp.diff_tensor.set(param='theta', value=param_vector[2], category='sim', sim_index=sim_index)
484 cdp.diff_tensor.set(param='phi', value=param_vector[3], category='sim', sim_index=sim_index)
485 diffusion_tensor.fold_angles(sim_index=sim_index)
486
487
488 param_index = param_index + 4
489
490
491 elif cdp.diff_tensor.type == 'ellipsoid':
492
493 cdp.diff_tensor.set(param='tm', value=param_vector[0], category='sim', sim_index=sim_index)
494 cdp.diff_tensor.set(param='Da', value=param_vector[1], category='sim', sim_index=sim_index)
495 cdp.diff_tensor.set(param='Dr', value=param_vector[2], category='sim', sim_index=sim_index)
496 cdp.diff_tensor.set(param='alpha', value=param_vector[3], category='sim', sim_index=sim_index)
497 cdp.diff_tensor.set(param='beta', value=param_vector[4], category='sim', sim_index=sim_index)
498 cdp.diff_tensor.set(param='gamma', value=param_vector[5], category='sim', sim_index=sim_index)
499 diffusion_tensor.fold_angles(sim_index=sim_index)
500
501
502 param_index = param_index + 6
503
504
505 elif model_type == 'diff' or model_type == 'all':
506
507 if cdp.diff_tensor.type == 'sphere':
508
509 cdp.diff_tensor.set(param='tm', value=param_vector[0])
510
511
512 param_index = param_index + 1
513
514
515 elif cdp.diff_tensor.type == 'spheroid':
516
517 cdp.diff_tensor.set(param='tm', value=param_vector[0])
518 cdp.diff_tensor.set(param='Da', value=param_vector[1])
519 cdp.diff_tensor.set(param='theta', value=param_vector[2])
520 cdp.diff_tensor.set(param='phi', value=param_vector[3])
521 diffusion_tensor.fold_angles()
522
523
524 param_index = param_index + 4
525
526
527 elif cdp.diff_tensor.type == 'ellipsoid':
528
529 cdp.diff_tensor.set(param='tm', value=param_vector[0])
530 cdp.diff_tensor.set(param='Da', value=param_vector[1])
531 cdp.diff_tensor.set(param='Dr', value=param_vector[2])
532 cdp.diff_tensor.set(param='alpha', value=param_vector[3])
533 cdp.diff_tensor.set(param='beta', value=param_vector[4])
534 cdp.diff_tensor.set(param='gamma', value=param_vector[5])
535 diffusion_tensor.fold_angles()
536
537
538 param_index = param_index + 6
539
540
541 if model_type != 'diff':
542
543 if spin:
544 loop = [spin]
545 else:
546 loop = spin_loop(spin_id)
547
548
549 for spin in loop:
550
551 if not spin.select:
552 continue
553
554
555 for j in range(len(spin.params)):
556
557 if spin.params[j] == 'local_tm':
558 if sim_index == None:
559 spin.local_tm = param_vector[param_index]
560 else:
561 spin.local_tm_sim[sim_index] = param_vector[param_index]
562
563
564 elif spin.params[j] == 's2':
565 if sim_index == None:
566 spin.s2 = param_vector[param_index]
567 else:
568 spin.s2_sim[sim_index] = param_vector[param_index]
569
570
571 elif spin.params[j] == 's2f':
572 if sim_index == None:
573 spin.s2f = param_vector[param_index]
574 else:
575 spin.s2f_sim[sim_index] = param_vector[param_index]
576
577
578 elif spin.params[j] == 's2s':
579 if sim_index == None:
580 spin.s2s = param_vector[param_index]
581 else:
582 spin.s2s_sim[sim_index] = param_vector[param_index]
583
584
585 elif spin.params[j] == 'te':
586 if sim_index == None:
587 spin.te = param_vector[param_index]
588 else:
589 spin.te_sim[sim_index] = param_vector[param_index]
590
591
592 elif spin.params[j] == 'tf':
593 if sim_index == None:
594 spin.tf = param_vector[param_index]
595 else:
596 spin.tf_sim[sim_index] = param_vector[param_index]
597
598
599 elif spin.params[j] == 'ts':
600 if sim_index == None:
601 spin.ts = param_vector[param_index]
602 else:
603 spin.ts_sim[sim_index] = param_vector[param_index]
604
605
606 elif spin.params[j] == 'rex':
607 if sim_index == None:
608 spin.rex = param_vector[param_index]
609 else:
610 spin.rex_sim[sim_index] = param_vector[param_index]
611
612
613 elif spin.params[j] == 'r':
614 if sim_index == None:
615 spin.r = param_vector[param_index]
616 else:
617 spin.r_sim[sim_index] = param_vector[param_index]
618
619
620 elif spin.params[j] == 'csa':
621 if sim_index == None:
622 spin.csa = param_vector[param_index]
623 else:
624 spin.csa_sim[sim_index] = param_vector[param_index]
625
626
627 else:
628 raise RelaxError("Unknown parameter.")
629
630
631 param_index = param_index + 1
632
633
634 if model_type != 'diff':
635
636 if spin:
637 loop = [spin]
638 else:
639 loop = spin_loop(spin_id)
640
641
642 for spin in loop:
643
644 if not spin.select:
645 continue
646
647
648 if sim_index == None:
649
650 if 's2' not in spin.params and 's2f' in spin.params and 's2s' in spin.params:
651 spin.s2 = spin.s2f * spin.s2s
652
653
654 if 's2f' not in spin.params and 's2' in spin.params and 's2s' in spin.params:
655 if spin.s2s == 0.0:
656 spin.s2f = 1e99
657 else:
658 spin.s2f = spin.s2 / spin.s2s
659
660
661 if 's2s' not in spin.params and 's2' in spin.params and 's2f' in spin.params:
662 if spin.s2f == 0.0:
663 spin.s2s = 1e99
664 else:
665 spin.s2s = spin.s2 / spin.s2f
666
667
668 else:
669
670 if 's2' not in spin.params and 's2f' in spin.params and 's2s' in spin.params:
671 spin.s2_sim[sim_index] = spin.s2f_sim[sim_index] * spin.s2s_sim[sim_index]
672
673
674 if 's2f' not in spin.params and 's2' in spin.params and 's2s' in spin.params:
675 if spin.s2s_sim[sim_index] == 0.0:
676 spin.s2f_sim[sim_index] = 1e99
677 else:
678 spin.s2f_sim[sim_index] = spin.s2_sim[sim_index] / spin.s2s_sim[sim_index]
679
680
681 if 's2s' not in spin.params and 's2' in spin.params and 's2f' in spin.params:
682 if spin.s2f_sim[sim_index] == 0.0:
683 spin.s2s_sim[sim_index] = 1e99
684 else:
685 spin.s2s_sim[sim_index] = spin.s2_sim[sim_index] / spin.s2f_sim[sim_index]
686
687
688 -def linear_constraints(num_params, model_type=None, spin=None, spin_id=None, scaling_matrix=None):
689 """Set up the model-free linear constraint matrices A and b.
690
691 Standard notation
692 =================
693
694 The order parameter constraints are::
695
696 0 <= S2 <= 1
697 0 <= S2f <= 1
698 0 <= S2s <= 1
699
700 By substituting the formula S2 = S2f.S2s into the above inequalities, the additional two
701 inequalities can be derived::
702
703 S2 <= S2f
704 S2 <= S2s
705
706 Correlation time constraints are::
707
708 te >= 0
709 tf >= 0
710 ts >= 0
711
712 tf <= ts
713
714 te, tf, ts <= 2 * tm
715
716 Additional constraints used include::
717
718 Rex >= 0
719 0.9e-10 <= r <= 2e-10
720 -300e-6 <= CSA <= 0
721
722
723 Rearranged notation
724 ===================
725
726 The above inequality constraints can be rearranged into::
727
728 S2 >= 0
729 -S2 >= -1
730 S2f >= 0
731 -S2f >= -1
732 S2s >= 0
733 -S2s >= -1
734 S2f - S2 >= 0
735 S2s - S2 >= 0
736 te >= 0
737 tf >= 0
738 ts >= 0
739 ts - tf >= 0
740 Rex >= 0
741 r >= 0.9e-10
742 -r >= -2e-10
743 CSA >= -300e-6
744 -CSA >= 0
745
746
747 Matrix notation
748 ===============
749
750 In the notation A.x >= b, where A is an matrix of coefficients, x is an array of parameter
751 values, and b is a vector of scalars, these inequality constraints are::
752
753 | 1 0 0 0 0 0 0 0 0 | | 0 |
754 | | | |
755 |-1 0 0 0 0 0 0 0 0 | | -1 |
756 | | | |
757 | 0 1 0 0 0 0 0 0 0 | | 0 |
758 | | | |
759 | 0 -1 0 0 0 0 0 0 0 | | -1 |
760 | | | |
761 | 0 0 1 0 0 0 0 0 0 | | S2 | | 0 |
762 | | | | | |
763 | 0 0 -1 0 0 0 0 0 0 | | S2f | | -1 |
764 | | | | | |
765 |-1 1 0 0 0 0 0 0 0 | | S2s | | 0 |
766 | | | | | |
767 |-1 0 1 0 0 0 0 0 0 | | te | | 0 |
768 | | | | | |
769 | 0 0 0 1 0 0 0 0 0 | . | tf | >= | 0 |
770 | | | | | |
771 | 0 0 0 0 1 0 0 0 0 | | ts | | 0 |
772 | | | | | |
773 | 0 0 0 0 0 1 0 0 0 | | Rex | | 0 |
774 | | | | | |
775 | 0 0 0 0 -1 1 0 0 0 | | r | | 0 |
776 | | | | | |
777 | 0 0 0 0 0 0 1 0 0 | | CSA | | 0 |
778 | | | |
779 | 0 0 0 0 0 0 0 1 0 | | 0.9e-10 |
780 | | | |
781 | 0 0 0 0 0 0 0 -1 0 | | -2e-10 |
782 | | | |
783 | 0 0 0 0 0 0 0 0 1 | | -300e-6 |
784 | | | |
785 | 0 0 0 0 0 0 0 0 -1 | | 0 |
786
787
788 @param num_params: The number of parameters in the model.
789 @type num_params: int
790 @keyword model_type: The model type, one of 'all', 'diff', 'mf', or 'local_tm'.
791 @type model_type: str
792 @keyword spin: The spin data container. If this argument is supplied, then the
793 spin_id argument will be ignored.
794 @type spin: SpinContainer instance
795 @keyword spin_id: The spin identification string.
796 @type spin_id: str
797 @keyword scaling_matrix: The diagonal, square scaling matrix.
798 @type scaling_matrix: numpy diagonal matrix
799 """
800
801
802 upper_time_limit = 1
803
804
805 A = []
806 b = []
807 zero_array = zeros(num_params, float64)
808 i = 0
809 j = 0
810
811
812 if model_type != 'mf' and diffusion_tensor.diff_data_exists():
813
814 if cdp.diff_tensor.type == 'sphere':
815
816 A.append(zero_array * 0.0)
817 A.append(zero_array * 0.0)
818 A[j][i] = 1.0
819 A[j+1][i] = -1.0
820 b.append(0.0 / scaling_matrix[i, i])
821 b.append(-200.0 * 1e-9 / scaling_matrix[i, i])
822 i = i + 1
823 j = j + 2
824
825
826 elif cdp.diff_tensor.type == 'spheroid':
827
828 A.append(zero_array * 0.0)
829 A.append(zero_array * 0.0)
830 A[j][i] = 1.0
831 A[j+1][i] = -1.0
832 b.append(0.0 / scaling_matrix[i, i])
833 b.append(-200.0 * 1e-9 / scaling_matrix[i, i])
834 i = i + 1
835 j = j + 2
836
837
838 if cdp.diff_tensor.spheroid_type == 'prolate':
839 A.append(zero_array * 0.0)
840 A[j][i] = 1.0
841 b.append(0.0 / scaling_matrix[i, i])
842 i = i + 1
843 j = j + 1
844
845
846 i = i + 2
847
848
849 elif cdp.diff_tensor.spheroid_type == 'oblate':
850 A.append(zero_array * 0.0)
851 A[j][i] = -1.0
852 b.append(0.0 / scaling_matrix[i, i])
853 i = i + 1
854 j = j + 1
855
856
857 i = i + 2
858
859 else:
860
861 i = i + 3
862
863
864 elif cdp.diff_tensor.type == 'ellipsoid':
865
866 A.append(zero_array * 0.0)
867 A.append(zero_array * 0.0)
868 A[j][i] = 1.0
869 A[j+1][i] = -1.0
870 b.append(0.0 / scaling_matrix[i, i])
871 b.append(-200.0 * 1e-9 / scaling_matrix[i, i])
872 i = i + 1
873 j = j + 2
874
875
876 A.append(zero_array * 0.0)
877 A[j][i] = 1.0
878 b.append(0.0 / scaling_matrix[i, i])
879 i = i + 1
880 j = j + 1
881
882
883 A.append(zero_array * 0.0)
884 A.append(zero_array * 0.0)
885 A[j][i] = 1.0
886 A[j+1][i] = -1.0
887 b.append(0.0 / scaling_matrix[i, i])
888 b.append(-1.0 / scaling_matrix[i, i])
889 i = i + 1
890 j = j + 2
891
892
893 i = i + 3
894
895
896 if model_type != 'diff':
897
898 if spin:
899 loop = [spin]
900 else:
901 loop = spin_loop(spin_id)
902
903
904 for spin in loop:
905
906 if not spin.select:
907 continue
908
909
910 old_i = i
911
912
913 for l in range(len(spin.params)):
914
915 if spin.params[l] == 'local_tm':
916 if upper_time_limit:
917
918 A.append(zero_array * 0.0)
919 A.append(zero_array * 0.0)
920 A[j][i] = 1.0
921 A[j+1][i] = -1.0
922 b.append(0.0 / scaling_matrix[i, i])
923 b.append(-200.0 * 1e-9 / scaling_matrix[i, i])
924 j = j + 2
925 else:
926
927 A.append(zero_array * 0.0)
928 A[j][i] = 1.0
929 b.append(0.0 / scaling_matrix[i, i])
930 j = j + 1
931
932
933 elif match('s2', spin.params[l]):
934
935 A.append(zero_array * 0.0)
936 A.append(zero_array * 0.0)
937 A[j][i] = 1.0
938 A[j+1][i] = -1.0
939 b.append(0.0 / scaling_matrix[i, i])
940 b.append(-1.0 / scaling_matrix[i, i])
941 j = j + 2
942
943
944 if spin.params[l] == 's2':
945 for m in range(len(spin.params)):
946 if spin.params[m] == 's2f' or spin.params[m] == 's2s':
947 A.append(zero_array * 0.0)
948 A[j][i] = -1.0
949 A[j][old_i+m] = 1.0
950 b.append(0.0)
951 j = j + 1
952
953
954 elif match('t[efs]', spin.params[l]):
955
956 A.append(zero_array * 0.0)
957 A[j][i] = 1.0
958 b.append(0.0 / scaling_matrix[i, i])
959 j = j + 1
960
961
962 if spin.params[l] == 'ts':
963 for m in range(len(spin.params)):
964 if spin.params[m] == 'tf':
965 A.append(zero_array * 0.0)
966 A[j][i] = 1.0
967 A[j][old_i+m] = -1.0
968 b.append(0.0)
969 j = j + 1
970
971
972 if upper_time_limit:
973 if not spin.params[l] == 'tf':
974 if model_type == 'mf':
975 A.append(zero_array * 0.0)
976 A[j][i] = -1.0
977 b.append(-2.0 * cdp.diff_tensor.tm / scaling_matrix[i, i])
978 else:
979 A.append(zero_array * 0.0)
980 A[j][0] = 2.0
981 A[j][i] = -1.0
982 b.append(0.0)
983
984 j = j + 1
985
986
987 elif spin.params[l] == 'rex':
988 A.append(zero_array * 0.0)
989 A[j][i] = 1.0
990 b.append(0.0 / scaling_matrix[i, i])
991 j = j + 1
992
993
994 elif spin.params[l] == 'r':
995
996 A.append(zero_array * 0.0)
997 A.append(zero_array * 0.0)
998 A[j][i] = 1.0
999 A[j+1][i] = -1.0
1000 b.append(0.9e-10 / scaling_matrix[i, i])
1001 b.append(-2e-10 / scaling_matrix[i, i])
1002 j = j + 2
1003
1004
1005 elif spin.params[l] == 'csa':
1006
1007 A.append(zero_array * 0.0)
1008 A.append(zero_array * 0.0)
1009 A[j][i] = 1.0
1010 A[j+1][i] = -1.0
1011 b.append(-300e-6 / scaling_matrix[i, i])
1012 b.append(0.0 / scaling_matrix[i, i])
1013 j = j + 2
1014
1015
1016 i = i + 1
1017
1018
1019 A = array(A, int8)
1020 b = array(b, float64)
1021
1022 return A, b
1023
1024
1026 """Return the units for the Rex parameter.
1027
1028 @return: The field strength dependent Rex units.
1029 @rtype: str
1030 """
1031
1032
1033 if not hasattr(cdp, 'frq_labels') or len(cdp.frq_labels) == 0:
1034 return ''
1035
1036
1037 return cdp.frq_labels[0] + ' MHz'
1038