все!
Я пытаюсь реализовать в Python алгоритм McCaskill (http://rna.informatik.uni -freiburg.de / Teaching / index.jsp? ToolName = McCaskill) и столкнулись с некоторыми проблемами с заполнением матрицы Q и Qbp. Мой код:
import sys
import numpy as np
def initialize_old(N, q = 0.):
DP = np.empty((N,N))
DP[:] = np.NAN
for i in range(N):
DP[i][i] = q
if i != 0:
DP[i][i-1] = q
return DP
min_loop_length = 1 #Minimal loop length
RT = 1 #'Normalized' temperature
Eij = -2 #Energy weight of base pair
Seq = 'GUUCC' #for instance consider such seq
L = len(Seq)
Q = initialize_old(L, 1.0)
Qbp = initialize_old(L)
p = initialize_old(L)
Q = np.c_[np.zeros(L),Q]
Q[0][0] = 1.0
Qbp = np.c_[np.zeros(L),Qbp]
def Eps_ij(i,j):
if (Seq[i].upper(),Seq[j].upper()) in (('A','U'),('U','A'),('C','G'),('G','C'),('G','U')):
Qbp[i][j] = Q[i+1][j-1] * float(np.exp(-Eij / RT))
return float(Qbp[i][j])
return 0.0
def McCaskill():
for k in range(1, L+1):
for i in range(L-k+1):
j = i + k -1
sum = 0.0
for h in range(i,j-1):
sum += Q[i][h-1] * Eps_ij(h,j)
Q[i][j] = Q[i][j-1] + sum
if __name__ == '__main__':
McCaskill()
print(Q)
print(Qbp)
И вывод был:
Q = [[ 1. 1. nan nan nan nan]
[ 0. 1. 1. 1. 1. nan]
[ 0. nan 1. 1. nan nan]
[ 0. nan nan 1. 1. nan]
[ 0. nan nan nan 1. 1.]]
Qbp = [[ 0. 0. 7.3890561 7.3890561 7.3890561 nan]
[ 0. 0. 0. nan nan nan]
[ 0. nan 0. 0. nan nan]
[ 0. nan nan 0. 0. nan]
[ 0. nan nan nan 0. 0. ]]
В связи с отсутствием каких-либо подробностей начальных условий и подробностей о диапазоне рекурсии я не могу понять, как заполнять оба массива шаг за шагом.
Буду очень благодарен за любые советы и идеи!
Заранее спасибо!
С наилучшими пожеланиями!