import re
import sys

"""Read in the _hmmparse.out(generated by HMMparse.py) and _true_coords.txt(generated by the simulator), 
compare and analyze the two groups of coordinates in these two files, and print out the result:
TP, TN, FP, FN, Sensitivity, Trims etc.
"""

def read_params(hmmparse, true_coords):
    fileptr1 = open(hmmparse)
    num_passed_TEs = 0    # store the number of found TEs
    chunks_list = []      # store the TEs coords
                          # such as [[950, 1000], ...]
    for line in fileptr1:
        line = line.rstrip()
        arr = re.split("\s+", line)   # split the line to an array
        
        if arr[0] == "The":           # get the number of passed TEs (first line)
            num_passed_TEs = int(arr[-1])

        elif arr[0][0] == 'c':        # for each TE (chunk+number)
            chunk_num = arr[0]
            line = fileptr1.next()
            line = line.rstrip()
            arr2 = re.split("\s+", line)
            chunks_list.append([int(arr2[0]), int(arr2[1])])

        else:
            continue
    fileptr1.close()
    chunks_list = sorted(chunks_list, key = lambda chunk: chunk[0])   # sort the detected intervals

    fileptr2 = open(true_coords)
    true_coords = []                  # store the true coordinates to a list
    for line in fileptr2:
        line = line.rstrip()
        arr = re.split("\s+", line)
        true_coords.append([int(arr[0]), int(arr[1])])

    #print "chunks_list: ", chunks_list
    #print "true_coords: ", true_coords
    return [chunks_list, true_coords]


def result_analyser(result, threshhold):
    chunks_list = result[0]
    new_chunks_list = []
    true_coords = result[1]
    miss_itv = []               # store the missed intervals when combining chunks_list
    i = 0
    if len(result[0])==1:       # if there's only one chunk
        new_chunks_list.append([chunks_list[0][0], chunks_list[0][1]])
    else:                       # if there are more than 1 chunks
        while i<(len(chunks_list)-1):     # combine the close intervals
            i2 = i+1
            flag = False
            while chunks_list[i2][0] - chunks_list[i2-1][1] <= 20:
                flag = True
                miss_itv.append([chunks_list[i2-1][1], chunks_list[i2][0]])
                if i2<(len(chunks_list)-1):
                    i2 = i2 + 1
                else:
                    break
            if i2==(len(chunks_list)-1) and flag == True:
                new_chunks_list.append([chunks_list[i][0], chunks_list[i2][1]])
            elif i2==(len(chunks_list)-1) and flag == False:
                new_chunks_list.append([chunks_list[i][0], chunks_list[i2-1][1]])
                new_chunks_list.append([chunks_list[i2][0], chunks_list[i2][1]])
            else:
                new_chunks_list.append([chunks_list[i][0], chunks_list[i2-1][1]])
            i = i2
    ##    print "length of new_chunks_list: ", len(new_chunks_list)
    ##    print "true coord: "+str(len(true_coords))
    ##    print "miss_itv:", miss_itv
        
    t = d = 0
    trims = []
    TP = FP = TN = FN = 0
    split_itvs = []       # when combining two intervals in the detected group, +1 split
    join_itvs = []        # when combining two intervals in the true group, +1 join
    debug_counter = 0
    debug_dic = {}
    for i in range(len(true_coords)):
        debug_dic[str(true_coords[i])]=" "
    while (t<len(true_coords) and d<len(new_chunks_list)):
        t2 = t
        d2 = d
        if new_chunks_list[d2][1]<true_coords[t2][0]:  # No overlap && detected chunk is ahead of true chunk
            FP = FP + 1
            d = d + 1    # move to the next detected chunk
        elif new_chunks_list[d2][0]>true_coords[t2][1]: # No overlap && detected chunk is behind true chunk
            FN = FN + 1
            t = t + 1    # move to the next true chunk
        #elif not(new_chunks_list[d2-1][1]<true_coords[t2-1][0] or new_chunks_list[d2-1][0]>true_coords[t2-1][1]): # if overlap
        else:                                                                                                      # if overlap
            if d2<len(new_chunks_list)-1:
                while not(new_chunks_list[d2+1][1]<true_coords[t2][0] or new_chunks_list[d2+1][0]>true_coords[t2][1]): # if d2+1 and t2 overlap
                    """ if true:        -------
                        detected:     ---  ---    """
                    split_itvs.append([new_chunks_list[d2][1], new_chunks_list[d2+1][0]])
                    d2 = d2 + 1
            if t2<len(true_coords)-1:
                while not(new_chunks_list[d2][1]<true_coords[t2+1][0] or new_chunks_list[d2][0]>true_coords[t2+1][1]): # if d2 and t2+1 overlap
                    """ if true:        ---  ---
                        detected:     ---------    """
                    join_itvs.append([true_coords[t2][1], true_coords[t2+1][0]])
                    t2 = t2 + 1
            
            if new_chunks_list[d2][1]>=true_coords[t][0] and new_chunks_list[d2][1]<=true_coords[t2][1] and \
               new_chunks_list[d][0]<=true_coords[t][0]:
                """ if true:        -----
                    detected:     ---      """
                if (new_chunks_list[d2][1]-true_coords[t][0])*1.0/(true_coords[t2][1]-true_coords[t][0])>=threshhold:
                    """if D covers more than 0.8T"""
                    TP = TP + 1
                    trims.append([true_coords[t][0]-new_chunks_list[d][0], new_chunks_list[d2][1]-true_coords[t2][1]])
                    d = d2 + 1
                    t = t2 + 1
                else:
                    FP = FP + 1
                    d = d2 + 1
            elif new_chunks_list[d][0]>=true_coords[t][0] and new_chunks_list[d2][1]<=true_coords[t2][1]:
                """ if true:        -----
                    detected:        ---      """
                if (new_chunks_list[d2][1]-new_chunks_list[d][0])*1.0/(true_coords[t2][1]-true_coords[t][0])>=threshhold:
                    """if D covers more than 0.8T """
                    TP = TP + 1
                    trims.append([true_coords[t][0]-new_chunks_list[d][0], new_chunks_list[d2][1]-true_coords[t2][1]])
                    d = d2 + 1
                    t = t2 + 1
                else:
                    FP = FP + 1
                    d = d2 + 1
            elif new_chunks_list[d2][1]>=true_coords[t2][1] and new_chunks_list[d][0]<=true_coords[t][0]:
                """ if true:        -----
                    detected:     ---------
                    In this case D covers all T, hence TP++ """
                TP = TP + 1
                trims.append([true_coords[t][0]-new_chunks_list[d][0], new_chunks_list[d2][1]-true_coords[t2][1]])
                d = d2 + 1
                t = t2 + 1
            elif new_chunks_list[d][0]>=true_coords[t][0] and new_chunks_list[d][0]<=true_coords[t2][1] and \
                 new_chunks_list[d2][1]>true_coords[t2][1]:
                """ if true:      -----
                    detected:        ------      """
                if (true_coords[t2][1]-new_chunks_list[d][0])*1.0/(true_coords[t2][1]-true_coords[t][0])>=threshhold:
                    """if D covers more than 0.8T"""
                    TP = TP + 1
                    trims.append([true_coords[t][0]-new_chunks_list[d][0], new_chunks_list[d2][1]-true_coords[t2][1]])
                    d = d2 + 1
                    t = t2 + 1
                else:
                    FP = FP + 1
                    d = d2 + 1
        if t<len(true_coords) and d==len(new_chunks_list):    # if the pointers has gone through the detected group but not the true group, add to FN
            FN = FN + len(true_coords) - t
        if t==len(true_coords) and d<len(new_chunks_list):    # if the pointers has gone through the true group but not the detected group, add to FP
            FP = FP + len(new_chunks_list) - d
    print "TP =",TP, "FP =",FP, "TN =",TN, "FN =",FN
    print "Sensitivity is:", TP*1.0/((TP+FN)*1.0)
    #print "Trims:", trims
    return trims

