算法优化：动态规划加速，货物运输问题，四边形不等式, 从O(n^2)到O(n^3)

货物运输问题

递归方程为：

更为一般形式的递归方程

看起来是不是像可以使用分治的策略实现，但是min里面子问题太多了，只能使用动态规划的招了。

i,j是什么含义了？动态规划里i,j都是指的是问题规模，对应到货物运输问题指的是什么了？我们从数学上理解i,j是指数组arr = [1,2,3,4,5,6,7,8]两边的标号，或者子问题对应的起始和终止标号，类似与快排里面的数组的标号，实际问题抽象为数组之后，就比较容易理解，在实际问题里面i,j就非常不好理解，i,j是指的是子问题的起始标号和终止标号。

不是动态规划里[i,j]指的是问题规模吗？这里怎么是子问题的起始标号和终止标号，起始标号和终止标号实际上对应的就是问题的规模，j-i就是问题规模，也就是使用了二维表示了一维的规模，快排里也可以不适用quick_sort(arr,p,q),直接使用quick_sort(left_or_right_arr)内部直接使用的left_or_right_arr的0和end。要想原地操作就需要借助于p,q。

假如一个实际问题已经抽象成了一个数组问题，就可以直接从数组的思路考虑，暂时忘记实际的问题，从建立的模型考虑，这样就可以回到熟悉的套路上来。

现在考虑如何实现：

m[i,j] = min{m[i,k-1] + m[k,j] } + w[i,j] i<k<=j；i ==j ,m[i,j] =0

代码实现如下，这种动态规划，根据约束调价i<j,他扫描是沿着对角线方向扫描，对角线从左往右平移，直到扫描到右上角最后一个点，就是最终的解，时间复杂度为O(n^3)

# -*- coding: utf-8 -*-
import numpy as np
# 任意点到第一个点的和
def sum_(arr,n):    
    b = [0]*(n+1)    
    for i in range(1,n+1):
        b[i] = b[i-1] + arr[i-1]        
    return b
# 任意两点之间的数组和
def weight_i_j(b,i,j):
    return b[j+1]-b[i]
    
#  m[i,j] = min{m[i,k-1] + m[k,j] } + w[i,j]  i<k<=j,0<i<j<n+1
# i ==j ,m[i,j] =0

    
def DynamicProgramming(arr):
    n = len(arr)
    b = sum_(arr,n)
    print(b)
    
    # 初始化，这个问题初始化为i==j时为0，也就是对角线上均为0
    m = np.zeros((n,n),int)
    solution = np.full((n,n),-1)
    
    # 扫描的时候注意了，这种类型的问题，扫描顺序不是从左到右，从上到下
    # 扫描的基准是j-i的间距，沿着对角线往右上扫描，这里的i就是为j-i的长度,扫描的最后就是
    # m[0,n-1],这个就是最后的解
    # i <j 所以间距为从1 到 n-1
    for r in range(1,n):
        # 沿着对角线往下走
        for i in range(n-r):
            j = i + r        
            # 首先取第一个k = i+1
            m[i][j] = m[i][i] + m[i+1][j]
            # 用于解的追踪
            solution[i][j] = i
            
            # 比较求最小的
            for k in range(i+2,j+1):
                temp = m[i][k-1] + m[k][j]
                
                if temp < m[i][j]:
                    m[i][j] = temp
                    solution[i][j] = k-1
            # min{m[i,k-1] + m[k,j] } + w[i,j]       
            m[i][j] += weight_i_j(b,i,j)
    
    return m,solution
                
   
arr = [1000,22,13,4,5000,86,57,18]
print(arr)
m,s = DynamicProgramming(arr)
print(m)
print(end='\n')
print(s)

runfile('D:/share/test/dp.py', wdir='D:/share/test')
[1000, 22, 13, 4, 5000, 86, 57, 18]
[0, 1000, 1022, 1035, 1039, 6039, 6125, 6182, 6200]
[[    0  1022  1070  1095  7134 12306 12563 12692]
 [    0     0    35    56  5095 10220 10420 10531]
 [    0     0     0    17  5034 10137 10337 10448]
 [    0     0     0     0  5004 10094 10294 10405]
 [    0     0     0     0     0  5086  5286  5397]
 [    0     0     0     0     0     0   143   236]
 [    0     0     0     0     0     0     0    75]
 [    0     0     0     0     0     0     0     0]]

[[-1  0  0  0  3  3  3  3]
 [-1 -1  1  1  3  4  4  4]
 [-1 -1 -1  2  3  4  4  4]
 [-1 -1 -1 -1  3  4  4  4]
 [-1 -1 -1 -1 -1  4  4  4]
 [-1 -1 -1 -1 -1 -1  5  5]
 [-1 -1 -1 -1 -1 -1 -1  6]
 [-1 -1 -1 -1 -1 -1 -1 -1]]

四边形不等式

结论

回到实际的问题

代码实现如下,实现的时候，只要把原来的i<k<=j,修改成s[i,j-1]<=k<=s[i+1,j]就可以了,时间复杂度为O(n^2)，加速算法有点问题，s初始化都为0，那k不是只能取0吗？程序还怎么走？

def DynamicProgrammingSpeed(arr):
    n = len(arr)
    b = sum_(arr,n)
    print(b)
    
    # 初始化，这个问题初始化为i==j时为0，也就是对角线上均为0
    m = np.zeros((n,n),int)
    solution = np.full((n,n),0)

    
    # 扫描的时候注意了，这种类型的问题，扫描顺序不是从左到右，从上到下
    # 扫描的基准是j-i的间距，沿着对角线往右上扫描，这里的i就是为j-i的长度,扫描的最后就是
    # m[0,n-1],这个就是最后的解
    # i <j 所以间距为从1 到 n-1
    for r in range(1,n):
        # 沿着对角线往下走
        for i in range(n-r):
            j = i + r        
            
            # 利用四边形不等式，把范围缩小为s[i][j-1]到s[i+1][j]
            i1 = solution[i][j-1]
            j1 = solution[i+1][j]
            # 首先取第一个k = i1
            m[i][j] = m[i][i1-1] + m[i1][j]
            # 用于解的追踪
            solution[i][j] = i1-1
            
            # 比较求最小的
            for k in range(i1+1,j1+1):
                temp = m[i][k-1] + m[k][j]
                
                if temp < m[i][j]:
                    m[i][j] = temp
                    solution[i][j] = k-1
            # min{m[i,k-1] + m[k,j] } + w[i,j]       
            m[i][j] += weight_i_j(b,i,j)
    
    return m,solution