dials.integrate overlaps_filter options broken

[dagewa]

Reported on the support list. We don't have a multi lattice data at hand to test (I think?) but the issue can be demonstrated with any old dataset

dials.import $(dials.data get -q insulin)
dials.find_spots imported.expt
dials.index imported.expt strong.refl space_group=I23
dials.integrate indexed.expt indexed.refl\
  integration.overlaps_filter.foreground_foreground.enable=True\
  integration.overlaps_filter.foreground_background.enable=True

Ends with

Traceback (most recent call last):
  File "/home/fcx32934/sw/cctbx/build/../modules/dials/command_line/integrate.py", line 734, in <module>
    run()
  File "/home/fcx32934/sw/cctbx/conda_base/lib/python3.8/contextlib.py", line 75, in inner
    return func(*args, **kwds)
  File "/home/fcx32934/sw/cctbx/build/../modules/dials/command_line/integrate.py", line 712, in run
    experiments, reflections, report = run_integration(
  File "/home/fcx32934/sw/cctbx/build/../modules/dials/command_line/integrate.py", line 585, in run_integration
    reflections = integrator.integrate()
  File "/home/fcx32934/sw/cctbx/modules/dials/algorithms/integration/integrator.py", line 1322, in integrate
    self.reflections, self.experiments = self.finalize_reflections(
  File "/home/fcx32934/sw/cctbx/modules/dials/algorithms/integration/integrator.py", line 524, in _finalize_rotation
    reflections, experiments = _finalize(reflections, experiments, params)
  File "/home/fcx32934/sw/cctbx/modules/dials/algorithms/integration/integrator.py", line 512, in _finalize
    overlaps_filter.remove_foreground_foreground_overlaps()
  File "/home/fcx32934/sw/cctbx/modules/dials/algorithms/integration/overlaps_filter.py", line 160, in remove_foreground_foreground_overlaps
    f.remove_foreground_foreground_overlaps()
  File "/home/fcx32934/sw/cctbx/modules/dials/algorithms/integration/overlaps_filter.py", line 136, in remove_foreground_foreground_overlaps
    self.create_referenced_mask(self.code_fgd, "foreground")
  File "/home/fcx32934/sw/cctbx/modules/dials/algorithms/integration/overlaps_filter.py", line 67, in create_referenced_mask
    shoebox = obs["shoebox"]
KeyError: 'Please report this error to dials-support@lists.sourceforge.net: shoebox'

[graeme-winter]

I'll find you something

[dagewa]

If there are multi-lattice images with overlaps on Zenodo, I think that would make a good addition for dials.data. We probably only need a fairly narrow wedge to test success (or not) of whatever our overlap handling facilities are.

[graeme-winter]

@dagewa what do you know https://zenodo.org/record/10820

Like, we did this 7 years back 😀👍

Fortune favours the prepared

[dagewa]

Nice. Possibly a bit big for dials data though for writing a new test?

[dagewa]

Guess we can take one of the tarballs, split to a smaller one and make a new Zenodo upload of just a small sweep

[dagewa]

I'll look into doing that with a reduced version of semisynthetic_multilattice_data_6.tar.bz2

[dagewa]

Right, I have integrated the 6 lattices in the semisynthetic set acegik, but according to flags, there are no overlaps:

>>> from dials.array_family import flex
>>> rt=flex.reflection_table.from_file("integrated.refl")
>>> rt.get_flags(rt.flags.overlapped_fg).count(True)
0
>>> rt.get_flags(rt.flags.overlapped_bg).count(True)
0

Does anybody understand the overlap handling in DIALS?

[dagewa]

These really should be overlaps, I'm sure:

[dagewa]

Here they are on the image:

So, these are very clearly overlapped reflections, but DIALS does not count them as such.

[graeme-winter]

I for one am starting to smell a little (well justified) scope creep with this issue.

[dagewa]

You do have a point. I was trying to figure out overlaps_filter is supposed to work, so that it could be tested, but so far I remain in the dark.

[dagewa]

Ok, the first problem is that shoeboxes get deleted here

dials/algorithms/integration/processor.py

Lines 561 to 563 in a2cb71b

	# Delete the shoeboxes
	if self.params.debug.separate_files or not self.params.debug.output:
	del self.reflections["shoebox"]

But even if I stop that happening, it goes on to crash:

Traceback (most recent call last):
  File "/home/fcx32934/sw/cctbx/build/../modules/dials/command_line/integrate.py", line 734, in <module>
    run()
  File "/home/fcx32934/sw/cctbx/conda_base/lib/python3.8/contextlib.py", line 75, in inner
    return func(*args, **kwds)
  File "/home/fcx32934/sw/cctbx/build/../modules/dials/command_line/integrate.py", line 712, in run
    experiments, reflections, report = run_integration(
  File "/home/fcx32934/sw/cctbx/build/../modules/dials/command_line/integrate.py", line 585, in run_integration
    reflections = integrator.integrate()
  File "/home/fcx32934/sw/cctbx/modules/dials/algorithms/integration/integrator.py", line 1322, in integrate
    self.reflections, self.experiments = self.finalize_reflections(
  File "/home/fcx32934/sw/cctbx/modules/dials/algorithms/integration/integrator.py", line 524, in _finalize_rotation
    reflections, experiments = _finalize(reflections, experiments, params)
  File "/home/fcx32934/sw/cctbx/modules/dials/algorithms/integration/integrator.py", line 512, in _finalize
    overlaps_filter.remove_foreground_foreground_overlaps()
  File "/home/fcx32934/sw/cctbx/modules/dials/algorithms/integration/overlaps_filter.py", line 160, in remove_foreground_foreground_overlaps
    f.remove_foreground_foreground_overlaps()
  File "/home/fcx32934/sw/cctbx/modules/dials/algorithms/integration/overlaps_filter.py", line 136, in remove_foreground_foreground_overlaps
    self.create_referenced_mask(self.code_fgd, "foreground")
  File "/home/fcx32934/sw/cctbx/modules/dials/algorithms/integration/overlaps_filter.py", line 75, in create_referenced_mask
    if (shoebox.mask[posn_in_shoebox] & test_code) == test_code:
IndexError: Please report this error to dials-support@lists.sourceforge.net: Index out of range.

[stale]

This issue has been automatically marked as stale because it has not had recent activity. The label will be removed automatically if any activity occurs. Thank you for your contributions.

[biochem-fan]

Removing the stale label because I see many datasets with multiple lattices. See also @2555 and @2296.

@dagewa
If it takes time to fix the overlap filter in dials.integrate, how about writing a script to load integrated reflections and remove spots with close predicted xyz? Was this done before? Otherwise I will try myself. The algorithm is very similar to RELION's duplicate particle remover for single particle analysis I wrote several years ago.

[dagewa]

I don't think a script like that was tried before. I won't have time to look at it unfortunately, but I would be interested in your results.

[biochem-fan]

My rudimentary attempt did not improve the map artifact; it can be due to a bug in my script, a bad threshold for duplicates and/or the artifact being caused by some other pathology. At the moment I don't have time to pursue this further, but for the record I put my script here.

from dials.array_family import flex
import numpy as np
import scipy.spatial
import sys

if len(sys.argv) != 4:
    sys.stderr.write("Usage: dials.python remove_overlaps.py input1.refl input2.refl output.refl\n")
    exit(-1)

_, fn_in1, fn_in2, fn_out = sys.argv

d1 = flex.reflection_table.from_file(fn_in1)
d2 = flex.reflection_table.from_file(fn_in2)

a1 = np.array(d1['xyzcal.px'])
a2 = np.array(d2['xyzcal.px'])
t1 = scipy.spatial.KDTree(a1)
t2 = scipy.spatial.KDTree(a2)
dups = t1.query_ball_tree(t2, r = 5)

fine = [len(x) == 0 for x in dups]
good = d1.select(flex.bool(fine))

print("Number of reflections in Input 1:", len(d1))
print("Number of reflections in Input 1:", len(d2))
print("Numer of written reflections:", len(good))
print("Number of overlapped and removed reflections:", len(d1) - len(good))

good.as_file(fn_out)