# -*- coding: utf-8 -*-
"""
Created on Sun Oct 22 07:42:21 2023

@author: james
"""
#used to find the sequentially allocated positions for each of the RTAs in input from NCACX
from contextlib import redirect_stdout

ncacxin_filename = "ncacx_ambi.txt"#input RTAs from NCACX as the initial reference
ncacxout_filename = "ncacx4nextrd.txt"#output sequential assignments for RTAs in NCACX
position_filename = "position.txt"
total_col = 10

#function to extract only the number part of a residue type asgn
def extract_number(rtyp):
    #import re
    kresnum=''
    for i in rtyp:
        if i.isdigit():
            kresnum = "".join([kresnum, i])
    return kresnum
#end of extract number function

#print('Read reference ncacx input from ',ncacxin_filename,'.')
with open(ncacxin_filename, mode='r') as fp:
    count = 0
    print(fp)
    for line in fp:
        count += 1
        #Read in the first row, how many residue type asgn in ncacx and
        #how many fre columns in each residue type asgn
        if count == 1:
            line1=line.split( )
            reference_RTAnum = int(line1[0])                
            #nfreq_nca = int(line1[1])#actual number of freq used for match finding
            reference_CSnum = int(line1[1])#maximum number of freq input present, the max is 7 (including amide nitrogen and C', this means up to C-e is allowed in input, sufficient to differentiate any residue types)

#print('Read output ncacx input from ',ncacxout_filename,'.')
with open(ncacxout_filename, mode='r') as fp:
    count = 0
    print(fp)
    for line in fp:
        count += 1
        #Read in the first row, how many residue type asgn in ncacx and
        #how many fre columns in each residue type asgn
        if count == 1:
            line1=line.split( )
            out_RTAnum = int(line1[0])                
            #nfreq_nca = int(line1[1])#actual number of freq used for match finding
            out_CSnum = int(line1[1])#maximum number of freq input present, the max is 7 (including amide nitrogen and C', this means up to C-e is allowed in input, sufficient to differentiate any residue types)

valid_CS = [[0 for i in range(reference_CSnum)] for j in range(reference_RTAnum+1)]
reference_CS =[[1e6 for i in range(reference_CSnum)] for j in range(reference_RTAnum+1)]
out_CS =[[0 for i in range(reference_CSnum+1)] for j in range(out_RTAnum+1)]
position =[[0 for i in range(2)] for j in range(out_RTAnum+1)]   

print('Read reference ncacx input from ',ncacxin_filename,'.')
with open(ncacxin_filename, mode='r') as fp:
    count = 0
    print(fp)
    for line in fp:
        count += 1
        #Read in the first row, how many residue type asgn in ncacx and
        #how many fre columns in each residue type asgn
        if count > 1:
            kline=line.split()
            num_valid = 0
            for i in range(reference_CSnum):
                reference_CS[count-1] [i]= float(kline[i]) 
                if float(kline[i]) < 1000:
                    num_valid += 1
            valid_CS[count-1] = num_valid
                
print('Read out ncacx input from ',ncacxout_filename,'.')
with open(ncacxout_filename, mode='r') as fp:
    count = 0
    print(fp)
    for line in fp:
        count += 1
        #Read in the first row, how many residue type asgn in ncacx and
        #how many fre columns in each residue type asgn
        if count > 1:
            kline=line.split()
            for i in range(reference_CSnum):
                out_CS[count-1][i] = float(kline[i])
            out_CS[count-1][reference_CSnum] = extract_number(kline[total_col-1])
                
with open(position_filename, mode='w') as f: 
    with redirect_stdout(f): 
        for i in range(reference_RTAnum+1):
            position[i][0] = i
            find = 0 
            for j in range(out_RTAnum+1):
                if find == 0:
                    flag = 0
                    if reference_CS[i][0] == out_CS[j][0] and reference_CS[i][1] == out_CS[j][1] and reference_CS[i][2] == out_CS[j][2] and reference_CS[i][3] == out_CS[j][3]:#if reference_CS[i][k] == out_CS[j][k]: #and reference_CS[j][k] != 1000000:
                       #find = 1     #flag += 1 
                    #if flag == valid_CS[i]:
                        position[i][1] = out_CS[j][reference_CSnum]
                        find = 1 
            #if i == 11:
                #print('position 11 is', position[i])
            if find == 0 or position[i][1] == '':
                position[i][1] = 1000000 
        
        count = 0
        for i in range(reference_RTAnum+1):
            print(position[i][0],' ',position[i][1])
            
            if int(position[i][0]) == int(position[i][1]) and position[i][0] != 0:
                count += 1 
        print(count,' residues are correctly assigned ')
print(count,' residues are correctly assigned ')  
        
                 
   

