scs.py

from fasta import readFASTA
import numpy as np
import string
import operator as op

with open("rosalind_scsp.txt","rU") as f:
    A=f.readline().strip()
    B=f.readline().strip()


#A='ATCTGAT'
#B='TGCATA'
nA=len(A)
nB = len(B)


dp=np.zeros((nB,nA))
pointerx=np.zeros((nB,nA),dtype=int)
pointery=np.zeros((nB,nA),dtype=int)
savedp=0
savedi=0
savedj=0

if A[0]==B[0]:
    dp[0][0]=1

for i in range(1,nA):
    pointerx[0][i]=-1
    if B[0]==A[i]:
        dp[0][i]=1
        savedi=i
        savedp=1
        pointerx[0][i]=-1
        pointery[0][i]=-1
    else:
        dp[0][i]=dp[0][i-1]
        pointerx[0][i]=-1

for j in range(1,nB):
    pointery[j][0]=-1
    if B[j]==A[0]:
        dp[j][0]=1
        savedj=j
        savedp=1
        pointerx[j][0]=-1
        pointery[j][0]=-1
    else:
        dp[j][0]=dp[j-1][0]
        pointery[j][0]=-1

for i in range(1,nA):
    for j in range(1,nB):
        if A[i]==B[j]:
           dp[j][i]=dp[j-1][i-1]+1
           (dx,dy)=(-1,-1)
        else:
            t=[(dp[j-1][i],(-1,0)),(dp[j][i-1],(0,-1))]
            (t,(dx,dy))=max(t,key=lambda x: op.getitem(x,0))
            dp[j][i]=dp[j+dx][i+dy]       
        pointerx[j][i],pointery[j][i]=dx,dy
        if dp[j][i]>savedp:
            savedp=dp[j][i]
            savedi,savedj=i,j

            
answer=[]
locs=[]
i,j = savedi,savedj
while True:
##   print i,j,pointerx[j][i],pointery[j][i],dp[j][i]
    if pointerx[j][i]==pointery[j][i] and pointerx[j][i]==-1:
        answer.append(A[i])
        locs.append((j,i))
    i,j=i+pointery[j][i],j+pointerx[j][i]  
    if i<0 or j<0:
        break
    
print ''.join(reversed(answer))
print locs[::-1]  
lastu=0
lastv=0
final=[]
for x in locs[::-1]:
    (u,v)=x
    final.append(B[lastu:u])
    final.append(A[lastv:v])
    final.append(B[u])
    print u,v,final
    lastu=u+1
    lastv=v+1
final.append(A[lastu:])
final.append(B[lastv:])
#print final
#print A
#print B
print ''.join(final)