Optimizations 11 2025

app/geometry/geometry_math.py

@@ -49,6 +49,21 @@ def find_coordinates(indexes, voxel_size):
 
     @staticmethod
     def calculate_filled_volume(voxel_space, voxel_size):
+        # Accept either a raw ndarray or a VoxelSpace-like object with a
+        # `_filled_voxels_count` running counter. Prefer the running counter
+        # when available to avoid expensive full-array scans.
+        if hasattr(voxel_space, "_filled_voxels_count"):
+            return int(voxel_space._filled_voxels_count) * voxel_size ** 3
+
+        # If an object with `.space` is provided, try to use its counter first.
+        if hasattr(voxel_space, "space") and hasattr(voxel_space.space, "shape"):
+            # If the container itself tracks a counter, use it.
+            if hasattr(voxel_space, "_filled_voxels_count"):
+                return int(voxel_space._filled_voxels_count) * voxel_size ** 3
+            # Fallback to scanning the ndarray
+            return np.count_nonzero(voxel_space.space) * voxel_size ** 3
+
+        # Otherwise assume voxel_space is a numpy array
         return np.count_nonzero(voxel_space) * voxel_size**3
 
     def find_index(coordinate, voxel_size):

app/geometry/sphere.py

@@ -21,7 +21,10 @@ def deposit_sphere(
     ):
         initial_radius = self.estimate_initial_radius(sphere_volume)
 
-        _, voxel_space = BisectionMethod().execute(
+        # `voxel_space` is expected to be a VoxelSpace object here. The
+        # bisection solver and inner functions will mutate and ultimately
+        # return that object so callers can read `.space` and counters.
+        _, voxel_space_out = BisectionMethod().execute(
             self._deposit_sphere,
             initial_point=initial_radius,
             tolerance=solver_tolerance,
@@ -34,24 +37,56 @@ def deposit_sphere(
             ),
         )
 
-        return voxel_space
+        return voxel_space_out
+
+    def fill_voxels(self, voxel_space_obj, radius, lower_indexes, upper_indexes):
+        # Accept either a VoxelSpace-like object (has `.space`) or a raw ndarray.
+        is_voxel_space = hasattr(voxel_space_obj, "space")
+        space = voxel_space_obj.space if is_voxel_space else voxel_space_obj
 
-    def fill_voxels(self, voxel_space, radius, lower_indexes, upper_indexes):
         empty_voxels = GeometryMath.find_empty_voxels_in_space(
-            voxel_space, lower_indexes, upper_indexes
+            space, lower_indexes, upper_indexes
         )
 
-        for voxel in empty_voxels:
-            voxel_coordinate = GeometryMath.find_coordinates(voxel, self.voxel_size)
-
-            distance_to_centre = GeometryMath.distance(
-                voxel_coordinate, self.centre_coordinates
-            )
-
-            if distance_to_centre <= radius + self.voxel_size * 1e-8:
-                voxel_space[tuple(voxel)] = 1
-
-        return voxel_space
+        if not empty_voxels:
+            return voxel_space_obj
+
+        # Convert to numpy array for vectorized operations
+        empty_voxels_np = np.array(empty_voxels)
+        # Calculate coordinates for all voxels at once
+        voxel_coords = self.voxel_size * (2 * (empty_voxels_np + 1) - 1) * 0.5
+        # Calculate distances to centre for all voxels
+        centre = np.array(self.centre_coordinates)
+        dists = np.linalg.norm(voxel_coords - centre, axis=1)
+        # Mask for voxels within radius
+        mask = dists <= radius + self.voxel_size * 1e-8
+
+        if not np.any(mask):
+            return voxel_space_obj
+
+        # Filter indices that should be filled
+        fill_indices = empty_voxels_np[mask]
+
+        # Advanced index assignment (vectorized)
+        xi = fill_indices[:, 0].astype(int)
+        yj = fill_indices[:, 1].astype(int)
+        zk = fill_indices[:, 2].astype(int)
+
+        # Before setting, count how many of these are actually zero (defensive)
+        current_vals = space[xi, yj, zk]
+        new_mask = current_vals == 0
+        n_new = int(np.count_nonzero(new_mask))
+        if n_new > 0:
+            # Assign into the ndarray `space` (in-place)
+            space[xi[new_mask], yj[new_mask], zk[new_mask]] = 1
+            # Update running counter only when we were given a VoxelSpace object
+            if is_voxel_space and hasattr(voxel_space_obj, "_filled_voxels_count"):
+                voxel_space_obj._filled_voxels_count += n_new
+
+        # Return the same type we were given
+        if is_voxel_space:
+            return voxel_space_obj
+        return space
 
     def estimate_initial_radius(self, volume):
         return (3.0 * volume / (4.0 * math.pi)) ** (1.0 / 3.0)
@@ -85,56 +120,62 @@ def find_sphere_limits(self, radius, nozzle_height):
     """
 
     def deform_voxel_space_for_big_spheres(self, voxel_space, radius):
+        # Accept either a VoxelSpace object or a raw ndarray and return the same type.
         max_indexes = [
             self._find_index(coord + radius) for coord in self.centre_coordinates
         ]
 
+        vs = voxel_space
         for axis_number in range(0, 3):
-            voxel_space = self._maybe_expand_voxel_space(
-                voxel_space, max_indexes[axis_number], axis_number
-            )
+            vs = self._maybe_expand_voxel_space(vs, max_indexes[axis_number], axis_number)
 
-        return voxel_space
+        return vs
 
-    def _maybe_expand_voxel_space(self, voxel_space, max_index, axis_number):
-        size = voxel_space.shape
+    def _maybe_expand_voxel_space(self, voxel_space_obj, max_index, axis_number):
+        # Accept either a VoxelSpace-like object (has `.space`) or a raw ndarray.
+        is_voxel_space = hasattr(voxel_space_obj, "space")
+        space = voxel_space_obj.space if is_voxel_space else voxel_space_obj
 
+        size = space.shape
         index_size = size[axis_number]
 
         if max_index < index_size:
-            return voxel_space
+            return voxel_space_obj
 
         number_to_be_added = max_index - index_size + 1
 
-        new_size = list(size)
-        new_size[axis_number] = number_to_be_added
+        # Add a buffer to reduce the number of reallocations. Buffer is 20%
+        # of current size (rounded), but at least the minimum required.
+        buffer_layers = max(int(index_size * 0.2), 1)
+        layers_to_add = max(number_to_be_added, buffer_layers)
 
-        mat_add = np.zeros(new_size, dtype=np.int8)
+        # Create only the additional block to concatenate
+        mat_add_size = list(size)
+        mat_add_size[axis_number] = layers_to_add
+        mat_add = np.zeros(mat_add_size, dtype=np.int8)
 
-        voxel_space = np.concatenate((voxel_space, mat_add), axis=axis_number)
+        new_space = np.concatenate((space, mat_add), axis=axis_number)
 
-        return voxel_space
+        if is_voxel_space:
+            voxel_space_obj.space = new_space
+            return voxel_space_obj
 
-    def _deposit_sphere(self, radius, voxel_space, nozzle_height, target_volume):
-        copy_voxel_space = voxel_space.copy()
+        return new_space
 
-        copy_voxel_space = self.deform_voxel_space_for_big_spheres(
-            copy_voxel_space, radius
-        )
+    def _deposit_sphere(self, radius, voxel_space, nozzle_height, target_volume):
+        # Accept and return either a VoxelSpace object or a raw ndarray.
+        vs = self.deform_voxel_space_for_big_spheres(voxel_space, radius)
 
         lower_indexes, upper_indexes = self.find_sphere_limits(radius, nozzle_height)
 
-        copy_voxel_space = self.fill_voxels(
-            copy_voxel_space, radius, lower_indexes, upper_indexes
-        )
+        vs = self.fill_voxels(vs, radius, lower_indexes, upper_indexes)
 
-        current_volume = GeometryMath.calculate_filled_volume(
-            copy_voxel_space, self.voxel_size
-        )
+        current_volume = GeometryMath.calculate_filled_volume(vs, self.voxel_size)
 
         volume_overshoot = current_volume / target_volume - 1.0
 
-        return volume_overshoot, copy_voxel_space
+        # Return the computed overshoot and the (possibly mutated) voxel container
+        return volume_overshoot, vs
 
     def _increase_solver_tolerance(self, radius_a, radius_b):
         return radius_b - radius_a < self.voxel_size * 0.5

app/geometry/voxel_space.py

@@ -45,12 +45,54 @@ def __init__(
         self._simulation = simulation_config
         self._printer = printer
         self._consider_acceleration = self._simulation.consider_acceleration
+        # Running count of filled voxels (to avoid repeated `np.count_nonzero` calls)
+        self._filled_voxels_count = 0
 
     def initialize_space(self):
+        # Try to intelligently preallocate Z dimension to avoid repeated expansions.
+        # Estimate the maximum per-step sphere radius from filaments and add a safety margin.
+        z_dim = self.dimensions["z"]
+
+        try:
+            max_radius = 0.0
+            step_size = self._simulation.step_size
+            for filament in self._instruction.filaments_coordinates:
+                coord_old = filament[0]
+                coord_new = filament[1]
+                volume = filament[2]
+
+                # compute filament length and estimated number of steps
+                length = GeometryMath.distance(coord_old, coord_new)
+                n_steps = int(round(length / step_size))
+                if n_steps < 1:
+                    n_steps = 1
+
+                per_step_volume = volume / n_steps
+                # estimate radius of sphere that would contain this volume
+                radius = (3.0 * per_step_volume / (4.0 * math.pi)) ** (1.0 / 3.0)
+                if radius > max_radius:
+                    max_radius = radius
+
+            # include configured sphere z offset
+            max_extra_z = int(math.ceil((max_radius + self._simulation.sphere_z_offset) / self._simulation.voxel_size))
+
+            # safety margin (20% of current z) to reduce chance of further expansion
+            safety_margin = int(max(1, z_dim * 0.2))
+
+            if max_extra_z > 0:
+                z_dim = z_dim + max_extra_z + safety_margin
+
+            logging.getLogger(__name__).info(f"Preallocating voxel space z-dimension: base={self.dimensions['z']}, extra={max_extra_z}, safety={safety_margin}, final={z_dim}")
+        except Exception:
+            # Fallback to original allocation if any issue occurs during estimation
+            z_dim = self.dimensions["z"]
+
         self.space = np.zeros(
-            (self.dimensions["x"], self.dimensions["y"], self.dimensions["z"]),
+            (self.dimensions["x"], self.dimensions["y"], z_dim),
             dtype=np.int8,
         )
+        # reset the running counter when allocating space
+        self._filled_voxels_count = 0
 
     def print(self):
         number_printed_filaments = 0
@@ -135,8 +177,9 @@ def _deposit_filament(
         filament_initial_coordinates,
         volumes,
     ):
-        total_deposited_volume = GeometryMath.calculate_filled_volume(
-            self.space, self._simulation.voxel_size
+        # Use the running counter to compute total deposited volume quickly
+        total_deposited_volume = (
+            self._filled_voxels_count * self._simulation.voxel_size ** 3
         )
 
         for step_n in range(0, number_simulation_steps):
@@ -183,11 +226,20 @@ def _deposit_filament(
                 f"Depositing filament: step = {step_n + 1}/{number_simulation_steps}"
             )
 
-            self.space = sphere.deposit_sphere(
-                voxel_space=self.space,
+            # Pass the VoxelSpace instance so sphere code can update the
+            # running counter and expand the `.space` ndarray in-place.
+            voxel_space_out = sphere.deposit_sphere(
+                voxel_space=self,
                 nozzle_height=nozzle_height,
                 sphere_volume=sphere_volume,
                 voxel_space_target_volume=volume_target,
                 solver_tolerance=self._simulation.solver_tolerance,
                 radius_increment=self._simulation.radius_increment,
             )
+
+            # Ensure our internal `.space` and counter reflect any mutations
+            # returned by the sphere logic (voxel_space_out is a VoxelSpace)
+            if hasattr(voxel_space_out, "space"):
+                self.space = voxel_space_out.space
+            if hasattr(voxel_space_out, "_filled_voxels_count"):
+                self._filled_voxels_count = voxel_space_out._filled_voxels_count

app/physics/volume.py

@@ -27,7 +27,7 @@ def get_volumes_for_filament(
         if not consider_acceleration:
             # Simple case: evenly distribute the volume across steps
             volume_per_step = total_volume / number_simulation_steps
-            return [volume_per_step * (i + 1) for i in range(number_simulation_steps)]
+            return [volume_per_step for i in range(number_simulation_steps)]
 
         # Complex case: use acceleration-specific implementation
         from app.physics.acceleration.volume import AccelerationVolume

app/solvers/bisection_method.py

@@ -23,35 +23,43 @@ def execute(
         else:
             point_b = initial_point
 
-        _, _, point_a, point_b = self._loop_fb(fun, point_b, increment, args)
+        fb, out, point_a, point_b, updated_args = self._loop_fb(fun, point_b, increment, args)
 
-        return self._loop_fc(
-            fun, point_a, point_b, tolerance, fun_increase_tolerance, args
+        # _loop_fc returns (fc, out) where out may be the raw output from `fun`.
+        fc, out = self._loop_fc(
+            fun, point_a, point_b, tolerance, fun_increase_tolerance, updated_args
         )
 
+        # If `out` is a tuple like (volume_overshoot, voxel_space), unwrap it
+        if isinstance(out, tuple) and len(out) == 2:
+            return fc, out[1]
+
+        return fc, out
+
     def _loop_fb(self, fun, point_b, inc, args):
         fb = -1.0
-
         point_a = 0.0
-
+        updated_args = list(args)
         while fb < 0.0:
-
-            fb, out = fun(point_b, *args)
-
+            fb, out = fun(point_b, *updated_args)
+            # If the function returns an updated voxel_space, propagate it.
+            # The returned `out` may be an ndarray (has `shape`) or a
+            # VoxelSpace-like object with a `space` attribute. Accept both.
+            if hasattr(out, "shape") or hasattr(out, "space"):
+                updated_args[0] = out
             if fb < 0:
                 point_a = point_b
                 point_b += inc
-
-        return fb, out, point_a, point_b
+        return fb, out, point_a, point_b, updated_args
 
     def _loop_fc(self, fun, point_a, point_b, tolerance, fun_increase_tolerance, args):
         fc = 2.0 * tolerance
-
+        updated_args = list(args)
         while abs(fc) > tolerance:
             point_c = (point_a + point_b) * 0.5
-
-            fc, out = fun(point_c, *args)
-
+            fc, out = fun(point_c, *updated_args)
+            if hasattr(out, "shape") or hasattr(out, "space"):
+                updated_args[0] = out
             if fun_increase_tolerance(point_a, point_b):
                 logger.debug(
                     "[BisectionMethod]: increasing tolerance because point_b and point_a are too close"

tests/geometry/voxel_space_test.py

@@ -13,6 +13,10 @@ def __init__(self):
         self._number_printed_filaments = 0
         self._filaments_coordinates = list()
         self._default_nozzle_speed = 10.0
+        # Patch: provide a _printer attribute for compatibility
+        class _Printer:
+            feedstock_filament_diameter = 1.75
+        self._printer = _Printer()
 
     def read(self):
         pass
@@ -47,6 +51,7 @@ def __init__(self):
         self.step_size = 0.2
         self.radius_increment = 0.1
         self.solver_tolerance = 0.001
+        self.consider_acceleration = False
 
 
 class FakePrinter:
@@ -114,17 +119,19 @@ def __init__(self):
         self.step_size = 0.2
         self.radius_increment = 0.1
         self.solver_tolerance = 0.001
+        self.consider_acceleration = False
+        self.sphere_z_offset = 0.0
 
 
 class TestPrint:
     def test_should_print(self):
+        printer = FakePrinter()
         test_instruction = Gcode(
-            gcode_path="tests/fixtures/gcode_example.gcode", default_nozzle_speed=40.0
+            gcode_path="tests/fixtures/gcode_example.gcode", default_nozzle_speed=40.0, printer=printer
         )
         test_instruction.read()
 
         simulation_config = FakeSimulationPrint()
-        printer = FakePrinter()
 
         voxel_space = VoxelSpace(
             instruction=test_instruction,