def median(pool):
    copy = sorted(pool)
    size = len(copy)
    if size == 0:
        return 0
    elif size % 2 == 1:
        return copy[(size - 1) / 2]
    else:
        return (copy[size/2 - 1] + copy[size/2])*1.0 / 2

def print_trim_criteria(trims):
    avr_sq_trim = left_median = right_median = left_posi_med = left_neg_med = right_posi_med = right_neg_med = 0.0
    if len(trims)!=0:
        avr_sq_trim = sum((start_trim*start_trim+end_trim*end_trim) for [start_trim, end_trim] in trims)*1.0/len(trims)
        #print "average square trim:", avr_sq_trim
        left_median = median([x for [x, y] in trims])
        right_median = median([y for [x, y] in trims])
        left_posi_med = median([x for [x, y] in trims if x > 0])
        left_neg_med = median([x for [x, y] in trims if x < 0])
        right_posi_med = median([y for [x, y] in trims if y > 0])
        right_neg_med = median([y for [x, y] in trims if y < 0])
    print "average square trim:", avr_sq_trim
    print "left median:", left_median
    print "right median:", right_median
    print "left median for positive numbers:", left_posi_med
    print "right median for positive numbers:", right_posi_med
    print "left median for negative numbers:", left_neg_med
    print "right median for negative numbers:", right_neg_med

def no_TE_found(hmmparse):      # If there's no detected TE in hmmparse result file, return True
    fp = open(hmmparse)
    lines = fp.readlines()
    lines[0] = lines[0].rstrip()
    #print lines[0]
    if re.split('\s+', lines[0])[-1] == '0':
        return True
    else:
        return False


if __name__=="__main__":
    proc_id = sys.argv[1]
    if no_TE_found(proc_id + "_hmmparse.out"):
        fp = open(proc_id + "_true_coords.txt")
        lines = fp.readlines()
        print "TP = 0 FP = 0 TN = 0 FN =",len(lines)
        print "Sensitivity is: 0.0"
        print "average square trim: 0.0"
        print "left median: 0.0"
        print "right median: 0.0"
        print "left median for positive numbers: 0.0"
        print "right median for positive numbers: 0.0"
        print "left median for negative numbers: 0.0"
        print "right median for negative numbers: 0.0"
    else:
        result = read_params(proc_id + "_hmmparse.out", proc_id + "_true_coords.txt")
        #print result[0]
        threshhold = 0.8
        trims = result_analyser(result, threshhold)
        print_trim_criteria(trims)
        #print result
