装箱问题：将长方体装入具有约束的更大长方体
Packing problem : Fit cuboids into a larger cuboid with constraints
我有一个代表掩码的 3D numpy 数组。这个掩码的“1”元素是没有进行计算的区域（下图中的蓝色区域）。掩码的“0”元素是要进行计算的区域。
每个计算都必须在最小边的长方体块上实现 Nmin。所以我试图找到长方体的最佳分布来填充这个掩码的“0”区域。
下面是掩码示例：mask.npy
及其 3D 表示：
mask = np.load('mask.npy')
fig = plt.figure()
ax1 = fig.add_subplot(projection='3d')
ax1.voxels(mask)
plt.show()
是一个packing problem。目标是将各种长方体打包成一个更大的具有一些约束的长方体（蓝色区域）。有很多可能的组合。
我认为解决这个问题的方法是先确定最大的“0”区域，然后重复这个过程，直到完全填满掩码。
我已经在 2D 中做过这种事情（参见 fdgrid），但我被困在这个 3D 问题上。如何找到“0”的最大块（长方体）？
编辑 1：
我有一个草稿解决方案，对于这个 50x50x50 网格来说有点残酷而且非常慢。
import time
def cuboid_list(mask, Nmin=3):
""" List of all cuboids than could fit in the domain.
List is formated as (volume, (Lx, Ly, Lz)). """
size = mask[mask==0].size
cuboids = [( i*j*k, (i, j, k)) for i, j, k in itertools.product(range(Nmin, nx), range(Nmin, ny), range(Nmin, nz)) if i*j*k <= size]
cuboids.sort(reverse=True)
return cuboids
def search_biggest(mask, cuboids):
"""Search biggest cuboid. """
t0 = time.perf_counter()
tested = 0
origin, size = None, None
for _, (Lx, Ly, Lz) in cuboids:
for o in itertools.product(range(nx-Lx), range(ny-Ly), range(nz-Lz)):
tested += 1
if (mask[o[0]:o[0]+Lx, o[1]:o[1]+Ly, o[2]:o[2]+Lz] == 0).all():
origin, size = o, (Lx, Ly, Lz)
break
if origin:
print(f'{tested} configurations tested in {time.perf_counter() - t0:.3f} s.')
break
return origin, size
# Load mask
mask = np.load('mask.npy')
# List of all cuboids
cuboids = cuboid_list(mask)
# Search biggest
sub = search_biggest(mask, cuboids)
在这种特殊情况下，这导致在 164.984 秒内测试了 2795610 个配置！最大的长方体有origin=(16, 16, 0)和size=(33, 33, 49).
为了使算法适用于更大的几何图形（例如 1024x1024x512 网格），我需要加快这个过程。
编辑 2：
def search_biggest(mask, cuboids):
"""Search biggest cuboid. """
# List of nonzero indexes
forbid = set([tuple(i) for i in np.argwhere(mask != 0)])
t0 = time.perf_counter()
tested = 0
origin, size = None, None
for _, (Lx, Ly, Lz) in cuboids:
# remove cuboid with origin in another cuboid and iterate
for o in set(itertools.product(range(nx-Lx), range(ny-Ly), range(nz-Lz))).difference(forbid):
tested += 1
if (mask[o[0]:o[0]+Lx, o[1]:o[1]+Ly, o[2]:o[2]+Lz] == 0).all():
origin, size = o, (Lx, Ly, Lz)
break
if origin:
print(f'{tested} configurations tested in {time.perf_counter() - t0:.3f} s.')
break
return origin, size
减少可能性导致在 ~2.6 秒内测试 32806 个配置。
有什么想法吗？
编辑 3：
def search_biggest(mask, cuboids):
"""Search biggest cuboid. """
# List of nonzero indexes
forbid = set([tuple(i) for i in np.argwhere(mask != 0)])
t0 = time.perf_counter()
tested = 0
origin, size = None, None
for _, (Lx, Ly, Lz) in cuboids:
tested += 1
tmp = [i for i in itertools.product(range(nx-Lx), range(ny-Ly), range(nz-Lz)) if i not in forbid]
tmp = [i for i in tmp if
all(x not in forbid for x in itertools.product(range(i[0], i[0]+Lx), range(i[1], i[1]+Ly), range(i[2], i[2]+Lz)))]
if tmp:
origin, size = tmp[0], (Lx, Ly, Lz)
print(f'{tested} configurations tested in {time.perf_counter() - t0:.3f} s.')
break
return origin, size
简单胜于复杂！ search_biggest 的替代实施导致在 1.017 秒内测试了 5991 个配置。 !
也许，我可以减少长方体列表来加快速度，但我认为这不是可行的方法。我尝试使用 500x500x500 网格并在一小时后停止脚本，但没有结果...
我还有另外两个可能的公式：
用大小为 Nmin x Nmin x Nmin 的基本长方体填充长方体，然后合并具有公共面的长方体。问题是在某些情况下，不可能仅用这种大小的长方体填充空白区域。
重点搜索平面 (x, y, z=0) 上的空矩形，然后沿 z 扩展它们，直到遇到约束。对 z $\in$ [1, nz-1] 重复此操作，然后合并可以合并的长方体。问题是我最终会得到一些比基本长方体还小的空白区域。
一直在寻找更好的方法来解决这个问题...
鉴于约束是长方体，我们可以通过跳过没有任何变化的 rows/columns/layers 来简化可能性检查。例如，二维蒙版可能看起来像这样，冗余 rows/columns 标记为：
v v
...11.
>...11.
1.....
>1.....
......
可以减少到
..1.
1...
....
我们只需要跟踪每个减少的单元格的真实长度，并用它来计算真实体积。
然而，这意味着我们需要采用不同的方法来寻找最大的长方体；我们不能（轻松地）pre-generate 全部并按音量顺序尝试，因为真实大小因我们测试的点而异。
下面的代码只是尝试所有可能的长方体并跟踪最大的长方体，而且看起来速度很快。
原始代码中还有一些off-by-one错误，我已避免或修复：range(Nmin, nx) -> range(Nmin, nx+1) in cuboid_list()，以及range(nx-Lx) -> range(nx-Lx+1) 在 search_biggest()。最大的长方体居然有size=(34, 34, 50).
import itertools
import time
import numpy as np
class ReducedSolver:
def __init__(self, mask):
self.splits = tuple(
np.unique(
# find the edges of the constraint regions in each axis
[0, mask.shape[ax], *np.where(np.diff(mask, axis=ax, prepend=0))[ax]]
)
for ax in range(3)
)
# extract exactly one value (the lowest corner) for each reduced region
self.mask = np.zeros([len(split) - 1 for split in self.splits], dtype=mask.dtype)
for i, x in enumerate(self.splits[0][:-1]):
for j, y in enumerate(self.splits[1][:-1]):
for k, z in enumerate(self.splits[2][:-1]):
self.mask[i, j, k] = mask[x, y, z]
# less readable:
# self.mask = mask[np.ix_(self.splits[0][:-1], self.splits[1][:-1], self.splits[2][:-1])]
# true sizes of each region
self.sizes = [np.diff(split) for split in self.splits]
def solve(self):
"""Return list of cuboids in the format (origin, size), using a greedy approach."""
nx, ny, nz = self.mask.shape
mask = self.mask.copy()
result = []
while np.any(mask == 0):
t0 = time.perf_counter()
tested = 0
max_volume = 0
# first corner of the cuboid
for x1, y1, z1 in itertools.product(range(nx), range(ny), range(nz)):
if mask[x1, y1, z1]:
continue
# opposite corner of the cuboid
for x2, y2, z2 in itertools.product(range(x1, nx), range(y1, ny), range(z1, nz)):
tested += 1
# slices are exclusive on the the end point
slc = (slice(x1, x2 + 1), slice(y1, y2 + 1), slice(z1, z2 + 1))
# np.any doesn't short-circuit
if any(np.nditer(mask[slc])):
continue
true_size = tuple(np.sum(size[s]) for size, s in zip(self.sizes, slc))
volume = np.prod(true_size)
if volume > max_volume:
max_volume = volume
origin = (
self.splits[0][x1],
self.splits[1][y1],
self.splits[2][z1],
)
# keep track of the region in the reduced mask with `slc`
biggest = ((origin, true_size), slc)
print(f"{tested} configurations tested in {(time.perf_counter() - t0)*1000:.2f} ms.")
result.append(biggest[0])
# mark this cuboid as filled
mask[biggest[1]] = 1
return result
# Load mask
mask = np.load("mask.npy")
# Solve on the simplified grid
reduced = ReducedSolver(mask)
sol = reduced.solve()
在我的机器上，使用 search_biggest() 在示例掩码上找到最大的长方体需要大约 1 秒，但使用我的代码需要大约 4 毫秒，而整个求解 16 个长方体需要大约 20 毫秒。
我有一个代表掩码的 3D numpy 数组。这个掩码的“1”元素是没有进行计算的区域（下图中的蓝色区域）。掩码的“0”元素是要进行计算的区域。
每个计算都必须在最小边的长方体块上实现 Nmin。所以我试图找到长方体的最佳分布来填充这个掩码的“0”区域。
下面是掩码示例：mask.npy
及其 3D 表示：
mask = np.load('mask.npy')
fig = plt.figure()
ax1 = fig.add_subplot(projection='3d')
ax1.voxels(mask)
plt.show()
是一个packing problem。目标是将各种长方体打包成一个更大的具有一些约束的长方体（蓝色区域）。有很多可能的组合。
我认为解决这个问题的方法是先确定最大的“0”区域，然后重复这个过程，直到完全填满掩码。
我已经在 2D 中做过这种事情（参见 fdgrid），但我被困在这个 3D 问题上。如何找到“0”的最大块（长方体）？
编辑 1：
我有一个草稿解决方案，对于这个 50x50x50 网格来说有点残酷而且非常慢。
import time
def cuboid_list(mask, Nmin=3):
""" List of all cuboids than could fit in the domain.
List is formated as (volume, (Lx, Ly, Lz)). """
size = mask[mask==0].size
cuboids = [( i*j*k, (i, j, k)) for i, j, k in itertools.product(range(Nmin, nx), range(Nmin, ny), range(Nmin, nz)) if i*j*k <= size]
cuboids.sort(reverse=True)
return cuboids
def search_biggest(mask, cuboids):
"""Search biggest cuboid. """
t0 = time.perf_counter()
tested = 0
origin, size = None, None
for _, (Lx, Ly, Lz) in cuboids:
for o in itertools.product(range(nx-Lx), range(ny-Ly), range(nz-Lz)):
tested += 1
if (mask[o[0]:o[0]+Lx, o[1]:o[1]+Ly, o[2]:o[2]+Lz] == 0).all():
origin, size = o, (Lx, Ly, Lz)
break
if origin:
print(f'{tested} configurations tested in {time.perf_counter() - t0:.3f} s.')
break
return origin, size
# Load mask
mask = np.load('mask.npy')
# List of all cuboids
cuboids = cuboid_list(mask)
# Search biggest
sub = search_biggest(mask, cuboids)
在这种特殊情况下，这导致在 164.984 秒内测试了 2795610 个配置！最大的长方体有origin=(16, 16, 0)和size=(33, 33, 49).
为了使算法适用于更大的几何图形（例如 1024x1024x512 网格），我需要加快这个过程。
编辑 2：
def search_biggest(mask, cuboids):
"""Search biggest cuboid. """
# List of nonzero indexes
forbid = set([tuple(i) for i in np.argwhere(mask != 0)])
t0 = time.perf_counter()
tested = 0
origin, size = None, None
for _, (Lx, Ly, Lz) in cuboids:
# remove cuboid with origin in another cuboid and iterate
for o in set(itertools.product(range(nx-Lx), range(ny-Ly), range(nz-Lz))).difference(forbid):
tested += 1
if (mask[o[0]:o[0]+Lx, o[1]:o[1]+Ly, o[2]:o[2]+Lz] == 0).all():
origin, size = o, (Lx, Ly, Lz)
break
if origin:
print(f'{tested} configurations tested in {time.perf_counter() - t0:.3f} s.')
break
return origin, size
减少可能性导致在 ~2.6 秒内测试 32806 个配置。
有什么想法吗？
编辑 3：
def search_biggest(mask, cuboids):
"""Search biggest cuboid. """
# List of nonzero indexes
forbid = set([tuple(i) for i in np.argwhere(mask != 0)])
t0 = time.perf_counter()
tested = 0
origin, size = None, None
for _, (Lx, Ly, Lz) in cuboids:
tested += 1
tmp = [i for i in itertools.product(range(nx-Lx), range(ny-Ly), range(nz-Lz)) if i not in forbid]
tmp = [i for i in tmp if
all(x not in forbid for x in itertools.product(range(i[0], i[0]+Lx), range(i[1], i[1]+Ly), range(i[2], i[2]+Lz)))]
if tmp:
origin, size = tmp[0], (Lx, Ly, Lz)
print(f'{tested} configurations tested in {time.perf_counter() - t0:.3f} s.')
break
return origin, size
简单胜于复杂！ search_biggest 的替代实施导致在 1.017 秒内测试了 5991 个配置。 !
也许，我可以减少长方体列表来加快速度，但我认为这不是可行的方法。我尝试使用 500x500x500 网格并在一小时后停止脚本，但没有结果...
我还有另外两个可能的公式：
用大小为
Nmin x Nmin x Nmin的基本长方体填充长方体，然后合并具有公共面的长方体。问题是在某些情况下，不可能仅用这种大小的长方体填充空白区域。重点搜索平面 (x, y, z=0) 上的空矩形，然后沿 z 扩展它们，直到遇到约束。对 z $\in$ [1, nz-1] 重复此操作，然后合并可以合并的长方体。问题是我最终会得到一些比基本长方体还小的空白区域。
一直在寻找更好的方法来解决这个问题...
鉴于约束是长方体，我们可以通过跳过没有任何变化的 rows/columns/layers 来简化可能性检查。例如，二维蒙版可能看起来像这样，冗余 rows/columns 标记为：
v v
...11.
>...11.
1.....
>1.....
......
可以减少到
..1.
1...
....
我们只需要跟踪每个减少的单元格的真实长度，并用它来计算真实体积。 然而，这意味着我们需要采用不同的方法来寻找最大的长方体；我们不能（轻松地）pre-generate 全部并按音量顺序尝试，因为真实大小因我们测试的点而异。 下面的代码只是尝试所有可能的长方体并跟踪最大的长方体，而且看起来速度很快。
原始代码中还有一些off-by-one错误，我已避免或修复：range(Nmin, nx) -> range(Nmin, nx+1) in cuboid_list()，以及range(nx-Lx) -> range(nx-Lx+1) 在 search_biggest()。最大的长方体居然有size=(34, 34, 50).
import itertools
import time
import numpy as np
class ReducedSolver:
def __init__(self, mask):
self.splits = tuple(
np.unique(
# find the edges of the constraint regions in each axis
[0, mask.shape[ax], *np.where(np.diff(mask, axis=ax, prepend=0))[ax]]
)
for ax in range(3)
)
# extract exactly one value (the lowest corner) for each reduced region
self.mask = np.zeros([len(split) - 1 for split in self.splits], dtype=mask.dtype)
for i, x in enumerate(self.splits[0][:-1]):
for j, y in enumerate(self.splits[1][:-1]):
for k, z in enumerate(self.splits[2][:-1]):
self.mask[i, j, k] = mask[x, y, z]
# less readable:
# self.mask = mask[np.ix_(self.splits[0][:-1], self.splits[1][:-1], self.splits[2][:-1])]
# true sizes of each region
self.sizes = [np.diff(split) for split in self.splits]
def solve(self):
"""Return list of cuboids in the format (origin, size), using a greedy approach."""
nx, ny, nz = self.mask.shape
mask = self.mask.copy()
result = []
while np.any(mask == 0):
t0 = time.perf_counter()
tested = 0
max_volume = 0
# first corner of the cuboid
for x1, y1, z1 in itertools.product(range(nx), range(ny), range(nz)):
if mask[x1, y1, z1]:
continue
# opposite corner of the cuboid
for x2, y2, z2 in itertools.product(range(x1, nx), range(y1, ny), range(z1, nz)):
tested += 1
# slices are exclusive on the the end point
slc = (slice(x1, x2 + 1), slice(y1, y2 + 1), slice(z1, z2 + 1))
# np.any doesn't short-circuit
if any(np.nditer(mask[slc])):
continue
true_size = tuple(np.sum(size[s]) for size, s in zip(self.sizes, slc))
volume = np.prod(true_size)
if volume > max_volume:
max_volume = volume
origin = (
self.splits[0][x1],
self.splits[1][y1],
self.splits[2][z1],
)
# keep track of the region in the reduced mask with `slc`
biggest = ((origin, true_size), slc)
print(f"{tested} configurations tested in {(time.perf_counter() - t0)*1000:.2f} ms.")
result.append(biggest[0])
# mark this cuboid as filled
mask[biggest[1]] = 1
return result
# Load mask
mask = np.load("mask.npy")
# Solve on the simplified grid
reduced = ReducedSolver(mask)
sol = reduced.solve()
在我的机器上，使用 search_biggest() 在示例掩码上找到最大的长方体需要大约 1 秒，但使用我的代码需要大约 4 毫秒，而整个求解 16 个长方体需要大约 20 毫秒。