文章目录
前言
本文主要是通过讲动态规划的原理，（最简单的是切钢材），应用于NLP中，常见的是编辑距离，在生物学中应用的是Needleman-Wunsch算法；本文主要是针对学习过程进行记录笔记
一、动态规划原理
1、定义
1） 动态规划（dynamic programming) :将原问题划分为子问题，通过组合子问题的解来求解原问题。不同的子问题具有公共的子子问题（子问题的求解是递归进行的，将其划分为更小的子子问题）。在这种情况下，动态规划算法对每个子子问题只求解一次，将其解保存在一个表格中，从而无需每次求解一个子子问题时都重新计算，避免这种不必要的计算工作。
2） 动态规划方法通常用来求解最优化问题（optimization problem)。这种问题通常可以有很多可行解，每个解都有一个值，我们希望寻找具有最优值（最小值或最大值）的解，我们称这样的解为问题的一个最优解（an optimal solution), 而不是最优解（the optimal solution),因为有很多个解都达到最优值。
3）举例小故事—便于理解
A * "1+1+1+1+1+1+1+1 =？" *
A : "上面等式的值是多少"
B : *计算* "8!"
A *在上面等式的左边写上 "1+" *
A : "此时等式的值为多少"
B : *quickly* "9!"
A : "你怎么这么快就知道答案了"
A : "只要在8的基础上加1就行了"
A : "所以你不用重新计算因为你记住了第一个等式的值为8!动态规划算法也可以说是 '记住求过的解来节省时间'"
2、特点
1）针对机器学习算法而言：不需要大量数据和足够的特征值
2）相对于BFS和DFS、贪心算法而言：
1、BFS和DFS无法搜索过大的空间（时间、空间复杂度过大）
2、贪心算法会带来局部最优解而非全局最优解
3）
3、切钢材
简单解释如下图
3、代码实现
from collections import defaultdict
import time
original_prica = [1, 5, 8, 9, 10, 17, 17, 20, 24, 30, 35]
price_dic=defaultdict(int)
for i,p in enumerate(original_prica):
price_dic[i+1]=p
print(price_dic)
def get_pr(n):
return max([price_dic[n]]+[get_pr(int(n-i))+get_pr(int(i)) for i in range(1,(n+1)//2)])
"""
n=3:递归过程，以3为举例
max([price_dic[3]]+[get_pr(2)+get_pr(1)])===>max(8,6)=8
get_pr[2]==>return max([price_dic[2]]+[get_pr(1)+get_pr(1)])===>return max(5,2)=5
get_pr(1)==>max(price_dic(1)=1
"""
增加装饰器----加快递归速度：
#建立一个装饰器（记录已经形成的dic)
def decorate(fun):
already_compuate=defaultdict()
def wrape(arg):
if arg in already_compuate:
result=already_compuate[arg]
else:
result=fun(arg)
already_compuate[arg]=result
return result
return wrape
solution={}
@decorate
def cut_mat(n):
max_price,max_split=max([(price_dic[n],0)]+[(cut_mat(i)+cut_mat(n-i),i) for i in range(1,n)],key=lambda x:x[0])
solution[n]=(n-max_split,max_split)#表示木材的长度：切成什么样的形式：8---->6+2
return max_price
display(cut_mat(15),solution)
###########结果：
45
{1: (1, 0),
2: (2, 0),
3: (3, 0),
4: (2, 2),
5: (3, 2),
6: (6, 0),
7: (6, 1),
8: (6, 2),
9: (6, 3),
10: (10, 0),
11: (11, 0),
12: (11, 1),
13: (11, 2),
14: (11, 3),
15: (13, 2)}
二、编辑距离
1、定义及原理
参考原文，并致谢该作者
编辑距离（Edit Distance），又称Levenshtein距离，是指两个字串之间，由一个转成另一个所需的最少编辑操作次数。许可的编辑操作包括将一个字符替换成另一个字符，插入一个字符，删除一个字符。一般来说，编辑距离越小，两个串的相似度越大。
例如将kitten一字转成sitting：
sitten （k→s）
sittin （e→i）
sitting （→g）
俄罗斯科学家Vladimir Levenshtein在1965年提出这个概念。
简单介绍如下图
2、代码实现
#建立一个装饰器（记录已经形成的dic)
def decorate(fun):
already_compuate=defaultdict()
@wraps
def wrape(*arg):
if arg in already_compuate:
result=already_compuate[arg]
else:
result=fun(*arg)
already_compuate[arg]=result
return result
return wrape
solution={}
solutio={}
@decorate
def edit_distance(str1,str2):
#如果两个字符中有一个为空字符，即空字符如何变成字符串或字符串如何变成空字符，则通过insert 或del 字符串的长度
if len(str1)==0:return len(str2)
if len(str2)==0:return len(str1)
#统计最后字符：
tail_str1=str1[-1]
tail_str2=str2[-1]
'''
利用动态规划的计算规则，减少字符串的长度比较：（i,j）<----（i-1,j）+1：del str1
(i,j) <---- (i,j-1) +1：insert str2
(i,j) <---- (i+1,j+1) +1(str1[-1]!=str2[-1]) or 0 (str1[-1]==str2[-1])
'''
candicate=[(edit_distance(str1[:-1],str2[:])+1,"del:{}".format(tail_str1)),(edit_distance(str1[:],str2[:-1])+1,"insert:{}".format(tail_str2)
)]
if tail_str1==tail_str2:
both_forward=(edit_distance(str1[:-1],str2[:-1])+0,'')
else: #直接将最后两个字符进行replace,即A--->B或B---->A：
both_forward=(edit_distance(str1[:-1],str2[:-1])+1,'sub')
candicate.append(both_forward)
min_distance,operation=min(candicate,key=lambda x:x[0])
solutio[(str1,str2)]=operation
return min_distance
####调用--并观察时间
start_time=time.time()
print(edit_distance("ABCRFT","AECFRETGF"))
print(solutio)
end_time=time.time()
print("时间是{0}".format(end_time-start_time))
三、Needleman-Wunsch 算法
1、定义
动态规划算法（dynamic programming algorithm)是最精确的一种序列比对算法，分为全局动态规划算法和局部动态规划算法。经典的全局动态规划算法由Needleman 和Wunsch 于1970年提出，成为Needleman-Wunsch算法，用于发现两条序列的全局水平上的相似性。经典的局部动态规划算法由Smith 和 Waterman 于1981年提出，成为Smith 和 Waterman于1981年提出，称为Smith-Waterman 算法，用于发现两条序列在局部水平上的相似性。
2、原理
比较两个序列的相似性，主要是建立两个矩阵，得分矩阵（score_matrix）和回溯矩阵（trace_back）,前者主要是表示序列A和序列B在碱基比对上相似的打分情况，后者主要是通过相应位置的情况建立indel or match or mismatch label。
1）得分矩阵–score_matrix
计算公式：
1、 当规定match=1,mismatch=-1,indel=-1时，对角线+1/-1，两侧（左侧或上侧）-1，所以找计算当前位置的分数的时候，一般建议找对角线，因为对角线一般较大，其次是其他
2、 计算的顺序是从左上角—>右下角，一般左上角直接赋值为0，（可以想象两端序列最初都有空白符，纪念空白符一样，即惩罚分数为0），第一行和第一列都是0，-1，-2，…（-1）*len(sequence),可以解释为其中一个序列为空白符后，转为另一个序列需要进行几次insert
import itertools
import copy
#seq1 = '*GCATGCU'
#seq2 = '*GATTACA'
seq1='*cat'
seq2='*actg'
seq1="*cat"
seq2="*actg"
#score_matrix:计算数值
#trace_back:回溯矩阵，用于进行回溯找序列
score_matrix=[[0 for i in range(len(seq1))] for row in range(len(seq2))]
trace_matrix=[[[] for i in range(len(seq1))] for row in range(len(seq2))]
for i in range(len(score_matrix[0])):
score_matrix[0][i]=i*(-1)
if i>0:
trace_matrix[0][i].append('left')#一定是append，赋值=会减少情况
for i in range(len(score_matrix)):
score_matrix[i][0]=i*(-1)
if i>0:
trace_matrix[i][0].append('up')
trace_matrix[0][0].append('done')
####output结果：
trace_matrix：
[[['done'], ['left'], ['left'], ['left']],
[['up'], [], [], []],
[['up'], [], [], []],
[['up'], [], [], []],
[['up'], [], [], []]]
import numpy as np
np.array(score_matrix)##因为list的2维展示不清楚，所以用np.array展示
array([[ 0, -1, -2, -3],
[-1, 0, 0, 0],
[-2, 0, 0, 0],
[-3, 0, 0, 0],
[-4, 0, 0, 0]])
###填充两个矩阵
score_dic={'match':1,'mismatch':-1,'gap':-1}
for i in range(1,len(score_matrix)):
for j in range(1,len(score_matrix[0])):
if seq1[j]==seq2[i]:
cell_scor=score_dic['match']
else:
cell_scor=score_dic['mismatch']
top_score=score_matrix[i-1][j]+score_dic['gap']
left_score=score_matrix[i][j-1]+score_dic['gap']
diag_score=score_matrix[i-1][j-1]+cell_scor
score=max(top_score,left_score,diag_score)
score_matrix[i][j]=score
if top_score==score:#一定要用if ,不能使用elif,else:否则少情况
trace_matrix[i][j].append('up')
if left_score==score:
trace_matrix[i][j].append('left')
if diag_score==score:
trace_matrix[i][j].append('diag')
#####################output结果：
trace_matrix：
[[['done'], ['left'], ['left'], ['left']],
[['up'], ['diag'], ['diag'], ['left']],
[['up'], ['diag'], ['up', 'left'], ['diag']],
[['up'], ['up'], ['diag'], ['diag']],
[['up'], ['up'], ['up', 'diag'], ['up']]]
score_matrix：
[[0, -1, -2, -3],
[-1, -1, 0, -1],
[-2, 0, -1, -1],
[-3, -1, -1, 0],
[-4, -2, -2, -1]]
3、计算方式和编辑距离相似又不相似：
相同之处： 都是通过前一个小规模计算后一个规模：D（i-1,j-1）or D(i,j-1) or D(i-1,j)—>D(i,j)；都是找到最有对齐（匹配方式）
不同之处： 编辑距离主要存在3种方式：insert(+1)、delete(+1)、replace(+0/1),其中数字表示存在几步，所以最终是求最小值，即一个单词变为另一单词最少的过程是多少；Needleman-Wunsch存在3种方式：indel(-1)、match(+1)、mismatch（-1），即insert和del进行了合并。
2）回溯矩阵
1、回溯矩阵，主要是矩阵大小和得分矩阵一样，如下图，主要是从右下角---->左上角进行回溯，就是找最大值（但是是找该位置是由前者（对角线，左边，上边）哪个得来的，下图种match=2，所以画×的地方很明显不能进行回溯，因为不是最优，现在网上很多blog在回溯过程不说清楚，所以本人前期也是走了很多弯路：总之一句话 谁计算出来的找谁
2、有的时候可能出现两种情况，说明最优途径不只是一条。
3、写序列的时候可以从左到右或从右到左，该位置是(-)或match 或mismatch,主要是看看该位置的箭头是怎样表示的
4、可以将回溯矩阵（表示为up，left ，diag)矩阵（如代码），更清楚
def sequence_finder(current_label,current_pointer):#回溯通过label(up or left or diag) 和pointer(坐标)操作
if current_label=="diag":
letter=[seq1[current_pointer[1]],seq2[current_pointer[0]]]
next_pointer=[current_pointer[0]-1,current_pointer[1]-1]
next_label=trace_matrix[next_pointer[0]][next_pointer[1]]
return letter,next_label,next_pointer
elif current_label=="left":
letter=[seq1[current_pointer[1]],'-']#左边序列插入空白符-
next_pointer=[current_pointer[0],current_pointer[1]-1]
next_label=trace_matrix[next_pointer[0]][next_pointer[1]]
return letter,next_label,next_pointer
else:#up
letter=['-',seq2[current_pointer[0]]]#顶端序列插入空白符-
next_pointer=[current_pointer[0]-1,current_pointer[1]]
next_label=trace_matrix[next_pointer[0]][next_pointer[1]]
return letter,next_label,next_pointer
def align_sequence_finder(current_label,current_pointer,sequence_ls):#目的在于匹配当前在trace_matrix的位置，因为trace_matrix表格中可能存在一个、两个、三个label
#当已经到达左上角--终点时:
if current_label[0]=='done':
sequence_ls=[sequence_ls[0][::-1],sequence_ls[1][::-1]]
return sequence_ls
elif len(current_label)==1:
letter,current_label,current_pointer=sequence_finder(current_label[0],current_pointer)##current_label[0]---》关键点
sequence_ls[0]+=letter[0]
sequence_ls[1]+=letter[1]
return align_sequence_finder(current_label,current_pointer,sequence_ls)#进行下一个递归
#['up', 'left'], [2, 2], ['t', 't']
elif len(current_label)==2:
label_1=current_label[0]
pointer_1=copy.deepcopy(current_pointer)#必须是深拷贝
sequence_ls_1=copy.deepcopy(sequence_ls)
letter1,label_1,pointer_1=sequence_finder(label_1,pointer_1)
sequence_ls_1[0]+=letter1[0]
sequence_ls_1[1]+=letter1[1]
label_2=current_label[1]
pointer_2=copy.deepcopy(current_pointer)#必须是深拷贝
sequence_ls_2=copy.deepcopy(sequence_ls)
letter2,label_2,pointer_2=sequence_finder(label_2,pointer_2)
sequence_ls_2[0]+=letter2[0]
sequence_ls_2[1]+=letter2[1]
return list(itertools.chain(align_sequence_finder(label_1,pointer_1,sequence_ls_1),
align_sequence_finder(label_2,pointer_2,sequence_ls_2)))
elif len(current_label)==3:
label_1=current_label[0]
pointer_1=copy.deepcopy(current_pointer)#必须是深拷贝
sequence_ls_1=copy.deepcopy(sequence_ls)
letter,label_1,pointer_1=sequence_finder(label_1,pointer_1)
sequence_ls_1[0]+=letter[0]
sequence_ls_1[1]+=letter[1]
label_2=current_label[1]
pointer_2=copy.deepcopy(current_pointer)#必须是深拷贝
sequence_ls_2=copy.deepcopy(sequence_ls)
letter,label_2,pointer_2=sequence_finder(label_2,pointer_2)
sequence_ls_2[0]+=letter[0]
sequence_ls_2[1]+=letter[1]
label_3=current_label[1]
pointer_3=copy.deepcopy(current_pointer)#必须是深拷贝
sequence_ls_3=copy.deepcopy(sequence_ls)
letter,label_3,pointer_3=sequence_finder(label_3,pointer_3)
sequence_ls_3[0]+=letter[0]
sequence_ls_3[1]+=letter[1]
return list(itertools.chain(align_sequence_finder(label_1,pointer_1,sequence_ls_1),
align_sequence_finder(label_2,pointer_2,sequence_ls_2),
align_sequence_finder(label_3,pointer_3,sequence_ls_3))
)
#回溯起点：
#[seq2_index.seq1_index]
start_label=trace_matrix[len(seq2)-1][len(seq1)-1]
start_pointer=[len(seq2)-1,len(seq1)-1]
start_letter=['','']
display(start_label,start_pointer,start_letter)
#调用方法
sequence_final_ls=align_sequence_finder(start_label,start_pointer,start_letter)
for i in range(0,len(sequence_final_ls),2):
print([sequence_final_ls[i],sequence_final_ls[i+1]])
#########################最终结果：
['ca-t-', '-actg']
['-cat-', 'ac-tg']
3、完整代码实现
import copy
import itertools
seq1="*cat"
seq2="*actg"
#score_matrix:计算数值
#trace_back:回溯矩阵，用于进行回溯找序列
score_matrix=[[0 for i in range(len(seq1))] for row in range(len(seq2))]
trace_matrix=[[[] for i in range(len(seq1))] for row in range(len(seq2))]
for i in range(len(score_matrix[0])):
score_matrix[0][i]=i*(-1)
if i>0:
trace_matrix[0][i].append('left')#一定是append，赋值=会减少情况
for i in range(len(score_matrix)):
score_matrix[i][0]=i*(-1)
if i>0:
trace_matrix[i][0].append('up')
trace_matrix[0][0].append('done')
score_dic={'match':1,'mismatch':-1,'gap':-1}
for i in range(1,len(score_matrix)):
for j in range(1,len(score_matrix[0])):
if seq1[j]==seq2[i]:
cell_scor=score_dic['match']
else:
cell_scor=score_dic['mismatch']
top_score=score_matrix[i-1][j]+score_dic['gap']
left_score=score_matrix[i][j-1]+score_dic['gap']
diag_score=score_matrix[i-1][j-1]+cell_scor
score=max(top_score,left_score,diag_score)
score_matrix[i][j]=score
if top_score==score:#一定要用if ,不能使用elif,else:否则少情况
trace_matrix[i][j].append('up')
if left_score==score:
trace_matrix[i][j].append('left')
if diag_score==score:
trace_matrix[i][j].append('diag')
def sequence_finder(current_label,current_pointer):#回溯通过label(up or left or diag) 和pointer(坐标)操作
if current_label=="diag":
letter=[seq1[current_pointer[1]],seq2[current_pointer[0]]]
next_pointer=[current_pointer[0]-1,current_pointer[1]-1]
next_label=trace_matrix[next_pointer[0]][next_pointer[1]]
return letter,next_label,next_pointer
elif current_label=="left":
letter=[seq1[current_pointer[1]],'-']#左边序列插入空白符-
next_pointer=[current_pointer[0],current_pointer[1]-1]
next_label=trace_matrix[next_pointer[0]][next_pointer[1]]
return letter,next_label,next_pointer
else:#up
letter=['-',seq2[current_pointer[0]]]#顶端序列插入空白符-
next_pointer=[current_pointer[0]-1,current_pointer[1]]
next_label=trace_matrix[next_pointer[0]][next_pointer[1]]
return letter,next_label,next_pointer
def align_sequence_finder(current_label,current_pointer,sequence_ls):#目的在于匹配当前在trace_matrix的位置，因为trace_matrix表格中可能存在一个、两个、三个label
#当已经到达左上角--终点时:
if current_label[0]=='done':
sequence_ls=[sequence_ls[0][::-1],sequence_ls[1][::-1]]
return sequence_ls
elif len(current_label)==1:
letter,current_label,current_pointer=sequence_finder(current_label[0],current_pointer)##current_label[0]---》关键点
sequence_ls[0]+=letter[0]
sequence_ls[1]+=letter[1]
return align_sequence_finder(current_label,current_pointer,sequence_ls)#进行下一个递归
#['up', 'left'], [2, 2], ['t', 't']
elif len(current_label)==2:
label_1=current_label[0]
pointer_1=copy.deepcopy(current_pointer)#必须是深拷贝
sequence_ls_1=copy.deepcopy(sequence_ls)
letter1,label_1,pointer_1=sequence_finder(label_1,pointer_1)
sequence_ls_1[0]+=letter1[0]
sequence_ls_1[1]+=letter1[1]
label_2=current_label[1]
pointer_2=copy.deepcopy(current_pointer)#必须是深拷贝
sequence_ls_2=copy.deepcopy(sequence_ls)
letter2,label_2,pointer_2=sequence_finder(label_2,pointer_2)
sequence_ls_2[0]+=letter2[0]
sequence_ls_2[1]+=letter2[1]
return list(itertools.chain(align_sequence_finder(label_1,pointer_1,sequence_ls_1),
align_sequence_finder(label_2,pointer_2,sequence_ls_2)))
elif len(current_label)==3:
label_1=current_label[0]
pointer_1=copy.deepcopy(current_pointer)#必须是深拷贝
sequence_ls_1=copy.deepcopy(sequence_ls)
letter,label_1,pointer_1=sequence_finder(label_1,pointer_1)
sequence_ls_1[0]+=letter[0]
sequence_ls_1[1]+=letter[1]
label_2=current_label[1]
pointer_2=copy.deepcopy(current_pointer)#必须是深拷贝
sequence_ls_2=copy.deepcopy(sequence_ls)
letter,label_2,pointer_2=sequence_finder(label_2,pointer_2)
sequence_ls_2[0]+=letter[0]
sequence_ls_2[1]+=letter[1]
label_3=current_label[1]
pointer_3=copy.deepcopy(current_pointer)#必须是深拷贝
sequence_ls_3=copy.deepcopy(sequence_ls)
letter,label_3,pointer_3=sequence_finder(label_3,pointer_3)
sequence_ls_3[0]+=letter[0]
sequence_ls_3[1]+=letter[1]
return list(itertools.chain(align_sequence_finder(label_1,pointer_1,sequence_ls_1),
align_sequence_finder(label_2,pointer_2,sequence_ls_2),
align_sequence_finder(label_3,pointer_3,sequence_ls_3))
)
#回溯起点：
#[seq2_index.seq1_index]
start_label=trace_matrix[len(seq2)-1][len(seq1)-1]
start_pointer=[len(seq2)-1,len(seq1)-1]
start_letter=['','']
display(start_label,start_pointer,start_letter)
#调用方法
sequence_final_ls=align_sequence_finder(start_label,start_pointer,start_letter)
for i in range(0,len(sequence_final_ls),2):
print([sequence_final_ls[i],sequence_final_ls[i+1]])
#########################最终结果：
['ca-t-', '-actg']
['-cat-', 'ac-tg']
扫一扫
837
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
