魔術磁力方塊是之前買給小朋友玩的玩具，裡面有個組好的 3x3 積木，是用七種不同顏色立體方塊組合成的，可以利用磁性吸來吸去，做出不同的造型。除此之外還有發一副牌，上面有 108 種只有外型沒有顏色的題目，難度分成兩顆至四顆星，可以挑戰看看怎麼利用指定的方塊堆成牌面上的形狀。
四顆星的題目雖然只有八張，但有些真的有點難啊！
其中有一題編號 103 的題目，小朋友組不出來跑來問這個要怎麼辦？這個形狀並沒有顯著的特徵，也就是沒有邊邊角角是一定要某種特定才排得進去的，就沒有辦法減少搜尋嘗試的次數。
嘗試了一下覺得不容易啊，那為什麼不來寫個程式幫忙找一下呢？於是就有了這個 tetris3d 的計畫。開發的過程有些有趣的地方覺得值得紀錄一下，所以就有了這篇。
但在這之前，想要先解決如何顯示的問題。
原因是因為到時候找到如何組起來之後，想要有一個能夠 表達互動式 3D 的介面，最好是跟觸碰螢幕配合地不錯，這樣才能夠拿個平板給小朋友看，自己滑滑轉轉看能不能造圖組起來。
在寫這個顯示的部分，其實也花了不少時間，托最近 Deep Learning 很熱門的福，其實有不少好用的 libraries 可以試試看，不再只有 matplot 了。挑了幾個過去有用過與沒有用過的試試，最後採用 plotly 這個來畫。
plotly 的 3D charts 有許多範例可以選擇，但似乎不知道要挑哪一種呈現比較好，由於初步的規劃是想用一連串的 3D 的 numpy array 運算來解這些問題，直覺地想要用 3D Scatter，反正只是在中間顯示點，夠表達就好了，不過畫出來放大縮小來看其實不太理想。
再繼續了解之後，最後決定使用 Isosurface 來畫圖，雖然範例乍看一下好像比較無關，但其實可以畫出想要的效果，還能夠設定透明度呢！
import plotly.graph_objects as go
import dash
import dash_core_components as dcc
import dash_html_components as html
import numpy as np
def main():
points = [(0, 0, 0), (1, 0, 0), (2, 2, 2)]
data = []
for x, y, z in points:
points = np.mgrid[x:x+1:2j, y:y+1:2j, z:z+1:2j]
values = np.ones(points[0].shape) * 0.5
data.append(go.Isosurface(
x=points[0].flatten(),
y=points[1].flatten(),
z=points[2].flatten(),
value=values.flatten(),
colorscale=[[0, 'red'], [1, 'red']],
showscale=False,
opacity=0.5))
fig = go.Figure(data=data)
fig.update_layout(margin=dict(l=0, r=0, b=0, t=32))
fig.update_layout(scene_aspectmode='cube')
fig.write_html('isosurface.html')
app = dash.Dash()
app.layout = html.Div([dcc.Graph(figure=fig)])
app.run_server(host='0.0.0.0', debug=True, use_reloader=True)
if __name__ == "__main__":
main()接著就是來想該如何寫程式組合方塊，找到可以符合題目的組合呢？
首先來表達題目跟可以組合的方塊，方塊有七種剛好以彩虹的顏色命名。由於不會有一個絕對的大小，所以用含有方塊的點為 1 的 三維座標組來建立，而題目也是類似的方式，只是選擇用灰色來表示。
# color in read bricks
VALID_GROUPS = {
'red': [(0, 1, 0), (0, 0, 1), (0, 1, 1), (0, 2, 1)],
'orange': [(0, 1, 0), (1, 0, 0), (1, 1, 0), (1, 1, 1)],
'green': [(0, 0, 0), (0, 1, 0), (1, 1, 0), (2, 1, 0)],
'blue': [(0, 1, 1), (1, 0, 0), (1, 1, 0), (1, 1, 1)],
'indigo': [(0, 0, 0), (0, 1, 0), (1, 0, 0), (1, 0, 1)],
'violet': [(0, 0, 0), (0, 1, 0), (0, 1, 1), (0, 2, 1)],
'yellow': [(0, 0, 0), (0, 0, 1), (1, 0, 0)],
}不過有顏色的積木方塊與題目方塊還有一個最大的不同，就是積木方塊可以 旋轉，由於又有三維 x, y, z, 所以可能可以任意轉來轉去。
# https://en.wikipedia.org/wiki/Rotation_matrix
# Rx = [[ 1, 0, 0],
# [ 0, cos(), -sin()],
# [ 0, sin(), cos()]]
# Ry = [[ cos(), 0, sin()],
# [ 0, 1, 0],
# [ -sin(), 0, cos()]]
# Rz = [[ cos(), -sin(), 0],
# [ sin(), cos(), 0],
# [ 0, 0, 1]]
ROT90_X = np.array(
[[ 1, 0, 0],
[ 0, 0, -1],
[ 0, 1, 0]])
ROT90_Y = np.array(
[[ 0, 0, 1],
[ 0, 1, 0],
[-1, 0, 0]])
ROT90_Z = np.array(
[[ 0, -1 , 0],
[ 1, 0, 0],
[ 0, 0, 1]])所以每種顏色的基礎方塊會再延伸出 4 4 4 = 64 個方塊，這邊就只有歸零原點，並沒有另外把因為對稱而重複的表達刪除。
有了這些積木方塊組之後，其實最簡單的想法就是把這些方塊到處移動移動，看看能不能組出題目的形狀，好像沒有太困難，就是把題目形狀的點減去顏色方塊點就可以了。不過座標點不太好直接作算數運算，所以要先座標點轉成某種 shape 大小的 numpy array。
@property
def points(self) -> np.ndarray:
''' returns absolute points with self.offset and self.rotate '''
cur_points = self._points
for ri, rot in enumerate([ROT90_X, ROT90_Y, ROT90_Z]):
for rt in range(self.rotate[ri]):
cur_points = cur_points.dot(rot)
cur_points -= cur_points.min(axis=0)
cur_points += self.offset
return cur_points.astype(np.int8)
def volume(self, shape: Scaler3D=None) -> np.ndarray:
'''
returns volume of points with self.offset and self.rotate
'''
if shape is None:
shape = (self.points.max(axis=0) + 1).tolist()
npy = np.zeros(shape, dtype=np.int8)
for p in self.points.tolist():
npy[p[0], p[1], p[2]] = 1
return npy然後找出全部積木減完後剩下的 numpy array 都為 0 的組合這就是一組答案了。
但實際實作的時候，發覺並沒有想像中簡單，以問題 103. 為例子好了，最大的值為 4，所以會想要把座標點對應到一個 shape 為 (6, 6, 6) 的 numpy array，因為顏色積木最大值是 3，這樣就不用考慮超過的處理。稍微算一下，其實要比對的組合式天文數字，至少一般電腦室用 CPU 是跑不動的。
簡化地算所有 permuation 的個數
(4 x 4 x 4) x (6 x 6 x 6) ^ 7 = 96479729228174488169059713024
顯然不太能用 brute force 的方式，而稍微思考一下之後，會發現在 ^ 7 的第一個方塊運算的時候，就可以濾掉很多不合理的例子，最基礎的就是減完之後有負的值在裡面，這代表方塊是排不進去的情況。
在越初期就能濾掉越多，就越能大幅降低比對的組合
所以那些顏色的排列是有一點順序的，基本上就是格數越多越複雜的方塊先把前面。這其實跟人腦手動解這些題目一樣，會先把難排的先擺進去，然後一但後面都擺不好了，就移動一下旋轉一下。
然而單純使用減完之後有負值的情況，就變得比較實際一點了，大致上是幾天內可以算出一個題目。嘗試了 multiple processes 等等的額外加速，似乎有機會在一個週末可以跑完，於是就耐心地等。
結果是 memory 爆掉了
為此還調整了 data type，有注意到的話上面的 points 是用 np.int8 儲存的，不過對一台幾年前的筆記型電腦來說，16 GiB 還是不夠放。
那再回頭思考看看有沒有辦法再更有效地降低組合，有什麼情況可以更早知道這組排列是不可能有結果的呢？
其實顏色積木有七個，每個至少都是三格四格，而它們有一個共通點，就是這些格子一定會連在一起！
所以只要剩下來的 numpy array 中，包含有如底下 isolated 的情況，就可以跳過了
這邊就沒有再繼續找相鄰有一點連在一起的情況，如果有需要是還可以依此類推的。
要怎麼找出這種 isolated 的點呢？對於某個是 1 的點，上下左右前後都不是 1 就是了，不過這樣直接寫起來太多 loop，其實還是不夠快，並不太有效率。
再想一下，啊這不就是熟悉不過的 Convolution 嗎？
會了方便解釋，用 2D 的圖示來表達：
邊邊角角就一樣是用 padding 來處理，Convolution 完後會得到一張一樣大的 numpy array，裡面就是鄰居的有點的個數。如果有玩過的話，啊這不就是 採地雷 的嗎！
最後只要將這個 array 跟原本的做 logical and，如果裡面有任何值 > 1，那就是有 isolated points 了，可以不用再繼續，趕緊算別的吧。
class Solver():
...
kernel = k = np.array([
[[0, 0, 0], [0, 1, 0], [0, 0, 0]],
[[0, 1, 0], [1, 0, 1], [0, 1, 0]],
[[0, 0, 0], [0, 1, 0], [0, 0, 0]],
], dtype=np.int8)
...
@staticmethod
def fit(sol, volbricks, stack_vol, kernel) -> typing.List[Solution]:
possible_solutions = []
delta_vol = sol.volume - stack_vol
left_mask = delta_vol.min(axis=(1, 2, 3)) >= 0
for vb in itertools.compress(volbricks, left_mask):
left_vol = sol.volume - vb.volume
convs = scipy.ndimage.convolve(
left_vol, kernel, mode='constant', cval=0)
# refactored version with np.logical_and
if np.count_nonzero(np.logical_and(left_vol > 0, convs == 0)) > 0:
continue
possible_solutions.append(Solution(
sol.volume - vb.volume, sol.bricks + [vb.brick]))
# orignal version with for loop
'''
is_possible = True
for p in brick.vol2pts(left_vol):
if convs[p[0], p[1], p[2]] in [0]:
is_possible = False
break
if is_possible:
possible_solutions.append(Solution(
sol.volume - vb.volume, sol.bricks + [vb.brick]))
'''
return possible_solutions其實原本不是使用 logical and 來做的，其實是寫一個 for loop，認為剩下的點數應該不會很多，這樣跑起來大約花 數十秒 就可以跑完了，是最有效的 optimization。不過透過 profile 分析，這個 for loop 其實佔不少時間，所以乾脆採用 logical and 直接一次算完再判斷反而是比較快，大約快上 三倍，所以在當時跑的 i5-8250U 上頂多 10 幾秒 就可以算完。
這樣已經是相當可以使用的了，會算出多組可能的組合排列，這裡就只挑第一組畫出來，並沒有去計算結果是不是有重複總共有幾種這樣。
而在後來買的 Apple M1 MacBook Air 上，用 brew 上的 arm64 numpy 來算，也是快大約 三倍，要是當時候是在這台開發的話，可能就不會再去做 logical and 的 optimization。
最後這就是當初燒腦也解不開的題目 103. 的解答之一
最後就是把其他的題目也都輸入進去試試看，發覺也都是可以的，像是 107. 看起來像是座沙發的樣子。
而當初因為是在自己的 private git repo 開發，所以有另外把 code 整理出來放到 github 上，https://github.com/yumaokao/tetris3d，但不知道什麼時候會來寫個 Readme.md，哈哈哈，就先這樣啦！