dials / dials

Diffraction Integration for Advanced Light Sources
https://dials.github.io
BSD 3-Clause "New" or "Revised" License
71 stars 49 forks source link

`dials.index` writes out reflections for a rejected experiment #2652

Open dagewa opened 5 months ago

dagewa commented 5 months ago

I found a weird error trying to process data from https://zenodo.org/records/10974780.

Taking just the first data set (exp_705), which has more than one lattice, then running these commands:

dials.import exp_705/frames/*.rodhypix geometry.goniometer.axis=0.0488498,-0.998806,0
dials.find_spots imported.expt d_max=10
dials.index imported.expt strong.refl\
    detector.fix=distance space_group=P212121\
    max_lattices=3 minimum_angular_separation=1
dials.integrate indexed.expt indexed.refl

leads to an arcane error message like

dials Internal Error: /home/fcx32934/sw/cctbx/modules/dials/src/dials/array_family/boost_python/flex_reflection_table.cc(510): DIALS_ASSERT(id[i] < num_expr) failure.

I think this is related to the fact that we try to find 3 lattices, but the third is rejected as being too similar to the first. However, I have been unable so far to reproduce this with a different data set.

dagewa commented 5 months ago

It looks like this is the problem:

>>> from dials.array_family import flex
>>> rt=flex.reflection_table.from_file("indexed.refl")
>>> list(set(rt["id"]))
[0, 1, 2, -1]
>>> (rt["id"] == 0).count(True)
1329
>>> (rt["id"] == 1).count(True)
1924
>>> (rt["id"] == 2).count(True)
349
>>> (rt["id"] == -1).count(True)
467

There are 349 reflections with an id of 2, but dials.index only wrote out 2 experiments:

>>> from dxtbx.model.experiment_list import ExperimentList
>>> el=ExperimentList.from_file("indexed.expt")
>>> len(el)
2
dagewa commented 5 months ago

Indeed, if I do rt.del_selected(rt["id"] == 2) then save the reflections, then integration works fine.

So, this is a dials.index problem

dagewa commented 5 months ago

Ok, I think I understand why this is a corner case now. For this data set, when we get here: https://github.com/dials/dials/blob/9af1ec1261392ac36abe27ad4a11985e055352cd/src/dials/algorithms/indexing/indexer.py#L614-L616 the value of the loop counter i_cycle is 1. That is, we have already done one round of refinement before the crystal model becomes "too similar" to an existing model. This means that the attribute self.refined_reflections has already been set in the previous round with id values up to 2.

The fix would be to reset these reflections to have an id of -1 just before breaking out of the loop here.

dagewa commented 5 months ago

For the other data set I looked at trying to reproduce this, the new crystal model was immediately "too similar" to an existing model and no refinement with the new model was done. So self.refined_reflections did not contain the id of the new model.