In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
In [ ]:
import matplotlib as pl
from msa2qubo import *
from qubo2msa import *
import numpy as np
In [ ]:
m2q = Msa2Qubo(input="./seq3.fa", output="t2/qmsa.qubo", P=1, verbose=False, delta=2.0, l0=1.0, l1=0.1, l2=10.0)
In [ ]:
m2q.run()
In [ ]:
import math
m2q.bvc.createW()
W = m2q.bvc.W()
In [ ]:
X = m2q.bvc.x()
print("seqs=",X)
print("# seqs=",len(X))
nr = 0
ny = 0
nz = 0
for k in range(0,len(X)-1):
for q in range(k+1, len(X)):
for i in range(0, len(X[k])):
for j in range(0, len(X[q])):
nr+=1
ny+=1
nz+=1
print("#R=",nr)
print("#Y=",ny)
print("#Z=",nz)
print("#Additional Vars=",nr+ny+nz)
In [ ]:
numPosVars = m2q.bvc.get_NbPositioningVars(intmode=True)
numGapVars = m2q.bvc.get_NbGapVars(intmode=True)
print(numPosVars)
print(numGapVars)
In [ ]:
e2m = np.zeros((numPosVars+numGapVars+nr+ny+nz,numPosVars+numGapVars+nr+ny+nz))
print(e2m.shape)
In [ ]:
solution_vector = np.zeros((1, numPosVars+numGapVars+nr+ny+nz))
In [ ]:
In [ ]:
# Move to the beggining of the matrix where the positions for R,Y and Z are
rMatPos = numPosVars+numGapVars
# Place the R*W variables on the R part of the matrix
pbK = 0
pbQ = 0
posI = 0
posJ = 0
for k in range(0,len(X)-1):
pbK += len(X[k])
for q in range(k+1, len(X)):
pbQ += len(X[q])
for i in range(0, len(X[k])):
posI = 0
posI = pbK+i
for j in range(0, len(X[q])):
posJ = 0
posJ = pbQ+j
matPos = rMatPos + pbK + pbQ + posI + posJ
e2m[matPos][matPos] = m2q.bvc.w(i,j,k,q)
In [ ]:
# Place the X*Y, X*Z variables in the matrix multiplied by l_2 and W_ijkq
In [ ]:
yMatPos = rMatPos + nr
zMatPos = yMatPos + ny
pbK = 0
pbQ = 0
posI = 0
posJ = 0
for k in range(0,len(X)-1):
for q in range(k+1, len(X)):
for i in range(0, len(X[k])):
# This loop would need to be extended for the binary version
posI = 0
posI = pbK+i
for j in range(0, len(X[q])):
# This loop would need to be extended for the binary version
posJ = 0
posJ = pbQ+j
YmatPos = yMatPos + pbK + pbQ + posI + posJ
ZmatPos = zMatPos + pbK + pbQ + posI + posJ
e2m[pbK+i][YmatPos] = -1. * m2q.bvc.l2() * m2q.bvc.w(i,j,k,q)
e2m[pbQ+j][YmatPos] = 2. * m2q.bvc.l2() * m2q.bvc.w(i,j,k,q)
e2m[pbQ+j][ZmatPos] = -1. * m2q.bvc.l2() * m2q.bvc.w(i,j,k,q)
pbQ += len(X[q])
pbK += len(X[k])
In [ ]:
# Place the bounds on Y and Z based on the McCormick envelopes
#YunderMin = 0
#YunderMax = X_ki + Rijkq*Lmax - Lmax
#YoverMin = X_ki
#YoverMax = Rijkq*Lmax
# This goes on the integer problem constraints matrix/array
In [ ]:
In [ ]:
m = m2q.bvc.m()
nbR = nr
nbY = ny*m
nbZ = nz*m
print("nbR", nbR)
print("nbY", nbY)
print("nbZ", nbZ)
In [ ]:
delta = 1
bVarRPos = m2q.bvc.get_NbPositioningVars() + m2q.bvc.get_NbGapVars()
bVarYPos = bVarRPos + nr # No need to transform R to binary as it only either (0,1)
# Since R is either 0,1 in the following loops need to take that into account...
bVarZPos = bVarYPos + ny*m
bRMatPos = bVarRPos
bYMatPos = bVarYPos
bZMatPos = bVarZPos
matDims = bVarZPos+nz*m
e2bm = np.zeros((matDims, matDims))
pbK = 0
pbQ = 0
posI = 0
posJ = 0
piK = 0
piQ = 0
x_k = 0
x_q = 0
for k in range(len(X)-1):
pbQ = 0
piQ = 0
size_iq = 0
size_bq = 0
x_q = x_k + len(X[k]) * m
for q in range(k+1, len(X)):
size_ii = 0
size_bi = 0
for i in range(len(X[k])):
x_ki = x_k + (i * m)
size_ij = 0
size_bj = 0
for j in range(len(X[q])):
wl2 = m2q.bvc.w(i=i,j=j,k=k,q=q) * m2q.bvc.l2()
x_qj = x_q + (j * m)
i_idx = piK + size_iq + size_ii + size_ij
b_idx = pbK + size_bq + size_bi + size_bj
r_kqij = bRMatPos + i_idx
e2bm[r_kqij][r_kqij] += wl2
for bi in range(0,m):
for bj in range(0, m):
x_kia = x_ki + bi
x_qja = x_qj + bi
y_kqija = bYMatPos+b_idx + bi
y_kqijb = bYMatPos+b_idx + bj
z_kqija = bZMatPos+b_idx + bi
z_kqijb = bZMatPos+b_idx + bj
# -Yijkq*Xki
e2bm[x_kia][y_kqijb] += wl2 * -1. * 2**(bi+bj)
# 3d * Yijkq
if bi == bj:
e2bm[y_kqija][y_kqijb] += wl2 * 3.*delta * 2**(bi+bj)
elif bi < bj:
e2bm[y_kqija][y_kqijb] += wl2 * 3.*delta * 2**(bi+bj)*2
# d * Rijqk*Xki
e2bm[x_kia][r_kqij] += wl2 * 1.*delta * 2**(bi+bj)
# -2d * Yijkq*Rijkq
e2bm[r_kqij][y_kqijb] += wl2 * -2.*delta * 2**(bi+bj)
# 2d * Yijkq*Xki
e2bm[x_kia][y_kqijb] += wl2 * 2.*delta * 2**(bi+bj)
# -4d * Yijkq*Xqj
e2bm[x_qja][y_kqijb] += wl2 * -4.*delta * 2**(bi+bj)
# Zijkq*Xqj
e2bm[x_qja][z_kqijb] += wl2 * 1. * 2**(bi+bj)
# -3d * Zijkq
if bi == bj:
e2bm[z_kqija][z_kqijb] += wl2 * -3.*delta * 2**(bi+bj)
elif bi < bj:
e2bm[z_kqija][z_kqijb] += wl2 * -3.*delta * 2**(bi+bj)*2
# -d * Rijkq*Xqj
e2bm[x_qja][r_kqij] += wl2 * -1.*delta * 2**(bi+bj)
# 2d * Zijqk*Xqj
e2bm[x_qja][z_kqijb] += wl2 * 2.*delta * 2**(bi+bj)
# 2d * Zijkq*Rijkq
e2bm[r_kqij][z_kqijb] += wl2 * 2.*delta * 2**(bi+bj)
size_ij += 1
size_bj += m
size_ii += size_ij
size_bi += size_bj
size_iq += size_ii
size_bq += size_bi
x_q += len(X[q]) * m
pbK += size_bq
piK += size_iq
x_k += len(X[k]) * m
In [ ]:
G = np.zeros((matDims,matDims,3))
G[e2bm!=0] = [0,0,0]
G[e2bm==0] = [1,1,1]
plt.figure(figsize=(10,10))
#plt.annotate(s="G", xy=((m2q.bvc.get_NbPositioningVars())*0.5,-1))
#plt.annotate(s="G", xy=(matDims,(m2q.bvc.get_NbPositioningVars())*0.5))
plt.axvline(m2q.bvc.get_NbPositioningVars()-0.5)
plt.axhline(m2q.bvc.get_NbPositioningVars()-0.5)
plt.axvline(m2q.bvc.get_NbPositioningVars()+m2q.bvc.get_NbGapVars()-0.5)
plt.axhline(m2q.bvc.get_NbPositioningVars()+m2q.bvc.get_NbGapVars()-0.5)
plt.annotate(s="R", xy=((bVarRPos+bVarYPos)*0.5,-1))
plt.annotate(s="R", xy=(matDims,(bVarRPos+bVarYPos)*0.5))
plt.axvline(bVarRPos-0.5)
plt.axhline(bVarRPos-0.5)
plt.annotate(s="Y", xy=((bVarYPos+bVarZPos)*0.5,-1))
plt.annotate(s="Y", xy=(matDims,(bVarYPos+bVarZPos)*0.5))
plt.axvline(bVarYPos-0.5)
plt.axhline(bVarYPos-0.5)
plt.annotate(s="Z", xy=((bVarZPos+matDims)*0.5,-1))
plt.annotate(s="Z", xy=(matDims+0.11,(bVarZPos+matDims)*0.5))
plt.axvline(bVarZPos-0.5)
plt.axhline(bVarZPos-0.5)
plt.imshow(G, cmap='Greys', interpolation='nearest')
plt.show()
In [ ]: