Implement multi-sweep joint indexing and add tests

src/dials/command_line/index.py

@@ -218,7 +218,7 @@ def _index_joint_indexing(experiments, reflections, params):
             input_expts = experiments
             known_crystal_models = None
 
-    for id_ in sorted(set(reflections["id"]), reverse=True):
+    for id_ in sorted(set(reflections["id"]).difference({-1}), reverse=True):
         del reflections.experiment_identifiers()[id_]
     reflections["original_id"] = copy.deepcopy(reflections["imageset_id"])
     reflections["id"] = copy.deepcopy(reflections["imageset_id"])
@@ -245,18 +245,42 @@ def _index_joint_indexing(experiments, reflections, params):
         unindexed_refl = idxr.unindexed_reflections.select(sel)
         unindexed_refl["id"] = flex.int(unindexed_refl.size(), n_id)
         del unindexed_refl["original_id"]
-        unindexed_refl.experiment_identifiers()[n_id] == identifier
+        unindexed_refl.clean_experiment_identifiers_map()
+        unindexed_refl.experiment_identifiers()[n_id] = identifier
         n_id += 1
         indexed_reflections.extend(unindexed_refl)
+        indexed_experiments.append(unindexed[i])
         # now get the indexed
         i_expts = idxr.refined_experiments.where(imageset=iset)
         identifiers = [idxr.refined_experiments[i].identifier for i in i_expts]
+        n_indexed_this = len(identifiers)
         refls = idxr.refined_reflections.select_on_experiment_identifiers(identifiers)
+        refls.reset_ids()  # number from 0
         # now reset the ids in the refls
+        for id_ in sorted(set(refls["id"]), reverse=True):
+            identifier = refls.experiment_identifiers()[id_]
+            del refls.experiment_identifiers()[id_]
+            refls.experiment_identifiers()[id_ + n_id] = identifier
+        refls["id"] += n_id
+        n_id += n_indexed_this
+        del refls["original_id"]
+        indexed_reflections.extend(refls)
+        # the refiner copies the input beam, detector and gonio, we want to share these
+        indexed_experiments[-1].detector = idxr.refined_experiments[i_expts[0]].detector
+        indexed_experiments[-1].beam = idxr.refined_experiments[i_expts[0]].beam
+        if idxr.refined_experiments[
+            i_expts[0]
+        ].goniometer:  # might be using the stills indexer which deletes the gonio
+            indexed_experiments[-1].goniometer = idxr.refined_experiments[
+                i_expts[0]
+            ].goniometer
+        for j in i_expts:
+            indexed_experiments.append(idxr.refined_experiments[j])
+    indexed_reflections.assert_experiment_identifiers_are_consistent(
+        indexed_experiments
+    )
 
-    pass
-    # FIXME
-    return indexed_experiments, refls
+    return indexed_experiments, indexed_reflections
 
 
 def _index_experiments(

tests/algorithms/indexing/test_index.py

@@ -712,17 +712,48 @@ def test_index_multi_lattice_multi_sweep(dials_data, tmp_path):
             "max_lattices=2",
             "joint_indexing=False",
             "n_macro_cycles=2",
+            "output.retain_unindexed_experiment=False",
+            "output.experiments=indexed_1.expt",
+            "output.reflections=indexed_1.refl",
         ],
         cwd=tmp_path,
         capture_output=True,
     )
     assert not result.returncode and not result.stderr
-    assert (tmp_path / "indexed.refl").is_file()
-    assert (tmp_path / "indexed.expt").is_file()
-    expts = load.experiment_list(tmp_path / "indexed.expt", check_format=False)
+    assert (tmp_path / "indexed_1.refl").is_file()
+    assert (tmp_path / "indexed_1.expt").is_file()
+    expts = load.experiment_list(tmp_path / "indexed_1.expt", check_format=False)
     assert len(expts) == 4
     assert len(expts.crystals()) == 4
-    refls = flex.reflection_table.from_file(tmp_path / "indexed.refl")
+    refls = flex.reflection_table.from_file(tmp_path / "indexed_1.refl")
+    refls.assert_experiment_identifiers_are_consistent(expts)
+    assert set(refls["imageset_id"]) == {0, 1}
+    assert refls.get_flags(refls.flags.indexed).count(True) >= n_indexed_first
+
+    # now run with keeping unindexed experiment existing model
+    result = subprocess.run(
+        [
+            shutil.which("dials.index"),
+            tmp_path / "indexed.expt",
+            tmp_path / "indexed.refl",
+            "max_lattices=2",
+            "joint_indexing=False",
+            "n_macro_cycles=2",
+            "output.retain_unindexed_experiment=True",
+            "output.experiments=indexed_2.expt",
+            "output.reflections=indexed_2.refl",
+        ],
+        cwd=tmp_path,
+        capture_output=True,
+    )
+    assert not result.returncode and not result.stderr
+    assert (tmp_path / "indexed_2.refl").is_file()
+    assert (tmp_path / "indexed_2.expt").is_file()
+    expts = load.experiment_list(tmp_path / "indexed_2.expt", check_format=False)
+    assert len(expts) == 6
+    assert expts.crystals()[0] is None
+    assert len(expts.crystals()) == 5
+    refls = flex.reflection_table.from_file(tmp_path / "indexed_2.refl")
     refls.assert_experiment_identifiers_are_consistent(expts)
     assert set(refls["imageset_id"]) == {0, 1}
     assert refls.get_flags(refls.flags.indexed).count(True) >= n_indexed_first
@@ -993,16 +1024,45 @@ def test_multi_lattice_multi_sweep_joint(dials_data, tmp_path):
             dials_data("l_cysteine_dials_output", pathlib=True) / "indexed.refl",
             "max_lattices=2",
             "minimum_angular_separation=0.001",
+            "output.retain_unindexed_experiment=False",
+            "output.experiments=indexed_1.expt",
+            "output.reflections=indexed_1.refl",
         ],
         cwd=tmp_path,
         capture_output=True,
     )
     assert not result.returncode and not result.stderr
-    assert (tmp_path / "indexed.expt").is_file()
-    assert (tmp_path / "indexed.refl").is_file()
+    assert (tmp_path / "indexed_1.expt").is_file()
+    assert (tmp_path / "indexed_1.refl").is_file()
 
-    expts = ExperimentList.from_file(tmp_path / "indexed.expt", check_format=False)
-    refls = flex.reflection_table.from_file(tmp_path / "indexed.refl")
+    expts = ExperimentList.from_file(tmp_path / "indexed_1.expt", check_format=False)
+    refls = flex.reflection_table.from_file(tmp_path / "indexed_1.refl")
     assert len(expts) == 8
     assert len(expts.crystals()) == 2
     refls.assert_experiment_identifiers_are_consistent(expts)
+
+    # now rerun but keeping unindexed experiments
+    result = subprocess.run(
+        [
+            shutil.which("dials.index"),
+            dials_data("l_cysteine_dials_output", pathlib=True) / "indexed.expt",
+            dials_data("l_cysteine_dials_output", pathlib=True) / "indexed.refl",
+            "max_lattices=2",
+            "minimum_angular_separation=0.001",
+            "output.retain_unindexed_experiment=True",
+            "output.experiments=indexed_2.expt",
+            "output.reflections=indexed_2.refl",
+        ],
+        cwd=tmp_path,
+        capture_output=True,
+    )
+    assert not result.returncode and not result.stderr
+    assert (tmp_path / "indexed_2.expt").is_file()
+    assert (tmp_path / "indexed_2.refl").is_file()
+
+    expts = ExperimentList.from_file(tmp_path / "indexed_2.expt", check_format=False)
+    refls = flex.reflection_table.from_file(tmp_path / "indexed_2.refl")
+    assert len(expts) == 12
+    assert expts.crystals()[0] is None
+    assert len(expts.crystals()) == 3  # None and two crystals
+    refls.assert_experiment_identifiers_are_consistent(expts)