# -*- coding: utf-8 -*-
"""
Created on Wed Apr 12 19:20:59 2023

@author: bzchen
"""
from contextlib import redirect_stdout
###############################################################################
#specify input file names 
###############################################################################
ncacx_filename = "Consistency450MillionMCassignRun.txt" #pls update the ncacx data file name
num_runs = 10 #ncocx_filename = "ncocx_b.txt" #pls update the ncacx data file name
res_num = 230
protein_seq = "RSV_seq.txt"
###############################################################################
asgn_filename = "Consistency.txt"
################################################################################
###############################################################################
#define some functions frequently used
###############################################################################
#function to extract only the letter part of a residue type asgn, not suitable 
#if you want to extract a single letter from an ambiguous residue type asgm
def extract_letter(rtyp):
    #import re
    kresname=''
    for i in rtyp:
        if i.isalpha():
            kresname = "".join([kresname, i])
    return kresname
#end of extract letter function definition

#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

#end of frequent used function definition
###############################################################################
###############################################################################
#read in protein sequence
###############################################################################

print('Read ncacx input from ',ncacx_filename,'.')
with open(ncacx_filename, "r") as file1:  
    #len(file1.readlines())
    consistency = [0 for i in range(0, res_num+1)]
    count = 0
    correct_asgn = 0
    #print(file1)
    for line in file1:
        kline=line.split()
        reference = int(kline[0])
        flag = 0 
        for i in range(0,num_runs):
            knum = int(kline[i])
            if knum != reference:
                flag = 1
        if flag == 0 and reference != 0:
            consistency[count] += 1
                
        #if consistency[count] == num_runs:
            correct_asgn += 1
            
        count += 1
print('consistantly asgned number of residues are ',correct_asgn)                

