北京大学Julia语言入门讲义第25章： Julia编程示例–动态规划

using Random
Random.seed!(101)
M = rand(0:9, 5, 6)

function matdp_exaust(M; DEBUG=true)
    m, n = size(M)
    paths = []

    # 从第i行开始的递归方法
    # 直接修改外层的paths
    # 输入为截止到第i行的部分路径path_in和i
    function path_remain(path_in, i)
        if(i==m) # 最后一行已解决
            push!(paths, copy(path_in))
            return
        else # i < m, 考虑下一行的三种不同情况
            # 正下方
            path = copy(path_in)
            path[(i+1):end] .= 0
            path[i+1] = path[i] # 正下方
            path_remain(path, i+1)
            # 左下方
            if path[i] > 1
                path = copy(path_in)
                path[(i+1):end] .= 0
                path[i+1] = path[i] - 1
                path_remain(path, i+1)
            end
            # 右下方
            if path[i] < n
                path = copy(path_in)
                path[(i+1):end] .= 0
                path[i+1] = path[i] + 1
                path_remain(path, i+1)
            end
        end # if i==n 
        # 函数不需要显式地返回值，通过外层变量paths记录遍历结果
    end # function path_remain 

    # 只要遍历第一行
    path = fill(0, m)
    for start=1:n 
        path[1] = start 
        path_remain(path, 1)
    end

    #println(paths)
    npaths = length(paths)
    sums = map(sum, [[M[i, p[i]] for i=1:m] for p in paths])
    id = argmin(sums)
    minen = sums[id]
    if(DEBUG)
        println("矩阵大小：", (m, n), " 路径数：", npaths)
    end

    return paths[id], minen 
end

function test_matdp(m=8, n=m, solver=matdp_exaust; DEBUG=true)
    Random.seed!(101)
    M = rand(0:9, m, n)
    bestpath, minen = solver(M; DEBUG=DEBUG)
    if(DEBUG)
    end
    if(DEBUG)
        println("矩阵大小：", (m, n))
        println("最小能量和：", minen)
        println("最小能量路径：")
        println(bestpath)
        M
    end
end
test_matdp(5, 6, matdp_exaust, DEBUG=true)

## 仅考虑一个最优解的递归方法
function matdp_recurse(M; DEBUG=true)
    m, n = size(M)

    # 递归函数。
    # 返回从(i,j)元素向下， 找到的最小能量， 和下一行应去的列号。
    function next_ij(M, i, j)
        m, n = size(M)

        ## 初始化情形
        if i == m
            return M[m, j], n+1
        end

        ## 递归求解情形
        # 考虑第i+1行的j-1, j, j+1对应的最小能量
        js = unique(max.(1, min.(n, (j-1):(j+1))))
        min_ens = [next_ij(M, i+1, jj)[1] for jj in js]
        ind = argmin(min_ens)
        jnext = js[ind]
        min_en = min_ens[ind] + M[i, j]
        return min_en, jnext
    end

    # 只要从第一行开始找到最优开始位置，然后依次下行就可以
    minens = [next_ij(M, 1, j)[1] for j=1:n]
    j1 = argmin(minens)
    minen = minens[j1] # 全局最小能量
    bestpath = [j1;]
    for i = 2:m 
        en, j = next_ij(M, i-1, bestpath[i-1])
        push!(bestpath, j)
    end
    
    return bestpath, minen 
end

## 用动态规划方法求解的函数
function matdp_dp(M; DEBUG=true)
    m, n = size(M)

    # 对M的每个元素，求出并记住最优路径
    MM = Array{Any}(undef, (m, n))
    # MM[i,j] 存储一条最优子路径和最优子路径的部分和

    # 最后一行，最优路径即本身
    for j = 1:n
        MM[m, j] = ([j,], M[m,j])
    end

    # 第(m-1)到第1行：
    for i=(m-1):-1:1
        for j = 1:n 
            # 考虑M[i,j]的最优路径
            # 只要比较MM[i+1, j-1], MM[i+1, j], MM[i+1, j+1]
            js = [max(1, j-1), j, min(n, j+1)]
            # js: 下一行的三个可能列号
            vs = [MM[i+1, jj][2] for jj in js]
            # vs: 下一行的三个可能路径的部分和
            jmin = js[argmin(vs)]
            # jmin: 下一行的最优列号
            vmin = MM[i+1, jmin][2]
            # vmin: 下一行的最优值
            MM[i, j] = ([j; MM[i+1,jmin][1]], M[i,j] + vmin)
        end # for j
    end # for i 

    # 比较第一行的n个最优值，找到全局最优路径
    best1 = [MM[1,j][2] for j=1:n]
    bestj = argmin(best1)
    bestpath = MM[1, bestj][1]
    besten = MM[1, bestj][2]

    return bestpath, besten
end

using BenchmarkTools
@benchmark test_matdp(10, 8, matdp_exaust, DEBUG=false)
@benchmark test_matdp(10, 8, matdp_recurse, DEBUG=false)
@benchmark test_matdp(10, 8, matdp_dp, DEBUG=false)

# MIT公开课用迭代器进行穷举的方法
struct Paths # 用来表示所有路径的类型
    m::Int # 矩阵行数
    n::Int # 矩阵列数
end

# 要自定义迭代器，必须定义Base.iterate(itr)和Base.iterate(itr, state)
# 用来初始化的函数，从矩阵最左边一列开始
Base.iterate(p::Paths) = fill(1, p.m), fill(1, p.m) 

# 用来定义迭代器长度的类型
Base.IteratorSize(::Type{Paths}) = SizeUnknown()

# 用来迭代进入下一个状态的函数
function Base.iterate(p::Paths, state)
    if state ≠ fill(p.n, p.m) # 当路径为最后一列时结束
        newstate = next(state, p.n)
        return newstate, newstate
    end
end

# 给定一条路径，求下一条路径，n为列数
function next(path, n)
    m = length(path) # 行数

    k = m # 从最后一行开始向上查找
    # 从现有路径的最后一行向上查找，直到找到一个可以向右移动一格的位置
    # 下面的while结构是如果第k行还没有到第一行，
    # 并且当前路径第k行已经在最后一列，
    # 或者当前路径第k行如果向右一列会与第k-1行相比过于靠右无法连接，
    # 则第k行不变，考虑上一行即第k-1行
    while  k≥2 && ( path[k]==n || path[k]+1 > path[k-1]+1 )
        k -= 1
    end 
    
    # 上面的循环找到了适当的行号k，
    # 此行的当前路径点可以右移且不会与第k-1行割裂，
    # 于是右移第k行的这个点，第1:(k-1)行保持不变
    path[k] +=1 

    # 前面用倒序从最后一行开始找到了可以右移且不会与前面行割裂的一个行号k，
    # 将当前路径第k行右移一个格，
    # 然后找到从第k行开始的最左边的一条向下路径
    for j = (k+1):m
        # 每一行依次向左下方
        path[j] = max(path[j-1]-1,1)
    end

    return(path)
end


# 用迭代器获得所有路径
# m: 矩阵行数
# n: 矩阵列数
function allpaths(m,n)
    v = Vector{Int}[] # v的每个元素是一个整数向量，为一条路径
    paths = Paths(m,n) # 一个迭代器
    for p ∈ paths
        push!(v, copy(p))
    end
    return v
end

# 返回最优路径和最优值的函数
function matdp_iter(M; DEBUG=true)
    (m, n) = size(M)
    paths = allpaths(m, n)

    npaths = length(paths)
    sums = map(sum, [[M[i, p[i]] for i=1:m] for p in paths])
    id = argmin(sums)
    minen = sums[id]
    if(DEBUG)
        println("矩阵大小：", n, " 路径数：", npaths)
    end

    return paths[id], minen

end # function matdp_iter