PYCON 89
Absences.py Guest on 14th May 2020 10:50:14 AM
  1. from cctbx import sgtbx, uctbx
  2. from cctbx import miller
  3. from cctbx.array_family import flex
  4. from cctbx import sgtbx
  5. from cctbx import crystal
  6. from libtbx import table_utils
  7. import math,sys
  8.  
  9.  
  10.  
  11. #                           name     hkl selector  condition
  12. absence_and_conditions = { "2_1 (a)" : [(None,0,0),  (1.0/2.0,0,0)],
  13.                            "2_1 (b)" : [(0,None,0),  (0,1.0/2.0,0)],
  14.                            "2_1 (c)" : [(0,0,None),  (0,0,1.0/2.0)],
  15.  
  16.                            "3_1 (c)" : [(0,0,None),  (0,0,1.0/3.0)],
  17.                            "3_2 (c)" : [(0,0,None),  (0,0,1.0/3.0)],
  18.  
  19.                            "4_1 (c)" : [(0,0,None),  (0,0,1.0/4.0)],
  20.                            "4_2 (c)" : [(0,0,None),  (0,0,2.0/4.0)],
  21.                            "4_3 (c)" : [(0,0,None),  (0,0,3.0/4.0)],
  22.  
  23.                            "4_1 (a)" : [(None,0,0),  (1.0/4.0,0,0)],
  24.                            "4_2 (a)" : [(None,0,0),  (2.0/4.0,0,0)],
  25.                            "4_3 (a)" : [(None,0,0),  (3.0/4.0,0,0)],
  26.  
  27.                            "4_1 (b)" : [(0,None,0),  (0,1.0/4.0,0)],
  28.                            "4_2 (b)" : [(0,None,0),  (0,2.0/4.0,0)],
  29.                            "4_3 (b)" : [(0,None,0),  (0,3.0/4.0,0)],
  30.  
  31.                            "6_1 (c)" : [(0,0,None),  (0,0,1.0/6.0)],
  32.                            "6_2 (c)" : [(0,0,None),  (0,0,2.0/6.0)],
  33.                            "6_3 (c)" : [(0,0,None),  (0,0,3.0/6.0)],
  34.                            "6_4 (c)" : [(0,0,None),  (0,0,4.0/6.0)],
  35.                            "6_5 (c)" : [(0,0,None),  (0,0,5.0/6.0)],
  36.  
  37.                            "b (a)"   : [(0,None,None),  (0,1.0/2.0,0)],
  38.                            "c (a)"   : [(0,None,None),  (0,0,1.0/2.0)],
  39.                            "n (a)"   : [(0,None,None),  (0,1.0/2.0,1.0/2.0)],
  40.                            "d (a)"   : [(0,None,None),  (0,1.0/4.0,1.0/4.0)],
  41.  
  42.                            "a (b)"   : [(None,0,None),  (1.0/2.0,0,0)],
  43.                            "c (b)"   : [(None,0,None),  (0,0,1.0/2.0)],
  44.                            "n (b)"   : [(None,0,None),  (1.0/2.0,0,1.0/2.0)],
  45.                            "d (b)"   : [(None,0,None),  (1.0/4.0,0,1.0/4.0)],
  46.  
  47.                            "a (c)"   : [(None,None,0),  (1.0/2.0,0,0)],
  48.                            "b (c)"   : [(None,None,0),  (0,1.0/2.0,0)],
  49.                            "n (c)"   : [(None,None,0),  (1.0/2.0,1.0/2.0,0)],
  50.                            "d (c)"   : [(None,None,0),  (1.0/4.0,1.0/4.0,0)]
  51.                          }
  52.  
  53.  
  54. absence_classes = { "along a 2" : ["2_1 (a)"],
  55.                     "along a"   : ["b (a)", "c (a)", "n (a)", "d (a)"],
  56.                     "along b 4" : ["4_1 (b)", "4_2 (b)", "4_3 (b)"],
  57.                     "along b 2" : ["2_1 (b)"],
  58.                     "along b"   : ["a (b)", "c (b)", "n (b)", "d (b)"],
  59.                     "along c"   : ["a (c)", "b (c)", "n (c)", "d (c)"],
  60.                     "along c 2" : ["2_1 (c)"],
  61.                     "along c 3" : ["3_1 (c)", "3_2 (c)"],
  62.                     "along c 4" : ["4_1 (c)", "4_2 (c)", "4_3 (c)"],
  63.                     "along c 6" : ["6_1 (c)", "6_2 (c)", "6_3 (c)", "6_4 (c)", "6_5 (c)"] }
  64.  
  65.  
  66. along_a   = absence_classes["along a"]
  67. along_a_2 = absence_classes["along a 2"]
  68.  
  69. along_b_4 = absence_classes["along b 4"]
  70. along_b   = absence_classes["along b"]
  71. along_b_2 = absence_classes["along b 2"]
  72.  
  73. along_c   = absence_classes["along c"]
  74. along_c_2 = absence_classes["along c 2"]
  75. along_c_3 = absence_classes["along c 3"]
  76. along_c_4 = absence_classes["along c 4"]
  77. along_c_6 = absence_classes["along c 6"]
  78.  
  79. absences_via_intensity_symmetry = {
  80. "P -1"       : []                      ,# l>0 or (l==0 and (h>0 or (h==0 and k>=0)))
  81.  
  82. "P 1 2/m 1"  : along_b+along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  83. "C 1 2/m 1"  : along_b+along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  84. "A 1 2/m 1"  : along_b+along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  85. "I 1 2/m 1"  : along_b+along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  86.  
  87. "P 1 1 2/m"  : along_c+along_c_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  88. "A 1 1 2/m"  : along_c+along_c_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  89. "B 1 1 2/m"  : along_c+along_c_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  90. "I 1 1 2/m"  : along_c+along_c_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  91.  
  92. "P 2/m 1 1"  : along_b+along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  93. "B 2/m 1 1"  : along_b+along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  94. "C 2/m 1 1"  : along_b+along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  95. "I 2/m 1 1"  : along_b+along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  96.  
  97. "P m m m"    : along_a+along_b+along_c+along_c_2 ,# h>=0 and k>=0 and l>=0
  98. "C m m m"    : along_a+along_b+along_c+along_c_2 ,# h>=0 and k>=0 and l>=0
  99. "A m m m"    : along_a+along_b+along_c+along_c_2 ,# h>=0 and k>=0 and l>=0
  100. "B m m m"    : along_a+along_b+along_c+along_c_2 ,# h>=0 and k>=0 and l>=0
  101. "F m m m"    : along_a+along_b+along_c+along_c_2 ,# h>=0 and k>=0 and l>=0
  102. "I m m m"    : along_a+along_b+along_c+along_c_2 ,# h>=0 and k>=0 and l>=0
  103.  
  104. "P 4/m"      : along_c_4+along_c ,# l>=0 and ((h>=0 and k>0) or (h==0 and k==0))
  105. "I 4/m"      : along_c_4+along_c         ,# l>=0 and ((h>=0 and k>0) or (h==0 and k==0))
  106.  
  107. "P 4/m m m"  : along_a+along_c_4+along_c         ,# h>=k and k>=0 and l>=0
  108. "I 4/m m m"  : along_a+along_c_4+along_c         ,# h>=k and k>=0 and l>=0
  109.  
  110. "P -3"       : along_b+along_c_3+along_c         ,# (h>=0 and k>0) or (h==0 and k==0 and l>=0)
  111. "R -3 :H"    : along_b+along_c_3+along_c        ,# (h>=0 and k>0) or (h==0 and k==0 and l>=0)
  112. "R -3 :R"    : along_b+along_c_3+along_c         ,# (h>=0 and k>0) or (h==0 and k==0 and l>=0)
  113.  
  114. "P -3 1 m"   : along_b+along_c_3+along_c+along_b_2        ,# h>=k and k>=0 and (k>0 or l>=0)
  115. "P -3 m 1"   : along_b+along_c_3+along_c+along_b_2        ,# h>=k and k>=0 and (h>k or l>=0)
  116. "R -3 m :H"  : along_b+along_c_3+along_c+along_b_2        ,# h>=k and k>=0 and (h>k or l>=0)
  117. "R -3 m :R"  : along_b+along_c_3+along_c+along_b_2        ,# h>=k and k>=0 and (h>k or l>=0)
  118.  
  119. "P 6/m"      : along_c+along_b+along_c_6         ,# l>=0 and ((h>=0 and k>0) or (h==0 and k==0))
  120. "P 6/m m m"  : along_c+along_b+along_c_6         ,# h>=k and k>=0 and l>=0
  121.  
  122. "P m -3"     : along_c+along_c_2+along_c_4+along_b_2                ,# h>=0 and ((l>=h and k>h) or (l==h and k==h))
  123. "F m -3"     : along_c+along_c_2+along_c_4+along_b_2                 ,# h>=0 and ((l>=h and k>h) or (l==h and k==h))
  124. "I m -3"     : along_c+along_c_2+along_c_4+along_b_2                 ,# h>=0 and ((l>=h and k>h) or (l==h and k==h))
  125.  
  126. "P m -3 m"   : along_b+along_b_4+along_b_2                 ,# k>=l and l>=h and h>=0
  127. "F m -3 m"   : along_b+along_b_4+along_b_2                 ,# k>=l and l>=h and h>=0
  128. "I m -3 m"   : along_b+along_b_4+along_b_2                 ,# k>=l and l>=h and h>=0
  129. }
  130.  
  131.  
  132. absences_via_intensity_symmetry_prot = {
  133. "P -1"       : []                      ,# l>0 or (l==0 and (h>0 or (h==0 and k>=0)))
  134.  
  135. "P 1 2/m 1"  : along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  136. "C 1 2/m 1"  : along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  137. "A 1 2/m 1"  : along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  138. "I 1 2/m 1"  : along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  139.  
  140. "P 1 1 2/m"  : along_c_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  141. "A 1 1 2/m"  : along_c_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  142. "B 1 1 2/m"  : along_c_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  143. "I 1 1 2/m"  : along_c_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  144.  
  145. "P 2/m 1 1"  : along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  146. "B 2/m 1 1"  : along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  147. "C 2/m 1 1"  : along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  148. "I 2/m 1 1"  : along_b_2       ,# k>=0 and (l>0 or (l==0 and h>=0))
  149.  
  150. "P m m m"    : along_a_2+along_b_2+along_c_2 ,# h>=0 and k>=0 and l>=0
  151. "C m m m"    : along_c_2 ,# h>=0 and k>=0 and l>=0
  152. "A m m m"    : along_a_2 ,# h>=0 and k>=0 and l>=0
  153. "B m m m"    : along_b_2 ,# h>=0 and k>=0 and l>=0
  154. "F m m m"    : along_a_2+along_b_2+along_c_2 ,# h>=0 and k>=0 and l>=0
  155. "I m m m"    : along_a_2+along_b_2+along_c_2 ,# h>=0 and k>=0 and l>=0
  156.  
  157. "P 4/m"      : along_c_4 ,# l>=0 and ((h>=0 and k>0) or (h==0 and k==0))
  158. "I 4/m"      : along_c_4         ,# l>=0 and ((h>=0 and k>0) or (h==0 and k==0))
  159.  
  160. "P 4/m m m"  : along_a_2+along_c_4         ,# h>=k and k>=0 and l>=0
  161. "I 4/m m m"  : along_a_2+along_c_4         ,# h>=k and k>=0 and l>=0
  162.  
  163. "P -3"       : along_c_3        ,# (h>=0 and k>0) or (h==0 and k==0 and l>=0)
  164. "R -3 :H"    : along_c_3        ,# (h>=0 and k>0) or (h==0 and k==0 and l>=0)
  165. "R -3 :R"    : along_c_3         ,# (h>=0 and k>0) or (h==0 and k==0 and l>=0)
  166.  
  167. "P -3 1 m"   : along_c_3        ,# h>=k and k>=0 and (k>0 or l>=0)
  168. "P -3 m 1"   : along_c_3        ,# h>=k and k>=0 and (h>k or l>=0)
  169. "R -3 m :H"  : along_c_3        ,# h>=k and k>=0 and (h>k or l>=0)
  170. "R -3 m :R"  : along_c_3        ,# h>=k and k>=0 and (h>k or l>=0)
  171.  
  172. "P 6/m"      : along_c_6         ,# l>=0 and ((h>=0 and k>0) or (h==0 and k==0))
  173. "P 6/m m m"  : along_c_6         ,# h>=k and k>=0 and l>=0
  174.  
  175. "P m -3"     : along_c_2+along_c_4+along_b_2                ,# h>=0 and ((l>=h and k>h) or (l==h and k==h))
  176. "F m -3"     : along_c_2+along_c_4+along_b_2                 ,# h>=0 and ((l>=h and k>h) or (l==h and k==h))
  177. "I m -3"     : along_c_2+along_c_4+along_b_2                 ,# h>=0 and ((l>=h and k>h) or (l==h and k==h))
  178.  
  179. "P m -3 m"   : along_b_4+along_b_2                 ,# k>=l and l>=h and h>=0
  180. "F m -3 m"   : along_b_4+along_b_2                 ,# k>=l and l>=h and h>=0
  181. "I m -3 m"   : along_b_4+along_b_2                 ,# k>=l and l>=h and h>=0
  182. }
  183.  
  184.  
  185. lattice_from_string = { "1 (0, 0, 0) 0,1/2,1/2"   : "A",
  186.                         "1 (0, 0, 0) 1/2,0,1/2"   : "B",
  187.                         "1 (0, 0, 0) 1/2,1/2,0"   : "C",
  188.                         "1 (0, 0, 0) 1/2,1/2,1/2" : "I",
  189.                       }
  190.  
  191.  
  192. symbol_from_string = { "-2 (1, 0, 0) 0,1/2,0" : "b (a)",
  193.                        "-2 (1, 0, 0) 0,0,1/2" : "c (a)",
  194.                        "-2 (1, 0, 0) 0,1/2,1/2" : "n (a)",
  195.                        "2 (1, 0, 0) 1/2,0,0" : "2_1 (a)",
  196.                        "-2 (0, 0, 1) 1/2,1/2,0" : "n (c)",
  197.                        "6 (0, 0, 1) 0,0,2/3" : "6_4 (c)",
  198.                        "-2 (1, 0, 0) 0,1/4,1/4" : "d (a)",
  199.                        "6 (0, 0, 1) 0,0,5/6" : "6_5 (c)",
  200.                        "6 (0, 0, 1) 0,0,1/2" : "6_3 (c)",
  201.                        "-2 (0, 1, 0) 1/2,0,0" : "a (b)",
  202.                        "2 (0, 0, 1) 0,0,1/2" : "2_1 (c)",
  203.                        "4 (0, 0, 1) 0,0,1/4" : "4_1 (c)",
  204.                        "-2 (0, 0, 1) 1/4,1/4,0" : "d (c)",
  205.                        "4 (1, 0, 0) 1/4,0,0" : "4_1 (a)",
  206.                        "3 (0, 0, 1) 0,0,1/3" : "3_1 (c)",
  207.                        "4 (0, 0, 1) 0,0,3/4" : "4_3 (c)",
  208.                        "-2 (0, 0, 1) 0,1/2,0" : "b (c)",
  209.                        "6 (0, 0, 1) 0,0,1/3" : "6_2 (c)",
  210.                        "-2 (0, 1, 0) 1/2,0,1/2" : "n (b)",
  211.                        "4 (1, 0, 0) 3/4,0,0" : "4_3 (a)",
  212.                        "-2 (0, 1, 0) 0,0,1/2" : "c (b)",
  213.                        "6 (0, 0, 1) 0,0,1/6" : "6_1 (c)",
  214.                        "-2 (0, 0, 1) 1/2,0,0" : "a (c)",
  215.                        "2 (0, 1, 0) 0,1/2,0" : "2_1 (b)",
  216.                        "-2 (0, 1, 0) 1/4,0,1/4" : "d (b)",
  217.                        "4 (0, 0, 1) 0,0,1/2" : "4_2 (c)",
  218.                        "4 (1, 0, 0) 2/4,0,0" : "4_2 (a)",
  219.                        "3 (0, 0, 1) 0,0,2/3" : "3_2 (c)" }
  220.  
  221. string_from_symbol ={"c (b)"  :  "-2 (0, 1, 0) 0,0,1/2",
  222.                      "n (b)"  :  "-2 (0, 1, 0) 1/2,0,1/2",
  223.                      "a (b)"  :  "-2 (0, 1, 0) 1/2,0,0",
  224.                      "a (c)"  :  "-2 (0, 0, 1) 1/2,0,0",
  225.                      "n (c)"  :  "-2 (0, 0, 1) 1/2,1/2,0",
  226.                      "b (c)"  :  "-2 (0, 0, 1) 0,1/2,0",
  227.                      "b (a)"  :  "-2 (1, 0, 0) 0,1/2,0",
  228.                      "n (a)"  :  "-2 (1, 0, 0) 0,1/2,1/2",
  229.                      "c (a)"  :  "-2 (1, 0, 0) 0,0,1/2",
  230.                      "2_1 (b)":   "2 (0, 1, 0) 0,1/2,0",
  231.                      "2_1 (c)":   "2 (0, 0, 1) 0,0,1/2",
  232.                      "2_1 (a)":   "2 (1, 0, 0) 1/2,0,0",
  233.                      "d (b)"  :  "-2 (0, 1, 0) 1/4,0,1/4",
  234.                      "d (c)"  :  "-2 (0, 0, 1) 1/4,1/4,0",
  235.                      "d (a)"  :  "-2 (1, 0, 0) 0,1/4,1/4",
  236.                      "4_1 (c)":   "4 (0, 0, 1) 0,0,1/4",
  237.                      "4_2 (c)":   "4 (0, 0, 1) 0,0,1/2",
  238.                      "4_3 (c)":   "4 (0, 0, 1) 0,0,3/4",
  239.                      "4_1 (a)":   "4 (1, 0, 0) 1/4,0,0",
  240.                      "4_2 (a)":   "4 (1, 0, 0) 2/4,0,0",
  241.                      "4_3 (a)":   "4 (1, 0, 0) 3/4,0,0",
  242.                      "3_1 (c)":   "3 (0, 0, 1) 0,0,1/3",
  243.                      "3_2 (c)":   "3 (0, 0, 1) 0,0,2/3",
  244.                      "6_1 (c)":   "6 (0, 0, 1) 0,0,1/6",
  245.                      "6_5 (c)":   "6 (0, 0, 1) 0,0,5/6",
  246.                      "6_4 (c)":   "6 (0, 0, 1) 0,0,2/3",
  247.                      "6_2 (c)":   "6 (0, 0, 1) 0,0,1/3",
  248.                      "6_3 (c)":   "6 (0, 0, 1) 0,0,1/2" }
  249.  
  250.  
  251. class conditions_for_operator(object):
  252.   def __init__(self, s ):
  253.     r_info = sgtbx.rot_mx_info( s.r() )
  254.     t_info = sgtbx.translation_part_info( s )
  255.  
  256.     self.type = r_info.type()
  257.     self.ev   = r_info.ev()
  258.     self.trans= t_info.intrinsic_part()
  259.  
  260.   def condition(self):
  261.     cond = []
  262.     # get conditions for hkl
  263.     if self.type < 0:
  264.       if ( list(self.trans.as_double()) ).count(0) <3 :
  265.         for ii in ev:
  266.           if ii == 0:
  267.             cond.append( None )
  268.           else:
  269.             cond.append( 0 )
  270.       else:
  271.         cond = (None,None,None)
  272.     if type > 0:
  273.       if ( list(self.trans.as_double()) ).count(0) <3 :
  274.         for ii in ev:
  275.           if ii != 0:
  276.             cond.append(None)
  277.           else:
  278.             cond.append( 0 )
  279.  
  280.     if len(cond) == 0:
  281.       cond = [ None,None,None ]
  282.  
  283.     return [ cond, list(t) ]
  284.  
  285.   def absence_type(self):
  286.     #get conditions for translations
  287.     t = self.trans.as_double()
  288.     id_string = str(self.type)+" "+str(self.ev)+" "+str(self.trans)
  289.     op_name = "Non absent"
  290.     if symbol_from_string.has_key( id_string ):
  291.       op_name = symbol_from_string[ id_string ]
  292.     return op_name
  293.  
  294. class absences(object):
  295.   def __init__(self, mult=2.0, threshold=0.95):
  296.     self.lib = absence_and_conditions
  297.     self.absence_classes = absences_via_intensity_symmetry_prot
  298.     self.mult = mult
  299.     self.threshold = threshold
  300.  
  301.   def check(self, abs_type, hkl, return_bool=False ):
  302.     # check if it is there
  303.     message = "The reflection with index %s is "%(str(hkl))
  304.     if self.lib.has_key( abs_type ):
  305.       mask, condition = self.lib[ abs_type ]
  306.       mc = self.check_mask( hkl, mask )
  307.       cc = self.check_condition(hkl, condition )
  308.       if mc:
  309.         if cc:
  310.           message += " PRESENT under the %s operator"%(abs_type)
  311.         else:
  312.           message += " ABSENT under the %s operartor"%(abs_type)
  313.       else:
  314.         message += " not of the type that is absent under %s"%(abs_type)
  315.     else:
  316.       message =  "No such type of absence"
  317.     if return_bool:
  318.       return mc,cc
  319.     else:
  320.       return message
  321.  
  322.   def check_mask(self, hkl, mask):
  323.     mask_none_count =0
  324.     for ii in mask:
  325.       if ii is None:
  326.         mask_none_count += 1
  327.     count = 0
  328.     result = False
  329.     for m,h in zip(mask,hkl):
  330.       if m == h:
  331.         count += 1
  332.     if count == 3-mask_none_count:
  333.       result = True
  334.     return result
  335.  
  336.   def check_condition(self, hkl, condition):
  337.     """magic"""
  338.     result = False
  339.     x = hkl[0]*condition[0] + hkl[1]*condition[1] + hkl[2]*condition[2]
  340.     y = math.cos( x*2.0*math.pi )
  341.     z = 0.5*(math.tanh(self.mult*y)+1)
  342.     if z > self.threshold:
  343.       result = True
  344.     return result
  345.  
  346. class analyze_absences(object):
  347.   def __init__(self, miller_array, isigi_cut=3, out=None):
  348.     if out is None:
  349.       out = sys.stdout
  350.     self.out = out
  351.     self.cut = isigi_cut
  352.     self.miller_array = miller_array.deep_copy()
  353.  
  354.     self.n_abs        = []
  355.     self.n_n_abs      = []
  356.     self.n_tot        = []
  357.     self.n_abs_viol   = []
  358.     self.n_n_abs_viol = []
  359.     self.n_tot_viol   = []
  360.  
  361.     self.isi_abs      = []
  362.     self.isi_n_abs    = []
  363.     self.isi_tot      = []
  364.  
  365.     self.i_abs        = []
  366.     self.i_n_abs      = []
  367.     self.i_tot        = []
  368.  
  369.     self.op_name      = []
  370.     self.score        = []
  371.     self.present      = []
  372.  
  373.     self.table_text = """
  374.  
  375.  
  376. Systematic absences
  377. -------------------
  378.  
  379. The following table gives informaton about systematic absences.
  380. For each operator, the reflections are split in three classes:
  381.   Absent    : Reflections that are absent for this operator.
  382.   Non Absent: Reflection of the same type (i.e. (0,0,l)) as above, but they should be present.
  383.   Complement: All other reflections.
  384. For each class, the <I/sigI> is reported, as well as the number of 'violations'. A 'violation'
  385. is designated as a reflection for which a I/sigI criterion is not met. The criteria are
  386.  
  387.   Absent violation     : I/sigI > %2.1f
  388.   Non Absent violation : I/sigI < %2.1f
  389.   Complement violation : I/sigI < %2.1f
  390.  
  391. Operators with low associated violations for *both* absent and non absent reflections, are likely to
  392. be true screw axis or glide planes. Both the number of violations and their percentages are given.
  393. The number of violations within the 'complement' class, can be used as a comparison for the number
  394. of violations in the non-absent class.
  395.  
  396. """%(self.cut, self.cut, self.cut)
  397.  
  398.     assert self.miller_array.sigmas() is not None
  399.     #we need to have this in the standard setting please
  400.     self.sg = sgtbx.space_group_info( group = self.miller_array.space_group() )
  401.     self.cb_op = self.sg.change_of_basis_op_to_reference_setting()
  402.     self.miller_array = self.miller_array.change_basis( self.cb_op )
  403.     self.abs_check = absences()
  404.     self.check_conditions()
  405.  
  406.   def show(self):
  407.     print >> self.out, self.table_text
  408.     print >> self.out, self.table
  409.     print >> self.out
  410.     print >> self.out
  411.  
  412.  
  413.  
  414.   def score_isigi(self,isig, absent=False, a=30.0):
  415.     tmp = 0.5*(1+math.tanh( a*(isig-self.cut) ) )
  416.     if not absent:
  417.       return abs(math.log(tmp+1e-8))
  418.     else:
  419.       return abs(math.log(1-tmp+1e-8) )
  420.  
  421.  
  422.   def propose(self, ops, thres=1):
  423.     in_sg_and_seemingly_correct       = []
  424.     not_in_sg_and_seemingly_incorrect = []
  425.     in_sg_and_seemingly_incorrect     = []
  426.     observed_but_not_in_sg            = []
  427.     in_sg_but_no_observations         = []
  428.  
  429.     for op, sc, n, n_n  in zip(self.op_name, self.score, self.n_abs, self.n_n_abs):
  430.       if n > 0: # observed
  431.         if sc <= thres: # correct
  432.           if op in ops: # in proposed sg
  433.             in_sg_and_seemingly_correct.append( op )
  434.           else: # correct but not in proposed sg
  435.             observed_but_not_in_sg.append( op )
  436.         else: # not correct
  437.           if op not in ops: # not in sg, and not correct
  438.             not_in_sg_and_seemingly_incorrect.append( op )
  439.       else: # not observed
  440.         if op in ops: #in proposed sg
  441.           in_sg_but_no_observations.append( op )
  442.  
  443.     """
  444.     print
  445.     print "in sg", ops
  446.     print "in_sg_and_seemingly_correct", in_sg_and_seemingly_correct
  447.     print "not_in_sg_and_seemingly_incorrect", not_in_sg_and_seemingly_incorrect
  448.     print "in_sg_and_seemingly_incorrect", in_sg_and_seemingly_incorrect
  449.     print "observed_but_not_in_sg", observed_but_not_in_sg
  450.     print "in_sg_but_no_observations", in_sg_but_no_observations
  451.     print
  452.     """
  453.  
  454.     pos = len(in_sg_and_seemingly_correct) + len(not_in_sg_and_seemingly_incorrect)
  455.     neg = len(in_sg_and_seemingly_incorrect) + len(observed_but_not_in_sg)
  456.     abstain = len(in_sg_but_no_observations)
  457.  
  458.     return pos,neg,abstain
  459.  
  460.  
  461.  
  462.   def check_conditions(self):
  463.     table_labels = ('Operator', 'absent under operator\n <I/sigI> (violations)','\nn absent',
  464.                                 'not absent under operator \n <I/sigI> (violations)', '\nn not absent',
  465.                                 'all other reflections \n <I/sigI> (violations)', '\nn compl', '\n Score' )
  466.     for  item in [0]: # absence_class in self.abs_check.absence_classes[ self.sg.group().crystal_system() ]:
  467.       table_rows = []
  468.       for condition in self.abs_check.absence_classes[
  469.               str( sgtbx.space_group_info(group=self.sg.group().build_derived_reflection_intensity_group(False))
  470.                  ) ] : # crystal_system() ]:
  471.         n_abs   = 0
  472.         n_n_abs = 0
  473.         n_tot   = 0
  474.         n_abs_viol   = 0
  475.         n_n_abs_viol = 0
  476.         n_tot_viol   = 0
  477.  
  478.         isi_abs     = 0
  479.         isi_n_abs   = 0
  480.         isi_tot = 0
  481.  
  482.         i_abs     = 0
  483.         i_n_abs   = 0
  484.         i_tot  = 0
  485.  
  486.         score = 0
  487.  
  488.         for hkl, i, sigi in zip(self.miller_array.indices(), self.miller_array.data(), self.miller_array.sigmas() ):
  489.           mc, cc = self.abs_check.check(condition,hkl, return_bool=True)
  490.           if mc: # mask checks out
  491.             if cc: # not absent
  492.               n_n_abs += 1
  493.               isi_n_abs += i/sigi
  494.               i_n_abs   += i
  495.               # should be present. flag if not significant
  496.               if i/sigi < self.cut:
  497.                 n_n_abs_viol += 1
  498.               score += self.score_isigi(i/sigi,absent=False)
  499.             else: #absent
  500.               n_abs += 1
  501.               isi_abs += i/sigi
  502.               i_abs   += i
  503.               # should be absent: flag if significant
  504.               if i/sigi > self.cut:
  505.                 n_abs_viol += 1
  506.               score += self.score_isigi(i/sigi,absent=True)
  507.           else:
  508.             n_tot +=1
  509.             isi_tot += i/sigi
  510.             i_tot += i
  511.             if i/sigi <  self.cut:
  512.               n_tot_viol += 1
  513.         if n_abs > 0:
  514.           isi_abs   = isi_abs/n_abs
  515.           i_abs   = i_abs/n_abs
  516.         if n_n_abs > 0:
  517.           isi_n_abs = isi_n_abs/n_n_abs
  518.           i_n_abs = i_n_abs/n_n_abs
  519.         if n_tot > 0:
  520.           isi_tot   = isi_tot/n_tot
  521.           i_tot   = i_tot/n_tot
  522.  
  523.  
  524.         self.n_abs.append(n_abs)
  525.         self.n_n_abs.append(n_n_abs)
  526.         self.n_tot.append(n_tot)
  527.         self.n_abs_viol.append(n_abs_viol)
  528.         self.n_n_abs_viol.append(n_n_abs_viol)
  529.         self.n_tot_viol.append(n_tot_viol)
  530.  
  531.         self.isi_abs.append(isi_abs)
  532.         self.isi_n_abs.append(isi_n_abs)
  533.         self.isi_tot.append(isi_tot)
  534.  
  535.         self.i_abs.append(i_abs)
  536.         self.i_n_abs.append(i_n_abs)
  537.         self.i_tot.append(i_tot)
  538.  
  539.         self.op_name.append( condition )
  540.         score = score/max(1,n_abs+n_n_abs)
  541.         self.score.append( score )
  542.  
  543.  
  544.         table_rows.append( [condition, str("%8.2f  (%i, %4.1f%s)"%(isi_abs,n_abs_viol, 100.0*float(n_abs_viol)/max(1,n_abs),'%' )),
  545.                                        str("%8.0f"%(n_abs)),
  546.                                        str("%8.2f  (%i, %4.1f%s)"%(isi_n_abs,n_n_abs_viol, 100.0*float(n_n_abs_viol)/max(1,n_n_abs),'%' )),
  547.                                        str("%8.0f"%(n_n_abs)),
  548.                                        str("%8.2f  (%i, %4.1f%s)"%(isi_tot,n_tot_viol, 100.0*float(n_tot_viol)/max(1,n_tot),'%' )),
  549.                                        str("%8.0f"%(n_tot)),
  550.                                        str("%8.2e"%(abs(score)))] )
  551.  
  552.  
  553.       self.table = table_utils.format([table_labels]+table_rows,
  554.                                        comments=None,
  555.                                        has_header=True,
  556.                                        separate_rows=False,
  557.                                        prefix='| ',
  558.                                        postfix=' |')
  559.  
  560.  
  561.  
  562. class sgi_iterator(object):
  563.   def __init__(self,chiral=True, crystal_system=None, intensity_symmetry=None):
  564.     self.chiral = chiral
  565.     self.crystal_system = crystal_system
  566.     self.intensity_symmetry = intensity_symmetry
  567.     assert ( [self.crystal_system , self.intensity_symmetry] ).count(None) != 0
  568.  
  569.   def comparator(self, sgi):
  570.     if self.crystal_system is not None:
  571.       return (self.crystal_system is None) or (self.crystal_system == sgi.group().crystal_system())
  572.     else:
  573.       return (self.intensity_symmetry is None) or (self.intensity_symmetry == sgi.group().build_derived_reflection_intensity_group(False))
  574.  
  575.   def list(self):
  576.     for symbols in sgtbx.space_group_symbol_iterator():
  577.       sgi = sgtbx.space_group_info(group=sgtbx.space_group(space_group_symbols=symbols))
  578.       if self.comparator(sgi):
  579.         if (self.chiral is None) or (self.chiral == sgi.group().is_chiral()):
  580.           yield sgi
  581.  
  582.  
  583. class protein_space_group_choices(object):
  584.   def __init__(self, miller_array,threshold = 3, out=None,protein=True):
  585.     self.out = out
  586.     if self.out is None:
  587.       self.out = sys.stdout
  588.  
  589.     self.threshold = 3.0
  590.     self.miller_array = miller_array.deep_copy().f_sq_as_f().average_bijvoet_mates().f_as_f_sq().map_to_asu()
  591.  
  592.     self.absences_table = analyze_absences(self.miller_array,threshold,self.out)
  593.     self.absences_table.show()
  594.  
  595.  
  596.     self.sg_iterator = sgi_iterator(chiral = True,
  597.                                     intensity_symmetry = self.miller_array.space_group().build_derived_reflection_intensity_group(False) )
  598.  
  599.     self.sg_choices  = []
  600.     self.mean_i      = []
  601.     self.mean_isigi  = []
  602.     self.n           = []
  603.     self.violations  = []
  604.     self.abs_types   = []
  605.     self.tuple_score = []
  606.  
  607.     legend = [('space group', 'n absent', 'violations (n ; %)', '<Z>_absent', '<Z/sigZ>_absent', '+++', '---', '???')]
  608.     rows = []
  609.     score = []
  610.  
  611.  
  612.     for sg in self.sg_iterator.list():
  613.       xs = crystal.symmetry( unit_cell = self.miller_array.unit_cell(),
  614.                              space_group = sg.group() )
  615.       tmp_miller = self.miller_array.customized_copy( crystal_symmetry = xs )
  616.       these_absent_millers = tmp_miller.select(tmp_miller.sys_absent_flags().data() )
  617.  
  618.       if these_absent_millers.data().size() > 0:
  619.         tmp_mean_i = flex.mean( these_absent_millers.data() )
  620.         tmp_mean_isigi = flex.mean( these_absent_millers.data() / these_absent_millers.sigmas() )
  621.         tmp_n = these_absent_millers.data().size()
  622.         tmp_violations = flex.bool( these_absent_millers.data() / these_absent_millers.sigmas() > self.threshold ).count( True )
  623.       else:
  624.         tmp_mean_i = 0
  625.         tmp_mean_isigi = 0
  626.         tmp_n = 0
  627.         tmp_violations = 0
  628.  
  629.       to_be_checked = []
  630.       for s in sg.group():
  631.         #check if this is an operator that causes absences
  632.         tmp =  conditions_for_operator( s )
  633.         if tmp.absence_type() != "None":
  634.           if tmp.absence_type() in self.absences_table.op_name:
  635.             ii = self.absences_table.op_name.index( tmp.absence_type() )
  636.             if tmp.absence_type() not in to_be_checked:
  637.               to_be_checked.append( tmp.absence_type() )
  638.               tmp_score = self.absences_table.score[ ii ]
  639.       self.abs_types.append( to_be_checked )
  640.       tuple_score =  self.absences_table.propose( to_be_checked )
  641.       self.tuple_score.append( tuple_score )
  642.  
  643.       self.sg_choices.append(  sg )
  644.       self.mean_i.append( tmp_mean_i )
  645.       self.mean_isigi.append( tmp_mean_isigi )
  646.       self.n.append( tmp_n )
  647.       self.violations.append( tmp_violations )
  648.       rows.append( [str(sg), '%i'%(tmp_n), '%5i ; %5.1f '%(tmp_violations, 100.0*tmp_violations/max(tmp_n,1)),
  649.                     '%8.2f  '%(tmp_mean_i), '%8.2f  '%(tmp_mean_isigi), ' %i '%(tuple_score[0]), ' %i '%(tuple_score[1]), ' %i '%(tuple_score[2])  ] )
  650.  
  651.     self.table = table_utils.format( legend + rows,
  652.                                      comments=None,
  653.                                      has_header=True,
  654.                                      separate_rows=False,
  655.                                      prefix='| ',
  656.                                      postfix=' |')
  657.     tmp_rows = self.suggest_likely_candidates()
  658.     self.sorted_table = table_utils.format( legend + tmp_rows,
  659.                                      comments=None,
  660.                                      has_header=True,
  661.                                      separate_rows=False,
  662.                                      prefix='| ',
  663.                                      postfix=' |')
  664.     print >> self.out
  665.     print >> self.out, "Analyses of the absences table indicates a number of likely space group candidates,"
  666.     print >> self.out, "which are listed below. for each space group, the number of operators with absence"
  667.     print >> self.out, "classes arguing for the proposed space group are listed under the '+++' column."
  668.     print >> self.out, "The number of absences classes arguing against the proposed group are listed under '---'"
  669.     print >> self.out, "The number of absences associated with an operator in the proposed space group has not"
  670.     print >> self.out, "been found in the data, are listed under '???'."
  671.     print >> self.out
  672.     print >> self.out, self.sorted_table
  673.     print >> self.out
  674.     print >> self.out
  675.  
  676.  
  677.  
  678.   def suggest_likely_candidates( self, acceptable_violations = 1 ):
  679.     used = flex.bool( len(self.sg_choices), False )
  680.     order = []
  681.  
  682.     all_done = False
  683.     count = -1
  684.     while not all_done:
  685.       # find that spacegroup which is not used, has the largest number of absences, and a low level of acceptable violations
  686.       max_in_favour = -10
  687.       max_index = -10
  688.       count += 1
  689.       for ii in xrange( used.size() ):
  690.         if not used[ii]:
  691.           if self.tuple_score[ii][1] <= acceptable_violations :
  692.             if self.tuple_score[ii][0] > max_in_favour:
  693.               max_in_favour = self.tuple_score[ii][0]
  694.               max_index = ii
  695.       if max_in_favour < 0:
  696.         all_done = True
  697.       else:
  698.         used[ max_index ] =True
  699.         order.append( max_index )
  700.       if all_done:
  701.         if len(order)==0:
  702.           all_done=False
  703.           acceptable_violations += 1
  704.  
  705.  
  706.     sorted_rows = []
  707.     for ii in order:
  708.       sg             = self.sg_choices[ii]
  709.       tmp_n          = self.n[ii]
  710.       tmp_violations = self.violations[ii]
  711.       tmp_mean_i     = self.mean_i[ii]
  712.       tmp_mean_isigi = self.mean_isigi[ii]
  713.       tuple_score    = self.tuple_score[ii]
  714.  
  715.       sorted_rows.append( [str(sg), '%i'%(tmp_n), '%5i ; %5.1f '%(tmp_violations, 100.0*tmp_violations/max(tmp_n,1)),
  716.                            '%8.2f  '%(tmp_mean_i), '%8.2f  '%(tmp_mean_isigi),  ' %i '%(tuple_score[0]), ' %i '%(tuple_score[1]), ' %i '%(tuple_score[2])
  717.                           ])
  718.  
  719.     return sorted_rows
  720.  
  721.  
  722.  
  723.  
  724. def test():
  725.   tmp = absences()
  726.   assert tmp.check( "2_1 (c)", (0,0,1),True ) == (True, False)
  727.   assert tmp.check( "2_1 (c)", (0,0,4),True ) == (True, True)
  728.   assert tmp.check( "4_1 (a)", (4,0,0),True ) == (True, True)
  729.   assert tmp.check( "3_1 (c)", (0,0,3),True ) == (True, True)
  730.  
  731.   tmp = sgi_iterator(chiral = True, crystal_system = None,
  732.                      intensity_symmetry = sgtbx.space_group_info( "P222").group().build_derived_reflection_intensity_group(False)  )
  733.   sg_list = []
  734.   abs_list = []
  735.   for sg in tmp.list():
  736.     sg_list.append( str(sg) )
  737.     if str(sg)== "P 21 21 21":
  738.       for s in sg.group():
  739.         abs_list.append( conditions_for_operator( s ).absence_type() )
  740.   assert "2_1 (a)" in abs_list
  741.   assert "2_1 (b)" in abs_list
  742.   assert "2_1 (c)" in abs_list
  743.  
  744.  
  745.   assert "P 2 2 2" in sg_list
  746.   assert "P 21 2 2" in sg_list
  747.   assert "P 2 21 2" in sg_list
  748.   assert "P 2 2 21" in sg_list
  749.   assert "P 21 21 2" in sg_list
  750.   assert "P 21 2 21" in sg_list
  751.   assert "P 2 21 21" in sg_list
  752.   assert "P 21 21 21" in sg_list
  753.  
  754.   """
  755.   tmp = sgi_iterator(chiral=None)
  756.   pg  = []
  757.   asu = []
  758.   for sgi in tmp.list():
  759.     a = str( sgtbx.space_group_info( group = sgi.group().build_derived_reflection_intensity_group(False) ) )
  760.     b = str( sgi.reciprocal_space_asu().reference_as_string() )
  761.     if a not in pg:
  762.       pg.append( a )
  763.       asu.append( b )
  764.   for i,j in zip(pg,asu):
  765.     print i, "   ", j
  766.   """
  767.  
  768.   print "OK"
  769.  
  770. if __name__ == "__main__":
  771.   test()

Paste is for source code and general debugging text.

Login or Register to edit, delete and keep track of your pastes and more.

Raw Paste

Login or Register to edit or fork this paste. It's free.