"""Find the Low Speed Region along the human genome.

Cut genome region from N.
/zfssz3/SOLEXA_DATA/BC_RD_P2/USER/huangshujia/ifs5_BC_ONLINE_RD/project/201605_ODPS4WGS/20160626_find_lowspeed_region/bin

Author: Shujia Huang
Date: 2016-06-26
"""
import sys
import optparse

MAXSIZE = 4550000  # 4.6Mb
SLOWCHUNKSIZE = 500000  # 0.5 Mb


def main(opts):
    """main function
    """
    # accbed, accidx = load_accbed(opts.accbed_file)
    fa = get_fasta_size_from_faifile(opts.ref_fa_file)

    # get non N region
    call_region = get_nonN_region(opts.n_reg_file, opts.nsize, fa)
    slow_reg, s_idx = load_slow(opts.slow_region_file)
    for chrom in call_region:
        if chrom[:2] == "GL":
            continue

        for start_pos, end_pos in call_region[chrom]:
            reg = []
            for pos in range(start_pos, end_pos + 1, MAXSIZE):
                s = pos
                e = end_pos if pos + MAXSIZE >= end_pos else pos + MAXSIZE - 1

                flag = True
                ovlp_regs = []
                if chrom in s_idx:
                    for i in range(s_idx[chrom], len(slow_reg[chrom])):

                        if s > slow_reg[chrom][i][1]:
                            continue
                        if e < slow_reg[chrom][i][0]:
                            break

                        if flag:
                            flag = False
                            s_idx[chrom] = i

                        o_s, o_e = overlap_region(s, e, slow_reg[chrom][i][0], slow_reg[chrom][i][1])
                        ovlp_regs.append([o_s, o_e])

                if len(ovlp_regs) > 0:
                    pre_pos = s
                    for start, end in ovlp_regs:

                        if pre_pos < start:
                            reg.append([chrom, pre_pos, start - 1])

                        # reg.append([chrom, start, end])
                        for j in range(start - 1, end, SLOWCHUNKSIZE):
                            tmp_s = j
                            if j + SLOWCHUNKSIZE > end:
                                tmp_e = end
                            else:
                                tmp_e = j + SLOWCHUNKSIZE
                            reg.append([chrom, tmp_s + 1, tmp_e])

                        pre_pos = end + 1

                    if pre_pos <= e:
                        reg.append([chrom, pre_pos, e])

                else:
                    reg.append([chrom, s, e])

            for c, s, e in reg:
                print("%s\t%d\t%d" % (c, s, e))


def get_fasta_size_from_faifile(infile):
    """
    Get chromosome size from .fai file

    .fai format:
        1       249250621       52      60      61
        2       243199373       253404903       60      61
        3       198022430       500657651       60      61
        4       191154276       701980507       60      61
        5       180915260       896320740       60      61
        6       171115067       1080251307      60      61
    """
    fa = {}
    with open(infile) as I:
        for r in I:
            col = r.strip().split()
            fa[col[0]] = int(col[1])

    return fa


def load_accbed(infile):
    region, index = {}, {}
    with open(infile) as I:
        for r in I:
            col = r.strip("\n").split()
            if col[0] not in region:
                region[col[0]] = []
                index[col[0]] = 0

            region[col[0]].append([col[0], int(col[1]), int(col[2])])

    return region, index


def get_nonN_region(infile, thd_n_size, fa):
    """Read N region to get non-n region.
    """
    reg, pos = {}, {}
    pre_id, pre_end, curr_start = "", 0, 0
    with open(infile) as I:
        # dealing with N region file

        for r in I:

            col = r.strip().split()
            n_start, n_end = map(int, (col[1], col[2]))
            n_region_size = n_end - n_start + 1

            if col[0] not in reg:

                reg[col[0]] = []
                pos[col[0]] = 0

                if len(pre_id):
                    reg[pre_id].append([pos[pre_id] + 1, fa[pre_id]])

                pre_id, pre_end = col[0], n_end

            elif n_region_size < thd_n_size:
                # ignore this n region
                pass

            else:
                reg[pre_id].append([pos[pre_id] + 1, n_start - 1])
                pre_id, pre_end = col[0], n_end

            if n_region_size >= thd_n_size:
                pos[col[0]] = n_end

    for chrom, chromsize in fa.items():
        # loading all other non N regions
        if chrom not in reg:
            reg[chrom] = [1, chromsize]

    return reg


def load_slow(infile):
    """
    """
    reg, idx = {}, {}
    with open(infile) as I:
        for line in I:
            # Y  11234914  14674573
            if line[0] == "#":
                continue

            col = line.strip("\n").split()
            if col[0] not in reg:
                reg[col[0]] = []
                idx[col[0]] = 0

            reg[col[0]].append(map(int, [col[1], col[2]]))

    return reg, idx


def overlap_region(start1, end1, start2, end2):
    """
    """
    s, e = 0, 0
    if end1 <= end2 and start1 <= start2:
        s, e = start2, end1

    elif end1 <= end2 and start1 > start2:
        if end1 < start1:
            sys.stderr.write(f"[ERROR] Start1({start1}) > End1({end1})\n")
            sys.exit(1)
        s, e = start1, end1

    elif start1 <= start2 and end1 > end2:
        if end2 < start2:
            sys.stderr.write(f"[ERROR] Start2 ({start2}) > end2 ({end2})\n")
            sys.exit(1)
        s, e = start2, end2

    elif start1 <= end2 < end1:
        s, e = start1, end2

    else:
        sys.stderr.write("[ERROR] Region not right.\n")

    return s, e


if __name__ == "__main__":

    """
    loading all the command options.
    """
    usage = "\nUsage: python %prog [options] -f <fastafile> > Output"
    optp = optparse.OptionParser(usage=usage)
    optp.add_option("-f", "--fa", dest="ref_fa_file", metavar="STR",
                    help="The reference fasta file.", default="")
    # optp.add_option("-b", "--bed", dest="accbed_file", metavar="STR",
    #                help = "The accessible region bed file.", default="")
    optp.add_option("-s", "--slow", dest="slow_region_file", metavar="STR",
                    help="The alow region bed file.", default="")
    optp.add_option("-n", "--nbed", dest="n_reg_file", metavar="STR",
                    help="The N-base region bed file.", default="")
    optp.add_option("-N", "--nbase", dest="nsize", metavar="int",
                    help="The biggest continue N size could be allow "
                         "in one region.", default=1000)
    # optp.add_option("-w", "--win", dest="win", metavar="int",
    #                help = "The window size.", default=100000)

    opt, _ = optp.parse_args()
    if not opt.ref_fa_file:
        optp.error("Required [-f ref_fa_file]\n")

    if not opt.slow_region_file:
        optp.error("Required [-s slow_region_file]\n")
    # if not opt.accbed_file: optp.error("Required [-b accessible bed file.\n]")

    main(opt)
    sys.stderr.write(">> For the flowers bloom in the desert <<")
