Rework packing factor in Particle placement. (#762)

* Rework on packing factor in Particle placement.
Add tests and documentation.

* Implement feedback from review.

* Update tests/test_placement.py

---------

Co-authored-by: Robin De Schepper <robin.deschepper93@gmail.com>

bsb/placement/particle.py

@@ -187,6 +187,20 @@ def __init__(self, track_displaced=False, scaffold=None, strat=None):
         self.strat = strat
 
     def fill(self, voxels, particles, check_pack=True):
+        """
+        Fill a list of voxels with Particles.
+
+        :param bsb.voxels.VoxelSet voxels: List of voxels in which to place the particles.
+        :param List[dict] particles: List of dictionary for each particle to place.
+            Each dictionary needs to contain the "name" (str), "radius" (float) of particle.
+            It should also store the "count" of particle (int | List[int]) either in total or
+            for each voxel and "voxels" (List[int]) which gives in which voxel the placement will
+            append.
+        :param bool check_pack: If True, will check the packing factor before placing particles in
+            the voxels.
+        :raise PackingError: If check_pack is True and the resulting packing factor is greater than
+            0.4.
+        """
         # Amount of spatial dimensions
         self.dimensions = voxels.get_raw(copy=False).shape[1]
         # Extend list of particle types in the system
@@ -203,13 +217,13 @@ def fill(self, voxels, particles, check_pack=True):
                 voxels.as_spatial_coords(copy=False), voxels.get_size_matrix(copy=False)
             )
         )
-        pf = self.get_packing_factor()
-        if self.strat is not None:
-            strat_name = type(self.strat).__name__
-        else:
-            strat_name = "particle system"
-        msg = f"Packing factor {round(pf, 2)}"
         if check_pack:
+            pf = self.get_packing_factor()
+            if self.strat is not None:
+                strat_name = type(self.strat).__name__
+            else:
+                strat_name = "particle system"
+            msg = f"Packing factor {round(pf, 2)}"
             if pf > 0.4:
                 if pf > 0.64:
                     msg += " exceeds geometrical maximum packing for spheres (0.64)"
@@ -405,14 +419,29 @@ def deintersect(self, nearest_neighbours=None):
             nearest_neighbours = self.estimate_nearest_neighbours()
 
     def get_packing_factor(self, particles=None, volume=None):
+        """
+        Calculate the packing factor of the volume where particles will be placed.
+        It corresponds to the ratio of the sum of the particles' volume over the volume itself.
+
+        :param List[bsb.placement.particle.Particle] | None particles: List of Particle to place.
+            If None, it will use the ParticleSystem particle_types list.
+        :param float | None volume: Size of the volume in which the particles will be placed.
+            If None, it will use the total volume of the voxels of the ParticleSystem.
+        :return: Packing factor
+        :rtype: float
+        """
         if particles is None:
-            particles_volume = np.sum(
-                [p["count"] * sphere_volume(p["radius"]) for p in self.particle_types]
+            particles_volume = sum(
+                np.sum(p["count"]) * sphere_volume(p["radius"])
+                for p in self.particle_types
             )
         else:
-            particles_volume = np.sum([sphere_volume(p.radius) for p in particles])
+            particles_volume = sum(sphere_volume(p.radius) for p in particles)
         if volume is None:
-            volume = np.sum([np.prod(v.size) for v in self.voxels])
+            volume = sum(
+                sum(np.prod(v.size) for v in np.array(self.voxels)[np.array(p["voxels"])])
+                for p in self.particle_types
+            )
         return particles_volume / volume
 
     def _get_packing_factors(self, particles=None, volume=None):

tests/test_placement.py

@@ -1,3 +1,4 @@
+import re
 import unittest, os, sys, numpy as np, h5py
 
 from bsb import config
@@ -326,7 +327,7 @@ def test_parallel_arrays(self):
         self.assertAll(pos[:, 1] >= cfg.partitions.test_layer.data.ldc[1], "not in layer")
 
 
-class TestVoxelDensities(unittest.TestCase):
+class TestVoxelDensities(RandomStorageFixture, unittest.TestCase, engine_name="hdf5"):
     def test_particle_vd(self):
         cfg = Configuration.default(
             cell_types=dict(
@@ -342,7 +343,7 @@ def test_particle_vd(self):
                 )
             ),
         )
-        network = Scaffold(cfg)
+        network = Scaffold(cfg, self.storage)
         counts = network.placement.voxel_density.get_indicators()["test_cell"].guess(
             chunk=Chunk([0, 0, 0], [100, 100, 100]),
             voxels=network.partitions.test_part.vs,
@@ -354,6 +355,69 @@ def test_particle_vd(self):
         self.assertGreater(len(ps), 90)
         self.assertLess(len(ps), 130)
 
+    def _config_packing_fact(self):
+        return Configuration.default(
+            network={
+                "x": 20.0,
+                "y": 5.0,
+                "z": 20.0,
+                "chunk_size": [20, 20, 10],  # at least two chunks
+            },
+            partitions={
+                "first_layer": {"thickness": 5.0, "stack_index": 0},
+            },
+            cell_types=dict(
+                test_cell=dict(spatial=dict(radius=1.5, count=100)),
+                test_cell2=dict(
+                    spatial=dict(radius=2, relative_to="test_cell", count_ratio=0.05)
+                ),
+            ),
+            placement=dict(
+                test_place2=dict(
+                    strategy="bsb.placement.RandomPlacement",
+                    partitions=["first_layer"],
+                    cell_types=["test_cell2"],
+                ),
+                ch4_c25=dict(
+                    strategy="bsb.placement.ParticlePlacement",
+                    partitions=["first_layer"],
+                    cell_types=["test_cell"],
+                ),
+            ),
+        )
+
+    def test_packing_factor_error1(self):
+        cfg = self._config_packing_fact()
+        network = Scaffold(cfg, self.storage)
+        with self.assertRaises(PackingError) as exc:
+            network.compile(clear=True)
+        self.assertTrue(
+            re.match(
+                r"Packing factor .* exceeds geometrical maximum packing for spheres \(0\.64\).*",
+                str(exc.exception),
+            )
+        )
+
+    def test_packing_factor_error2(self):
+        cfg = self._config_packing_fact()
+        cfg.cell_types["test_cell"] = dict(spatial=dict(radius=1.3, count=100))
+        network = Scaffold(cfg, self.storage)
+        with self.assertRaises(PackingError) as exc:
+            network.compile(clear=True)
+        self.assertTrue(
+            re.match(
+                r"Packing factor .* too high to resolve with ParticlePlacement.*",
+                str(exc.exception),
+            )
+        )
+
+    def test_packing_factor_warning(self):
+        cfg = self._config_packing_fact()
+        cfg.cell_types["test_cell"] = dict(spatial=dict(radius=1, count=100))
+        network = Scaffold(cfg, self.storage)
+        with self.assertWarns(PackingWarning):
+            network.compile(clear=True)
+
 
 class VoxelParticleTest(Partition, classmap_entry="test"):
     vs = VoxelSet(