feat: add combine_phases arg to CIPHERGeometry

cipher_parse/cipher_input.py

@@ -333,6 +333,7 @@ def from_voronoi(
         num_phases=None,
         random_seed=None,
         is_periodic=False,
+        combine_phases=None,
     ):
         geometry = CIPHERGeometry.from_voronoi(
             num_phases=num_phases,
@@ -343,6 +344,7 @@ def from_voronoi(
             size=size,
             random_seed=random_seed,
             is_periodic=is_periodic,
+            combine_phases=combine_phases,
         )
 
         inp = cls(
@@ -366,6 +368,7 @@ def from_seed_voronoi(
         solution_parameters,
         random_seed=None,
         is_periodic=False,
+        combine_phases=None,
     ):
         return cls.from_voronoi(
             seeds=seeds,
@@ -378,6 +381,7 @@ def from_seed_voronoi(
             solution_parameters=solution_parameters,
             random_seed=random_seed,
             is_periodic=is_periodic,
+            combine_phases=combine_phases,
         )
 
     @classmethod
@@ -393,6 +397,7 @@ def from_random_voronoi(
         solution_parameters,
         random_seed=None,
         is_periodic=False,
+        combine_phases=None,
     ):
         return cls.from_voronoi(
             num_phases=num_phases,
@@ -405,6 +410,7 @@ def from_random_voronoi(
             solution_parameters=solution_parameters,
             random_seed=random_seed,
             is_periodic=is_periodic,
+            combine_phases=combine_phases,
         )
 
     @classmethod
@@ -418,13 +424,15 @@ def from_voxel_phase_map(
         outputs,
         solution_parameters,
         random_seed=None,
+        combine_phases=None,
     ):
         geometry = CIPHERGeometry(
             voxel_phase=voxel_phase,
             materials=materials,
             interfaces=interfaces,
             size=size,
             random_seed=random_seed,
+            combine_phases=combine_phases,
         )
         inp = cls(
             geometry=geometry,
@@ -445,6 +453,7 @@ def from_dream3D(
         solution_parameters,
         container_labels=None,
         phase_type_map=None,
+        combine_phases=None,
     ):
         default_container_labels = {
             "SyntheticVolumeDataContainer": "SyntheticVolumeDataContainer",
@@ -559,6 +568,7 @@ def from_dream3D(
             components=components,
             outputs=outputs,
             solution_parameters=solution_parameters,
+            combine_phases=combine_phases,
         )
 
     @property

cipher_parse/geometry.py

@@ -1,3 +1,5 @@
+import random
+
 from damask import Orientation
 import pyvista as pv
 import numpy as np
@@ -39,6 +41,7 @@ def __init__(
         time=None,
         increment=None,
         incremental_data_idx=None,
+        combine_phases=None,
     ):
         """
         Parameters
@@ -107,6 +110,22 @@ def __init__(
 
         self._ensure_phase_assignment(random_seed)
 
+        if combine_phases:
+            # combine multiple phases into the same phase to reduce the computational cost
+            self.voxel_phase = self.combine_phases_per_phase_type(
+                voxel_map,
+                materials,
+                combine_phases,
+            )
+            self.voxel_map = VoxelMap(
+                region_ID=voxel_phase,
+                size=size,
+                is_periodic=is_periodic,
+                quiet=quiet,
+            )
+            all_phases = self.present_phases
+            self._num_phases = all_phases.size
+
         self._phase_material = self._get_phase_material()
         self._validate_interfaces()
         self._check_interface_phase_pairs()
@@ -131,6 +150,93 @@ def __init__(
         self._misorientation_matrix = None
         self._misorientation_matrix_is_degrees = None
 
+    @staticmethod
+    def combine_phases_per_phase_type(voxel_map, materials, combine_phases):
+
+        print(f"combining phases according to {combine_phases}")
+
+        voxel_phase_new = voxel_map.region_ID
+        neighbours = voxel_map.neighbour_list
+
+        for mat in materials:
+            for pt_i in mat.phase_types:
+                print(f"{pt_i.name=}")
+                max_phases_i = combine_phases.get(pt_i.name)
+                if max_phases_i:
+
+                    n_phases = pt_i.phases.size
+                    group_size = int(np.ceil(n_phases / max_phases_i))
+
+                    print(f" n_phases: {n_phases}")
+                    print(f" max_phases_i: {max_phases_i}")
+                    print(f" group_size: {group_size}")
+
+                    voxel_phase_new = np.copy(voxel_phase_new)
+
+                    reassignment = {}
+                    reassignment_inv = {}
+                    possible_IDs_idx = set(range(max_phases_i, n_phases))
+                    possible_IDs = set(pt_i.phases[list(possible_IDs_idx)])
+                    count = 0
+                    kept_IDs = []
+                    for root_ID_idx in range(int(np.ceil(n_phases / group_size))):
+
+                        root_ID = pt_i.phases[root_ID_idx]
+                        kept_IDs.append(root_ID)
+
+                        neighbours_i = set(neighbours[:, neighbours[0] == root_ID][1])
+                        possible_IDs -= {root_ID}
+                        possible_IDs_i = possible_IDs - neighbours_i
+                        shared_IDs = [root_ID]
+                        count += 1
+
+                        for _ in range(1, group_size):
+                            if count >= n_phases:
+                                break
+                            try:
+                                sampled_ID = random.sample(possible_IDs_i, 1)[0]
+                            except ValueError:
+                                print(
+                                    f"No non-neighbouring samples left for root_ID: {root_ID}."
+                                )
+                                # allow touching phase IDs within this group:
+                                possible_IDs_i = possible_IDs - set(shared_IDs)
+
+                                if possible_IDs_i:
+                                    sampled_ID = random.sample(possible_IDs_i, 1)[0]
+                                else:
+                                    raise ValueError(
+                                        f"Cannot find non-neighbouring root IDs"
+                                    ) from None
+
+                            neighbours_sampled = set(
+                                neighbours[:, neighbours[0] == sampled_ID][1]
+                            )
+                            possible_IDs_i -= neighbours_sampled
+                            possible_IDs_i -= {sampled_ID}
+                            possible_IDs -= {sampled_ID}
+                            shared_IDs.append(sampled_ID)
+                            count += 1
+
+                        for i in shared_IDs:
+                            reassignment[i] = root_ID
+                            voxel_phase_new[voxel_phase_new == i] = root_ID
+
+                        reassignment_inv[root_ID] = shared_IDs
+
+                    pt_i.phases = kept_IDs  # modify phase type phases
+
+                    # reindex phases across all materials to maintain consecutive phase IDs:
+                    voxel_phase_new_flat = voxel_phase_new.reshape(-1)
+                    uniq, inv = np.unique(voxel_phase_new_flat, return_inverse=True)
+                    reindex = dict(zip(uniq, range(len(uniq))))
+                    voxel_phase_new = inv.reshape(voxel_phase_new.shape)
+                    for mat_j in materials:
+                        for pt_j in mat_j.phase_types:
+                            pt_j.phases = np.array([reindex[i] for i in pt_j.phases])
+
+        return voxel_phase_new
+
     def __eq__(self, other):
         # Note we don't check seeds (not stored in YAML file)
         if not isinstance(other, self.__class__):
@@ -456,7 +562,7 @@ def _get_phase_material(self):
                 "Not all phases are accounted for in the phase type definitions."
             )  # TODO: test raise
 
-        # check all phase indices form a consequtive range:
+        # check all phase indices form a consecutive range:
         num_phases_range = set(np.arange(self.num_phases))
         known_phases = set(np.hstack(all_phase_idx))
         miss_phase_idx = num_phases_range - known_phases
@@ -740,6 +846,7 @@ def from_voronoi(
         num_phases=None,
         random_seed=None,
         is_periodic=False,
+        combine_phases=None,
     ):
         if sum(i is not None for i in (seeds, num_phases)) != 1:
             raise ValueError(f"Specify exactly one of `seeds` and `num_phases`")
@@ -769,6 +876,7 @@ def from_voronoi(
             size=size,
             seeds=seeds,
             random_seed=random_seed,
+            combine_phases=combine_phases,
         )
 
     @classmethod