https://www.youtube.com/watch?v=DmT80OseAGs

sudoku.py

#!/usr/bin/env python3

from functools import reduce
import numpy as np

def box(n):
    return np.repeat(
        np.repeat(
            np.arange(1, n**2+1).reshape(n, n),
            n, 1
        ),
        n, 0
    )

def indices(n):
    nsq = n * n
    idx = np.arange(nsq)
    return np.array(np.meshgrid(idx, idx)).T.reshape(nsq, nsq, 2)

def rcb(n):
    return np.vstack([indices(n).T, box(n).reshape(n**2, n**2, 1).T]).T

def cmap(n):
    rcbs = rcb(n)
    return np.diagonal(np.equal.outer(rcbs, rcbs), axis1=2, axis2=5).max(axis=4).astype(int)

def avl(cont, puz, r, c):
    return np.delete(np.arange(len(puz) + 1), cont[r, c] * puz)

def at(puz, val, r, c):
    puz = puz.copy()
    puz[r, c] = val
    return puz

def pvec(cont, puz, r, c):
    # need explicit reshape and type for empty arrays
    return np.array([
        at(puz, x, r, c) for x in avl(cont, puz, r, c)],
        dtype=int,
    ).reshape(
        -1, *puz.shape
    )

def pvex(r, c, cont, puzs):
    return np.vstack([pvec(cont, puz, r, c) for puz in puzs])

def emt(puz):
    return np.argwhere(puz == 0)

def svec(puz):
    cont = cmap(round(len(puz) ** 0.5))
    return reduce(lambda puzs, idx: pvex(*idx, cont, puzs), emt(puz), [puz])


if __name__ == '__main__':
    # Let's try it!
    s44 = np.array([
        [0, 0, 0, 0],
        [0, 0, 2, 1],
        [3, 0, 0, 4],
        [0, 0, 0, 0],
    ])
    print(svec(s44))
    
    s99 = np.array([
        [0, 0, 1, 6, 9, 0, 5, 0, 0],
        [4, 0, 0, 2, 7, 0, 0, 0, 1],
        [0, 7, 0, 0, 0, 0, 0, 9, 0],
        [0, 0, 0, 0, 0, 0, 0, 3, 0],
        [0, 0, 0, 4, 3, 0, 0, 0, 7],
        [0, 0, 0, 7, 8, 0, 6, 0, 0],
        [0, 0, 6, 0, 0, 0, 8, 0, 5],
        [0, 2, 0, 1, 4, 0, 0, 6, 0],
        [0, 1, 0, 3, 5, 0, 0, 4, 0],
    ])
    print(svec(s99))