consistent asu diffmap

rsbooster/diffmaps/diffmap.py

@@ -70,7 +70,17 @@ def parse_arguments():
 
     return parser#.parse_args()
 
-
+def remove_HKL_duplicates(ds, ds_name):
+    num_duplicates = np.sum(ds.index.duplicated())
+    if num_duplicates > 0:
+        print(f"Warning: {ds_name} contains {num_duplicates} sets of duplicate Miller indices.")
+        print( "Only the first instance of each HKL will be retained.")
+        # useful diagnostic:
+        # print(ds.reset_index().loc[ds.index.duplicated(keep=False),:].sort_values(["H","K","L"]).head(6))
+        ds=ds.loc[~ds.index.duplicated(keep='first'),:]
+    ds.merged=True
+    return ds
+
 def main():
 
     # Parse commandline arguments
@@ -91,6 +101,14 @@ def main():
             f"{args.Phi} is not a phases column in {args.mtz2}. Try again."
         )
 
+    onmtz=onmtz.hkl_to_asu()
+    offmtz=offmtz.hkl_to_asu()
+    ref=ref.hkl_to_asu()
+
+    onmtz = remove_HKL_duplicates(onmtz,  'ON MTZ')
+    offmtz= remove_HKL_duplicates(offmtz, 'OFF MTZ')
+    ref   = remove_HKL_duplicates(ref,    'REF MTZ')
+
     diff = onmtz.merge(offmtz, on=["H", "K", "L"], suffixes=("_on", "_off"))
     diff["DF"] = diff["F_on"] - diff["F_off"]
     diff["SigDF"] = np.sqrt((diff["SigF_on"] ** 2) + (diff["SigF_off"] ** 2))
@@ -108,6 +126,7 @@ def main():
     # Useful for PyMOL
     diff["wDF"] = (diff["DF"] * diff["W"]).astype("SFAmplitude")
 
+
     if args.dmax is None:
         diff.write_mtz(args.outfile)
     else:

rsbooster/diffmaps/internaldiffmap.py

@@ -59,6 +59,13 @@ def parse_arguments():
         default=0.0,
         help="alpha value for computing difference map weights (default=0.0)",
     )
+    parser.add_argument(
+        "-u",
+        "--union",
+        # default=False,
+        action="store_true",
+        help="whether to return an MTZ containing also unmatched pairs of reflections (with DF, wDF=0)",
+    )
     parser.add_argument(
         "-d",
         "--dmax",
@@ -77,6 +84,16 @@ def parse_arguments():
 
     return parser#.parse_args()
 
+def remove_HKL_duplicates(ds, ds_name):
+    num_duplicates = np.sum(ds.index.duplicated())
+    if num_duplicates > 0:
+        print(f"Warning: {ds_name} contains {num_duplicates} sets of duplicate Miller indices.")
+        print( "Only the first instance of each HKL will be retained.")
+        # useful diagnostic:
+        # print(ds.reset_index().loc[ds.index.duplicated(keep=False),:].sort_values(["H","K","L"]).head(6))
+        ds=ds.loc[~ds.index.duplicated(keep='first'),:]
+    ds.merged=True
+    return ds
 
 def main():
 
@@ -88,7 +105,12 @@ def main():
     mtz = subset_to_FSigF(
         *args.inputmtz, {args.inputmtz[1]: "F", args.inputmtz[2]: "SigF"}
     )
+    mtz = mtz.hkl_to_asu()
+    mtz = remove_HKL_duplicates(mtz, 'MTZ of Interest')
+
     ref = rs.read_mtz(refmtz)
+    ref = ref.hkl_to_asu()
+    ref = remove_HKL_duplicates(ref, 'REF MTZ')
 
     # Canonicalize column names
     ref.rename(columns={phi_col: "Phi"}, inplace=True)
@@ -106,10 +128,15 @@ def main():
         sg = gemmi.SpaceGroup(args.spacegroup)
         op = sg.operations().sym_ops[isym]
     except ValueError:
-        op = gemmi.Operation(args.symop)
+        op = gemmi.Op(args.symop)
+    print(f"Applying sym op {op}")
 
+    merge_how='inner'
+    if args.union:
+        merge_how='outer'
+
     internal = mtz.merge(
-        mtz.apply_symop(op).hkl_to_asu(), on=["H", "K", "L"], suffixes=("1", "2")
+        mtz.apply_symop(op).hkl_to_asu(), how=merge_how, on=["H", "K", "L"], suffixes=("1", "2"), 
     )
     internal["DF"] = internal["F1"] - internal["F2"]
     internal["SigDF"] = np.sqrt((internal["SigF1"] ** 2) + (internal["SigF2"] ** 2))
@@ -127,6 +154,11 @@ def main():
     # Useful for PyMOL
     internal["wDF"] = (internal["DF"] * internal["W"]).astype("SFAmplitude")
 
+    if args.union:
+        internal["wDF"  ] = internal["wDF"  ].fillna(0)
+        internal["DF"   ] = internal[ "DF"  ].fillna(0)
+        internal["SigDF"] = internal["SigDF"].fillna(0)
+
     if args.dmax is None:
         internal.write_mtz(args.outfile)
     else:

rsbooster/utils/io.py

@@ -39,6 +39,7 @@ def subset_to_FSigF(mtzpath, data_col, sig_col, column_names_dict={}):
 
     # Run French-Wilson if intensities are provided
     if isinstance(mtz[data_col].dtype, rs.IntensityDtype):
+        print(f"Applying French-Wilson scaling to columns {data_col} and {sig_col}.")
         scaled = rs.algorithms.scale_merged_intensities(
             mtz, data_col, sig_col, mean_intensity_method="anisotropic"
         )