import copy


rule_file=open('super_rule_value.txt','r')
output=open('super_rule_HB_value.txt','w')


aa_participate_threshold=8




all=[]
rule=[]
ws=9
#----------------------------------------------------------------
#Blosum62 matrix
aa_participate_threshold=8
blosum62=[     [4,-1,-2,-2,0,-1,-1,0,-2,-1,-1,-1,-1,-2,-1,1,0,-3,-2,0],
               [-1,5,0,-2,-3,1,0,-2,0,-3,-2,2,-1,-3,-2,-1,-1,-3,-2,-3],
               [-2,0,6,1,-3,0,0,0,1,-3,-3,0,-2,-3,-2,1,0,-4,-2,-3],
               [-2,-2,1,6,-3,0,2,-1,-1,-3,-4,-1,-3,-3,-1,0,-1,-4,-3,-3],
               [0,-3,-3,-3,9,-3,-4,-3,-3,-1,-1,-3,-1,-2,-3,-1,-1,-2,-2,-1],
               [-1,1,0,0,-3,5,2,-2,0,-3,-2,1,0,-3,-1,0,-1,-2,-1,-2],
               [-1,0,0,2,-4,2,5,-2,0,-3,-3,1,-2,-3,-1,0,-1,-3,-2,-2],
               [0,-2,0,-1,-3,-2,-2,6,-2,-4,-4,-2,-3,-3,-2,0,-2,-2,-3,-3],
               [-2,0,1,-1,-3,0,0,-2,8,-3,-3,-1,-2,-1,-2,-1,-2,-2,2,-3],
               [-1,-3,-3,-3,-1,-3,-3,-4,-3,4,2,-3,1,0,-3,-2,-1,-3,-1,3],
               [-1,-2,-3,-4,-1,-2,-3,-4,-3,2,4,-2,2,0,-3,-2,-1,-2,-1,1],
               [-1,2,0,-1,-3,1,1,-2,-1,-3,-2,5,-1,-3,-1,0,-1,-3,-2,-2],
               [-1,-1,-2,-3,-1,0,-2,-3,-2,1,2,-1,5,0,-2,-1,-1,-1,-1,1],
               [-2,-3,-3,-3,-2,-3,-3,-3,-1,0,0,-3,0,6,-4,-2,-2,1,3,-1],
               [-1,-2,-2,-1,-3,-1,-1,-2,-2,-3,-3,-1,-2,-4,7,-1,-1,-4,-3,-2],
               [1,-1,1,0,-1,0,0,0,-1,-2,-2,0,-1,-2,-1,4,1,-3,-2,-2],
               [0,-1,0,-1,-1,-1,-1,-2,-2,-1,-1,-1,-1,-2,-1,1,5,-2,-2,0],
               [-3,-3,-4,-4,-2,-2,-3,-2,-2,-3,-2,-3,-1,1,-4,-3,-2,11,2,-3],
               [-2,-2,-2,-3,-2,-1,-2,-3,2,-1,-1,-2,-1,3,-3,-2,-2,2,7,-1],
               [0,-3,-3,-3,-1,-2,-2,-3,-3,3,1,-2,1,-1,-2,-2,0,-3,-1,4]
         ]

substitude=[19,10,9,12,13,17,18,7,0,14,15,16,4,8,1,11,5,6,2,3]

hydro=['hydro value',2.3, 2.2, 3.1, 1.1, 2.5, 1.5, 0.08, 0.67, 1.0, -0.29, -1.1, -0.75, 0.17, -1.7, -7.5, -4.6, -2.6, -2.9, -2.7, -3.0]
iup=['iup value',-0.7055, -0.4385, -0.515, -0.4765, -0.497, -0.257, 0.0825, 0.6675, -0.275, 1.117, 0.2965, 0.145, -0.1255, 0.135, -0.179, -0.0495, -0.055, -0.2745, 0.479, 0.4645]	   



#----------------Read in rules------------------

for line in rule_file.readlines():
    aa=[]
    count=0
    for i in line.split():
        if count==0:
            aa.append(i)
            count+=1
        else:
            aa.append(float(i))
    rule.append(aa)

#------------------------------------------------------------
print " "

 

#result compute: HSSP-BLOSUM62 Value

for rule_num in range(len(rule)):
    #print rule_num
    amino_acid=["V","L","I","M","F","W","Y","G","A","P","S","T","C","H","R","K","Q","E","N","D"]
    avg_HSSP_BLOSUM62=0

    for i in range(ws):
        amino_acid_participate=[]
        amino_acid_participate_value=[]

        for j in range(1,21):
            if rule[rule_num][i*20+j] >= aa_participate_threshold:
                amino_acid_participate.append(j)
                amino_acid_participate_value.append(rule[rule_num][i*20+j])

        if len(amino_acid_participate)==0:
            HSSP_BLOSUM62_result=0
            avg_HSSP_BLOSUM62+=0
            
        elif len(amino_acid_participate)==1:
            #if amino_acid_participate_value[0] >= 10:
            #    HSSP_BLOSUM62_result= blosum62[substitude[amino_acid_participate[0]-1]][substitude[amino_acid_participate[0]-1]]
            #else:
            #    HSSP_BLOSUM62_result= (blosum62[substitude[amino_acid_participate[0]-1]][substitude[amino_acid_participate[0]-1]]) / 2
            #avg_HSSP_BLOSUM62+=HSSP_BLOSUM62_result



            HSSP_BLOSUM62_result= blosum62[substitude[amino_acid_participate[0]-1]][substitude[amino_acid_participate[0]-1]]
            avg_HSSP_BLOSUM62+=HSSP_BLOSUM62_result
            ##HSSP_BLOSUM62_result=0
            ##avg_HSSP_BLOSUM62+=0

        else:
            upper=0
            lower=0
            for k in range(len(amino_acid_participate)):
                for l in range(k+1,len(amino_acid_participate)):
                    lower+=amino_acid_participate_value[k]*amino_acid_participate_value[l]
                    upper+=amino_acid_participate_value[k]*amino_acid_participate_value[l]*blosum62[substitude[amino_acid_participate[k]-1]][substitude[amino_acid_participate[l]-1]]

            HSSP_BLOSUM62_result=upper/lower                    
            avg_HSSP_BLOSUM62+=HSSP_BLOSUM62_result
            
        print "HSSP_BLOSUM62 value for position",i,"=",HSSP_BLOSUM62_result


    avg_HSSP_BLOSUM62/=ws
    print rule_num, '\t', avg_HSSP_BLOSUM62
    output.write(str(avg_HSSP_BLOSUM62)+'\n')