# Efficient way to fill up 3D space with cubes of different sizes?

I'm working on a script that will generate num_cubes number of cubes in a pre-defined 3D space. The cube dimensions will be random, using random.choice(cubeLengths), where cubeLengths = [1,2,4] (I will play with weights afterwards to make it more interesting).

I'd like to have it so that the cubes are as close together as possible without overlap, instead of having a lot of space in between them. I'd looked through this answer dealing with overlapping cubes that has a similar concept, and I was able to make that work, but I'm looking to make mine more compact, to achieve a specific look like this:

From the looks of things, it does not seem like his cubes overlap. Currently, I'm able to generate a single row of randomly sized cubes back-to-back, but then once I need to go to a second row, I get stuck. My code is obviously not optimized, nor is it the best, but here's what I have so far:

What could be the right approach to tackling this problem? Am I overthinking this? Or is this specific problem require a lot more math than it seems?

Here's a more modern solution to this problem, using Numpy and the built-in KDTree algorithm in Blender's mathutils. It generates random points for our cubes, then prevents overlaps by measureing the distance between each cube and its closest neighbor.

import bpy
import numpy as np
from mathutils import kdtree, Vector
import math

xmin, ymin, zmin = [ -5, -5, -5 ]
xmax, ymax, zmax = [ 5, 5, 5 ]
n = 500 # Number of cubes / density

x = xmin + np.random.random(n) * (xmax-xmin)
y = ymin + np.random.random(n) * (ymax-ymin)
z = zmin + np.random.random(n) * (zmax-zmin)
cube_centers =  np.array(list(zip(x,y,z)))

# For each point we'll find the nearest neighbor,
# and we'll set the radius of the current point to be the distance from the closest point * gap
gap = 1 / math.sqrt(2) # This value ensures no overlap - if you want a higher density and don't care about overlaps you can increase this value
kd = kdtree.KDTree(n) # A KDtree algorithm helps quickly calculate distances between points

for i, v in enumerate(cube_centers):
kd.insert(Vector(v), i)

kd.balance()

for cc in cube_centers:
co, index, dist = kd.find_n(Vector(cc), 2)[1]