Numoku solver: Solves i.e https://twitter.com/1to9puzzle/status/1255848687298740224

numoku.jl

using JuMP, ConstraintSolver, Combinatorics
const CS = ConstraintSolver

"""
    get_possible_splits(p, n, wished_sum, split_point)

p = possible numbers like 1:9 in a normal numoku
n = number of elements in a row
wished_sum = the sum of a row
split_point = index of where the left part should be a multiple of the right part or the other way around
`left_sum = sum(row[1:split_point])`
`right_sum = sum(row[split_point+1:end])`
Formally:
Return an Array of size (x,n) where each `row` fulfills `sum(row) == 30` and `left_sum = N*right_sum` or `N*left_sum = right_sum` where N is discrete
"""
function get_possible_splits(p, n, wished_sum, split_point)
    all_possible = collect(permutations(p, n))
    correct_sum = filter(v->sum(v) == wished_sum, all_possible)
    vecs = filter(v->sum(v[1:split_point]) % sum(v[split_point+1:end]) == 0 || sum(v[split_point+1:end]) % sum(v[1:split_point]) == 0, correct_sum)
    a = zeros(Int, (length(vecs), n))
    for j=1:length(vecs)
        a[j,:] .= vecs[j]
    end
    return a
end

function solve_numoku(grid, splits, block_height, block_width, wished_sum)
    height, width = size(grid)
    m = Model(CS.Optimizer)
    @variable(m, 1 <= x[1:height, 1:width] <= 9, Int)
    for i=1:height, j=1:width
        if grid[i,j] != 0
            @constraint(m, x[i,j] == grid[i,j])
        end
    end
    
    nblocks_height = convert(Int, height/block_height)
    nblocks_width = convert(Int, width/block_width)
    for i=1:nblocks_height, j=1:nblocks_width
        vec = [x[(i-1)*block_height+i1, (j-1)*block_height+j1] for i1 = 1:block_height, j1=1:block_width]
        @constraint(m, vec[:] in CS.AllDifferentSet())
    end
    for i=1:height
        @constraint(m, x[i, :] in CS.AllDifferentSet())
    end
    for j=1:width
        @constraint(m, x[:, j] in CS.AllDifferentSet())
    end

    split_tables_row = Vector{Array{Int, 2}}(undef, width-1)
    split_tables_col = Vector{Array{Int, 2}}(undef, height-1)
    
    for i=1:width-1
        split_tables_row[i] = get_possible_splits(1:9, width, wished_sum, i)
    end
    for i=1:height-1
        split_tables_col[i] = get_possible_splits(1:9, height, wished_sum, i)
    end

    
    for i=1:length(splits)
        split = splits[i]
        if split[3] == :r
            py = split[1]
            px = split[2]
            @constraint(m, x[py,:] in CS.TableSet(split_tables_row[px]))
        elseif split[3] == :d
            py = split[1]
            px = split[2]
            @constraint(m, x[:,px] in CS.TableSet(split_tables_col[py]))
        else 
            println("Only :r and :d splits are supported...")
        end
    end

    optimize!(m)
    com = JuMP.backend(m).optimizer.model.inner
    println(JuMP.termination_status(m))
    @show convert.(Int, JuMP.value.(x))
end

grid = zeros(Int, (6,6))
grid[2,2] = 5
grid[2,5] = 1
grid[5,2] = 8
grid[5,5] = 2

splits = Vector{Tuple{Int,Int,Symbol}}()
push!(splits, (1,1,:r))
push!(splits, (1,1,:d))
push!(splits, (1,2,:r))
push!(splits, (1,3,:d))
push!(splits, (1,4,:d))
push!(splits, (2,3,:r))
push!(splits, (2,4,:d))
push!(splits, (2,5,:d))
push!(splits, (2,5,:r))
push!(splits, (3,4,:d))
push!(splits, (4,1,:r))
push!(splits, (4,2,:r))
push!(splits, (4,3,:r))
push!(splits, (4,5,:r))
push!(splits, (4,5,:d))
push!(splits, (5,1,:d))
push!(splits, (5,1,:r))
push!(splits, (5,3,:d))
push!(splits, (5,4,:r))
push!(splits, (5,5,:r))
push!(splits, (5,6,:d))
push!(splits, (6,1,:r))
push!(splits, (6,5,:r))

@time solve_numoku(grid, splits, 3, 3, 30)