hifiasm/Purge_Dups.cpp

#define __STDC_LIMIT_MACROS
#include <stdint.h>
#include <stdlib.h>
#include "ksort.h"
#include "Purge_Dups.h"
#include "Overlaps.h"
#include "Correct.h"
#include "kthread.h"
#include "kdq.h"
#include "hic.h"
#include "partig.h"
#include "rcut.h"

KDQ_INIT(uint64_t)

uint8_t debug_enable = 0;

typedef struct {
    uint64_t weight;
    uint32_t x_beg_pos;
    uint32_t x_end_pos;
    uint32_t y_beg_pos;
    uint32_t y_end_pos;
    uint32_t index_beg;
    uint32_t index_end;
    long long score;
    uint8_t rev;
    asg_arc_t t;
}hap_candidates;

typedef struct {
    kvec_t(hap_candidates) a;
    uint64_t i;
}kvec_hap_candidates;


typedef struct {
    uint64_t* vote_counting;
    uint8_t* visit;
    kvec_t_u64_warp u_vecs;
    kvec_asg_arc_t_offset u_buffer;
    kvec_t_i32_warp u_buffer_tailIndex;
    kvec_t_i32_warp u_buffer_prevIndex;
    kvec_t_u8_warp u_buffer_flag;
    kvec_t_i32_warp u_buffer_beg;
    kvec_hap_candidates u_can;
}hap_alignment_struct;


typedef struct {
    hap_alignment_struct* buf;
    uint32_t num_threads;

    ma_ug_t *ug;
    asg_t *read_g;
    ma_hit_t_alloc* sources;
    ma_hit_t_alloc* reverse_sources;
    R_to_U* ruIndex;
    ma_sub_t *coverage_cut; 
    uint64_t* position_index;
    float Hap_rate;
    int max_hang;
    int min_ovlp;
    float chain_rate;
    hap_overlaps_list* all_ovlp;
    long long cov_threshold;
    hap_cov_t *cov;
}hap_alignment_struct_pip;

typedef struct {
    uint32_t baseBeg, baseEnd;
    uint32_t nodeBeg, nodeEnd;
    uint32_t h_lev_idx;
    uint32_t h_status, c_ug_id;
}p_node_t;

typedef struct {
    uint32_t beg;
    uint32_t occ;
}p_g_in_t;

typedef struct {
    ma_ug_t *ug;
    kvec_t(p_node_t) pg_het_node;
    asg_t *pg_het;
    asg_t *pg_h_lev;
    kvec_t(p_g_in_t) pg_h_lev_idx;
}p_g_t;

void print_peak_line(int c, int x, int exceed, int64_t cnt)
{
	int j;
	if (c >= 0) fprintf(stderr, "[M::%s] %5d: ", __func__, c);
	else fprintf(stderr, "[M::%s] %5s: ", __func__, "rest");
	for (j = 0; j < x; ++j) fputc('*', stderr);
	if (exceed) fputc('>', stderr);
	fprintf(stderr, " %lld\n", (long long)cnt);
}

void print_peak(long long* cov_buf, long long cov_buf_length, long long max_i)
{
    long long i;
    const long long hist_max = 100;
    // print histogram
	for (i = 0; i < cov_buf_length; ++i) 
    {
		long long x, exceed = 0;
		x = (int)((double)hist_max * cov_buf[i] / cov_buf[max_i] + .499);
		if (x > hist_max) exceed = 1, x = hist_max; // may happen if cnt[2] is higher
		if (i > max_i && x == 0) break;
		print_peak_line(i, x, exceed, cov_buf[i]);
	}
	{
		long long x, exceed = 0;
		long long rest = 0;
		for (; i < cov_buf_length; ++i) rest += cov_buf[i];
		x = (int)((double)hist_max * rest / cov_buf[max_i] + .499);
		if (x > hist_max) exceed = 1, x = hist_max;
		print_peak_line(-1, x, exceed, rest);
	}
}


void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo_peak_cov, 
long long* hom_peak, long long* het_peak, long long* k_mer_only, long long* coverage_only)
{
    long long i, start, err_i, max_i, max2_i, max3_i, topo_peak_i, max, max2, max3, topo_peak, min;

    i = start = err_i = max_i = max2_i = max3_i = topo_peak_i = -1;
    max = max2 = max3 = topo_peak = min = -1;

    ///cov_buf[0] is usually very large
    for (i = 1; i < cov_buf_length; ++i)
    {
        if(cov_buf[i] > cov_buf[i-1]) break;
    }
    err_i = i - 1;
    // find the global highest peak
	max_i = err_i + 1, max = cov_buf[max_i];
    for (i = max_i; i < cov_buf_length; ++i)
    {
        if (cov_buf[i] > max)
        {
            max = cov_buf[i];
            max_i = i;
        }
    }

    ///print_peak(cov_buf, cov_buf_length, max_i);

    // look for smaller peak on the low end
	max2 = -1; max2_i = -1;
	for (i = max_i - 1; i > err_i; --i) 
    {
		///at first, it should be a peak
		if (cov_buf[i] >= cov_buf[i-1] && cov_buf[i] >= cov_buf[i+1]) 
        {
			if (cov_buf[i] > max2)
            {
                max2 = cov_buf[i];
                max2_i = i;
            } 
		}
	}

    ///fprintf(stderr, "***max2: %lld, max2_i: %lld\n", max2, max2_i);

    if (max2_i != -1 && max2_i > err_i && max2_i < max_i) 
    {
		for (i = max2_i + 1, min = max; i < max_i; ++i)
        {
            if (cov_buf[i] < min) min = cov_buf[i];
        }
			
		///if the second peak is not significant
		if(max2 < max * 0.05 || min > max2 * 0.95) max2 = max2_i = -1;
	}
		

    // look for smaller peak on the high end
	max3 = -1; max3_i = -1;
    // we'd better use i < cov_buf_length - 1, since cov_buf[cov_buf_length-1] may have problem
	for (i = max_i + 1; i < cov_buf_length - 1; ++i) 
    {
		//at first, it should be a peak
		if (cov_buf[i] >= cov_buf[i-1] && cov_buf[i] >= cov_buf[i+1]) 
        {
			if (cov_buf[i] > max3)
            {
                max3 = cov_buf[i], max3_i = i;
            } 
		}
	}

    ///fprintf(stderr, "***max3: %lld, max3_i: %lld\n", max3, max3_i);
    
    //if found a peak
	if (max3 != -1 && max3_i > max_i) 
    {
		for (i = max_i + 1, min = max; i < max3_i; ++i)
        {
            if (cov_buf[i] < min) min = cov_buf[i];
        }
			
		if (max3 < max * 0.05 || min > max3 * 0.95 || max3_i > max_i * 3) max3 = max3_i = -1;
	}


    

    (*hom_peak) = (*het_peak) = -1;
    if (topo_peak_cov && (*topo_peak_cov) < cov_buf_length)
    {
        topo_peak_i = (*topo_peak_cov);
        topo_peak = cov_buf[topo_peak_i];
        if (topo_peak <= max * 0.05) topo_peak_i = topo_peak = -1;
    }
    
    long long k_mer_het, k_mer_hom, coverage_het, coverage_hom, alter_peak;
    k_mer_het = k_mer_hom = coverage_het = coverage_hom = alter_peak = -1;

    alter_peak = topo_peak_i;
    k_mer_het = asm_opt.het_cov;
    k_mer_hom = asm_opt.hom_cov;
    if(max3_i > 0)
    {
        coverage_het = max_i;
        coverage_hom = max3_i;
    }
    else
    {
        coverage_het = max2_i;
        coverage_hom = max_i;
    }

    if(k_mer_het != -1)
    {
        (*het_peak) = k_mer_het;
        (*hom_peak) = k_mer_hom;
        return;
    }
    else if(coverage_het != -1)
    {
        (*het_peak) = coverage_het;
        (*hom_peak) = coverage_hom;
        return;
    }
    else if(k_mer_hom > coverage_hom*1.5)
    {
        (*het_peak) = coverage_hom;
        (*hom_peak) = k_mer_hom;
        return;
    }
    else if(alter_peak != -1)
    {
        ///if peak is het, coverage peak is more reliable 
        if(coverage_hom >= alter_peak*0.8 && coverage_hom <= alter_peak*1.2) 
        {
            (*het_peak) = coverage_hom;
            return;
        }///if peak is homo, k-mer peak is more reliable 
        else if(k_mer_hom >= alter_peak*0.8*2 && k_mer_hom <= alter_peak*1.2*2)
        {
            (*hom_peak) = k_mer_hom;
            return;
        }
    }

    (*k_mer_only) = k_mer_hom;
    (*coverage_only) = coverage_hom;
    

    // fprintf(stderr, "max: %lld, max_i: %lld\n", max, max_i);
    // fprintf(stderr, "max2: %lld, max2_i: %lld\n", max2, max2_i);
    // fprintf(stderr, "max3: %lld, max3_i: %lld\n", max3, max3_i);
    // fprintf(stderr, "[M::%s] Heterozygous k-mer peak: %d\n", __func__, asm_opt.het_cov);
    // fprintf(stderr, "[M::%s] Homozygous k-mer peak: %d\n", __func__, asm_opt.hom_cov);
    // fprintf(stderr, "[M::%s] Heterozygous coverage peak: %lld\n", __func__, (*het_peak));
    // fprintf(stderr, "[M::%s] Homozygous coverage peak: %lld\n", __func__, (*hom_peak));
    // fprintf(stderr, "[M::%s] Alter coverage peak: %lld\n", __func__, topo_peak_i);
}




long long get_alter_peak(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, uint64_t* position_index, 
ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length)
{
    
    ma_utg_t* u = NULL;
    asg_t* nsg = ug->g;
    uint64_t v, j, k, qn, n_vtx = nsg->n_seq, primary_bases = 0, alter_bases = 0;
    uint32_t tn, is_Unitig;
    long long* cov_buf = NULL;
    ma_hit_t *h;
    cov_buf = (long long*)calloc(cov_buf_length, sizeof(long long));
    long long R_bases = 0, C_bases_primary = 0, C_bases_alter = 0, C_bases = 0;
    memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq);


    for (v = 0; v < n_vtx; ++v) 
    {
        if(nsg->seq[v].del) continue;
        if(nsg->seq[v].c == ALTER_LABLE) continue;
        u = &(ug->u.a[v]);
        if(u->m == 0) continue;
        for (k = 0; k < u->n; k++)
        {
            qn = u->a[k]>>33;
            position_index[qn] = 0;
            R_bases = coverage_cut[qn].e - coverage_cut[qn].s;
            primary_bases += R_bases;
        }
    }

    for (qn = 0; qn < read_g->n_seq; qn++)
    {
        if(position_index[qn] == 0) continue;
        if(read_g->seq[qn].del) continue;

        C_bases = C_bases_primary = C_bases_alter = 0;
        R_bases = coverage_cut[qn].e - coverage_cut[qn].s;
        alter_bases += R_bases;
        for (j = 0; j < (uint64_t)(sources[qn].length); j++)
        {
            h = &(sources[qn].buffer[j]);
            if(h->el != 1) continue;
            tn = Get_tn((*h));

            if(read_g->seq[tn].del == 1)
            {
                ///get the id of read that contains it 
                get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
                if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue;
            }

            if(position_index[tn] == 0)
            {
                C_bases_primary += Get_qe((*h)) - Get_qs((*h));
            }
            else
            {
                C_bases_alter += Get_qe((*h)) - Get_qs((*h));
            }
        }

        // if(qn == 1893151 || qn == 1929038)
        // {
        //     fprintf(stderr, "qn: %lu, C_bases_primary: %lld, C_bases_alter: %lld, C_bases: %lld\n",
        //     qn, C_bases_primary, C_bases_alter, C_bases);
        // }

        C_bases = C_bases_primary + C_bases_alter;
        if(C_bases_alter < C_bases * ALTER_COV_THRES) continue;

        C_bases = C_bases/R_bases;
        if(C_bases < 0 || C_bases >= cov_buf_length) continue;
        cov_buf[C_bases]++;
    }
    
    long long max_i = -1, max = -1;
    for (j = 0; (long long)j < cov_buf_length; ++j)
    {
        if (cov_buf[j] > max)
        {
            max = cov_buf[j];
            max_i = j;
        }
    }

    ///fprintf(stderr, "alter max_i: %lld, max: %lld\n", max_i, max);
    ///if(max_i < 5) max_i = max = -1;
    if(alter_bases < primary_bases * REAL_ALTER_THRES) max_i = max = -1;

    free(cov_buf);
    memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq);

    return max_i;
}

long long get_read_coverage_thres(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, uint64_t* position_index, 
ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, uint64_t n_read, long long cov_buf_length,
long long* k_mer_only, long long* coverage_only)
{
    uint64_t i, j;
    long long* cov_buf = NULL;
    ma_hit_t *h;
    cov_buf = (long long*)calloc(cov_buf_length, sizeof(long long));
    long long R_bases = 0, C_bases = 0;
    for (i = 0; i < n_read; ++i) 
    {
        C_bases = 0;
        R_bases = coverage_cut[i].e - coverage_cut[i].s;
        for (j = 0; j < (uint64_t)(sources[i].length); j++)
        {
            h = &(sources[i].buffer[j]);
            if(h->el != 1) continue;
            C_bases += Get_qe((*h)) - Get_qs((*h));
        }
        C_bases = C_bases/R_bases;
        if(C_bases < 0 || C_bases >= cov_buf_length) continue;
        cov_buf[C_bases]++;
    }

    long long alter_peak = -1, hom_peak = -1, het_peak = -1;
    if(position_index)
    {
        alter_peak = get_alter_peak(ug, read_g, ruIndex, position_index, sources, coverage_cut, 
        cov_buf_length);
    }
    
    get_read_peak(cov_buf, cov_buf_length, alter_peak == -1? NULL: &alter_peak, &hom_peak, &het_peak,
    k_mer_only, coverage_only);

    free(cov_buf);

    if(hom_peak != -1) return hom_peak*HOM_PEAK_RATE;
    if(het_peak != -1) return het_peak*HET_PEAK_RATE;
    return -1;
}










void init_hap_alignment_struct(hap_alignment_struct* x, uint32_t size)
{
    x->vote_counting = (uint64_t*)malloc(sizeof(uint64_t)*size);
    memset(x->vote_counting, 0, sizeof(uint64_t)*size);

    x->visit = (uint8_t*)malloc(sizeof(uint8_t)*size);
    memset(x->visit, 0, size);

    kv_init(x->u_vecs.a);
    kv_init(x->u_buffer.a);
    kv_init(x->u_buffer_tailIndex.a);
    kv_init(x->u_buffer_prevIndex.a);
    kv_init(x->u_buffer_beg.a);
    kv_init(x->u_buffer_flag.a);
    kv_init(x->u_can.a);
}

void destory_hap_alignment_struct(hap_alignment_struct* x)
{
    free(x->vote_counting);
    free(x->visit);
    kv_destroy(x->u_vecs.a);
    kv_destroy(x->u_buffer.a);
    kv_destroy(x->u_buffer_tailIndex.a);
    kv_destroy(x->u_buffer_prevIndex.a);
    kv_destroy(x->u_buffer_beg.a);
    kv_destroy(x->u_buffer_flag.a);
    kv_destroy(x->u_can.a);
}

void init_hap_alignment_struct_pip(hap_alignment_struct_pip* x, uint32_t num_threads, uint32_t n_seq,
ma_ug_t *ug, asg_t *read_g,  ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, ma_sub_t *coverage_cut, 
uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chain_rate, hap_overlaps_list* all_ovlp, hap_cov_t *cov)
{
    uint32_t i;
    x->num_threads = num_threads;
    x->buf = (hap_alignment_struct*)malloc(sizeof(hap_alignment_struct)*x->num_threads);
    for (i = 0; i < x->num_threads; i++)
    {
        init_hap_alignment_struct(&(x->buf[i]), n_seq);
    }

    x->ug = ug;
    x->read_g = read_g;
    x->sources = sources;
    x->reverse_sources = reverse_sources;
    x->ruIndex = ruIndex;
    x->coverage_cut = coverage_cut;
    x->position_index = position_index;
    x->Hap_rate = Hap_rate;
    x->max_hang = max_hang;
    x->min_ovlp = min_ovlp;
    x->chain_rate = chain_rate;
    x->all_ovlp = all_ovlp;
    x->cov = cov;
}


void destory_hap_alignment_struct_pip(hap_alignment_struct_pip* x)
{
    uint32_t i;
    for (i = 0; i < x->num_threads; i++)
    {
        destory_hap_alignment_struct(&(x->buf[i]));
    }

    free(x->buf);
}

void init_hap_overlaps_list(hap_overlaps_list* x, uint32_t num)
{
    uint32_t i = 0;
    x->num = num;
    x->x = (kvec_hap_overlaps*)malloc(sizeof(kvec_hap_overlaps)*x->num);
    for (i = 0; i < x->num; i++)
    {
        kv_init(x->x[i].a);
    }
}

void enable_debug_mode(uint32_t mode)
{
    debug_enable = mode;
}

void destory_hap_overlaps_list(hap_overlaps_list* x)
{
    uint32_t i = 0;
    for (i = 0; i < x->num; i++)
    {
        kv_destroy(x->x[i].a);
    }
    free(x->x);
}



inline void clean_visit_flag(uint8_t* visit, asg_t *read_g, R_to_U* ruIndex, uint32_t contigNum, 
ma_hit_t_alloc* x)
{
    uint32_t k, rId, is_Unitig, Hap_cId;

    if(x->length*2 > contigNum)
    {
        memset(visit, 0, contigNum);
    }
    else
    {
        for (k = 0; k < x->length; k++)
        {
            rId = Get_tn(x->buffer[k]);
        
            if(read_g->seq[rId].del == 1)
            {
                ///get the id of read that contains it 
                get_R_to_U(ruIndex, rId, &rId, &is_Unitig);
                if(rId == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[rId].del == 1) continue;
            }

            ///there are two cases: 
            ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id  
            ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1
            get_R_to_U(ruIndex, rId, &Hap_cId, &is_Unitig);
            if(is_Unitig == 0 || Hap_cId == (uint32_t)-1) continue;
            ///here rId is the id of the read coming from the different haplotype
            ///Hap_cId is the id of the corresponding contig (note here is the contig, instead of untig)
            visit[Hap_cId] = 0;
        }
    }
}

uint32_t prefilter(uint32_t x_pos, uint32_t y_pos, uint32_t xLen, uint32_t yLen, uint32_t dir, 
float Hap_rate, uint32_t seedOcc)
{
    uint32_t max_count = 0, min_count = 0;
    uint32_t /**xLeftBeg, **/xLeftLen, yLeftBeg, yLeftLen;
    uint32_t xRightBeg, xRightLen, yRightBeg, yRightLen;
    if(dir == 0)
    {
        /**xLeftBeg = 0;**/ xLeftLen = x_pos; xRightBeg = x_pos; xRightLen = xLen - xRightBeg;
        yLeftBeg = 0; yLeftLen = y_pos; yRightBeg = y_pos; yRightLen = yLen - yRightBeg;
    }
    else
    {
        /**xLeftBeg = 0;**/ xLeftLen = x_pos; xRightBeg = x_pos; xRightLen = xLen - xRightBeg;

        yLeftBeg = y_pos + 1; yLeftLen = yLen - yLeftBeg;
        yRightBeg = 0; yRightLen = y_pos + 1;
    }



    max_count = seedOcc;
    min_count = MIN(xLeftLen, yLeftLen) + MIN(xRightLen, yRightLen);
    if(min_count == 0) return NON_PLOID;
    if(max_count <= min_count*Hap_rate) return NON_PLOID;
    return PLOID;
}



inline uint64_t get_xy_pos(asg_t *read_g, asg_arc_t* t, uint32_t v_in_unitig, uint32_t w_in_unitig, 
uint32_t xUnitigLen, uint32_t yUnitigLen, uint64_t* position_index, uint8_t* rev)
{
    uint32_t x_pos, y_pos, x_dir = 0, y_dir = 0;
    uint64_t tmp;
    x_pos = y_pos = (uint32_t)-1;
    if((t->ul>>32)==v_in_unitig)///end pos
    {
        x_pos = (position_index[v_in_unitig>>1]>>32) + read_g->seq[v_in_unitig>>1].len - 1;
        x_dir = 0;
    }
    else if((t->ul>>32)==(v_in_unitig^1))///start pos
    {
        x_pos = (position_index[v_in_unitig>>1]>>32);
        x_dir = 1;
    }
    else
    {
        fprintf(stderr, "ERROR\n");
    }
    
    if(t->v == w_in_unitig)
    {
        y_pos = (position_index[w_in_unitig>>1]>>32) + t->ol - 1;
        y_dir = 0;
    }
    else if(t->v == (w_in_unitig^1))
    {
        y_pos = (position_index[w_in_unitig>>1]>>32) + read_g->seq[w_in_unitig>>1].len - t->ol;
        y_dir = 1;
    }
    else
    {
        fprintf(stderr, "ERROR\n");
    }

    (*rev) = x_dir^y_dir;
    if((*rev))
    {
        if(yUnitigLen <= y_pos)
        {
            y_pos = (uint32_t)-1;
        }
        else
        {
            y_pos = yUnitigLen - y_pos - 1;
        }
    } 

    if(x_pos>=xUnitigLen) x_pos = (uint32_t)-1;
    if(y_pos>=yUnitigLen) y_pos = (uint32_t)-1;
    
    tmp = x_pos; tmp = tmp << 32; tmp = tmp | y_pos;
    return tmp;
}


void print_debug_unitig(ma_utg_t *xReads, uint64_t* position_index, const char* infor)
{
    uint32_t k;
    fprintf(stderr, "\n%s: n = %u\n", infor, xReads->n);
    for (k = 0; k < xReads->n; k++)
    {
        fprintf(stderr, "(%u)v: %u, len: %u, index: %u, pos: %u\n", 
        k, (uint32_t)(xReads->a[k]>>32), (uint32_t)xReads->a[k], (uint32_t)(position_index[xReads->a[k]>>33]),
        (uint32_t)(position_index[xReads->a[k]>>33]>>32));
    }
}

void deduplicate_edge(kvec_asg_arc_t_offset* u_buffer)
{
    if(u_buffer->a.n == 0) return;
    long long i = u_buffer->a.n - 1, k, i_off;
    uint32_t v = u_buffer->a.a[i].x.ul>>33, m;

    for (; i >= 0; i--)
    {
        if((u_buffer->a.a[i].x.ul>>33) != v)
        {
            break;
        } 
    }

    ///fprintf(stderr, "u_buffer->a.n: %u, i: %lld\n", u_buffer->a.n, i);

    i = i + 1;
    for (m = i; i < (long long)u_buffer->a.n; i++)
    {
        if(u_buffer->a.a[i].x.del) continue;

        i_off = Cal_Off(u_buffer->a.a[i].Off);

        for (k = i + 1; k < (long long)u_buffer->a.n; k++)
        {
            if(u_buffer->a.a[k].x.del) continue;
            if(u_buffer->a.a[i].x.el != u_buffer->a.a[k].x.el) continue;
            if(i_off != Cal_Off(u_buffer->a.a[k].Off)) continue;
            u_buffer->a.a[k].x.del = 1;
            u_buffer->a.a[i].weight += u_buffer->a.a[k].weight;
        }
        u_buffer->a.a[m] = u_buffer->a.a[i];
        m++;
    }

    u_buffer->a.n = m;
    ///fprintf(stderr, "u_buffer->a.n: %u, i: %lld\n", u_buffer->a.n, i);
}

int cmp_hap_alignment(const void * a, const void * b)
{
    if((*(asg_arc_t_offset*)a).x.el > (*(asg_arc_t_offset*)b).x.el) return 1;
    if((*(asg_arc_t_offset*)a).x.el < (*(asg_arc_t_offset*)b).x.el) return -1;

    long long aOff = Cal_Off((*(asg_arc_t_offset*)a).Off);
    long long bOff = Cal_Off((*(asg_arc_t_offset*)b).Off);

    if(aOff > bOff) return 1;
    if(aOff < bOff) return -1;

    if(((*(asg_arc_t_offset*)a).Off>>32) > ((*(asg_arc_t_offset*)b).Off>>32)) return 1;
    if(((*(asg_arc_t_offset*)a).Off>>32) < ((*(asg_arc_t_offset*)b).Off>>32)) return -1;

    if((uint32_t)((*(asg_arc_t_offset*)a).Off) > (uint32_t)((*(asg_arc_t_offset*)b).Off)) return 1;
    if((uint32_t)((*(asg_arc_t_offset*)a).Off) < (uint32_t)((*(asg_arc_t_offset*)b).Off)) return -1;

    if((*(asg_arc_t_offset*)a).weight < (*(asg_arc_t_offset*)b).weight) return 1;
    if((*(asg_arc_t_offset*)a).weight > (*(asg_arc_t_offset*)b).weight) return -1;

    return 0;
}


int cmp_hap_alignment_chaining(const void * a, const void * b)
{
    if((*(asg_arc_t_offset*)a).x.el > (*(asg_arc_t_offset*)b).x.el) return 1;
    if((*(asg_arc_t_offset*)a).x.el < (*(asg_arc_t_offset*)b).x.el) return -1;

    if(((*(asg_arc_t_offset*)a).Off>>32) > ((*(asg_arc_t_offset*)b).Off>>32)) return 1;
    if(((*(asg_arc_t_offset*)a).Off>>32) < ((*(asg_arc_t_offset*)b).Off>>32)) return -1;

    if((uint32_t)((*(asg_arc_t_offset*)a).Off) > (uint32_t)((*(asg_arc_t_offset*)b).Off)) return 1;
    if((uint32_t)((*(asg_arc_t_offset*)a).Off) < (uint32_t)((*(asg_arc_t_offset*)b).Off)) return -1;

    return 0;
}

int cmp_hap_candidates(const void * a, const void * b)
{
    if((*(hap_candidates*)a).weight < (*(hap_candidates*)b).weight) return 1;
    if((*(hap_candidates*)a).weight > (*(hap_candidates*)b).weight) return -1;

    if((*(hap_candidates*)a).index_beg > (*(hap_candidates*)b).index_beg) return 1;
    if((*(hap_candidates*)a).index_beg < (*(hap_candidates*)b).index_beg) return -1;
    
    return 0;
}
inline long long get_hap_overlapLen(long long x_beg, long long x_end, long long xLen, 
long long y_beg, long long y_end, long long yLen, long long* n_x_beg, long long* n_x_end,
long long* n_y_beg, long long* n_y_end)
{
    if(x_beg <= y_beg)
    {
        y_beg = y_beg - x_beg;
        x_beg = 0;
    }
    else
    {
        x_beg = x_beg - y_beg;
        y_beg = 0;
    }

    long long x_right_length = xLen - x_end - 1;
    long long y_right_length = yLen - y_end - 1;


    if(x_right_length <= y_right_length)
    {
        x_end = xLen - 1;
        y_end = y_end + x_right_length;        
    }
    else
    {
        x_end = x_end + y_right_length;
        y_end = yLen - 1;
    }

    if(n_x_beg) (*n_x_beg) = x_beg;
    if(n_x_end) (*n_x_end) = x_end;
    if(n_y_beg) (*n_y_beg) = y_beg;
    if(n_y_end) (*n_y_end) = y_end;

    return x_end - x_beg + 1;
}



uint32_t classify_hap_overlap(long long xBeg, long long xEnd, long long xLen,
long long yBeg, long long yEnd, long long yLen, long long* r_xBeg, long long* r_xEnd, 
long long* r_yBeg, long long* r_yEnd)
{
    long long n_x_beg, n_x_end, n_y_beg, n_y_end;
    get_hap_overlapLen(xBeg, xEnd, xLen, yBeg, yEnd, yLen, &n_x_beg, &n_x_end, &n_y_beg, &n_y_end);
    if(r_xBeg) (*r_xBeg) = n_x_beg;
    if(r_xEnd) (*r_xEnd) = n_x_end;
    if(r_yBeg) (*r_yBeg) = n_y_beg;
    if(r_yEnd) (*r_yEnd) = n_y_end;
    if(n_x_beg == 0 && n_x_end == xLen - 1) return YCX;
    if(n_y_beg == 0 && n_y_end == yLen - 1) return XCY;
    if(n_y_beg == 0 && n_x_end == xLen - 1) return X2Y;
    if(n_x_beg == 0 && n_y_end == yLen - 1) return Y2X;
    return XCY;
}


uint64_t get_pair_hap_coverage(uint64_t* readIDs, uint32_t Len, ma_hit_t_alloc* sources, ma_sub_t* coverage_cut)
{
    uint32_t m, n, qn;
    ma_hit_t *h;
    uint64_t R_bases = 0, C_bases = 0;

    for (m = 0; m < Len; m++)
    {
        qn = readIDs[m]>>33;
        R_bases += coverage_cut[qn].e - coverage_cut[qn].s;
        for (n = 0; n < (uint64_t)(sources[qn].length); n++)
        {
            h = &(sources[qn].buffer[n]);
            C_bases += Get_qe((*h)) - Get_qs((*h));
        }
    }

    return C_bases/R_bases;
}


uint64_t get_pair_purge_coverage(ma_utg_t *xReads, long long xPosBeg, long long xPosEnd, 
ma_utg_t *yReads, long long yPosBeg, long long yPosEnd, uint32_t rev, asg_t *read_g, hap_cov_t *cov)
{
    long long offset, r_beg, r_end, i_beg, i_end, ovlp, IdxBeg, IdxEnd;
    uint64_t i, rId, uCov, uLen;
    ma_utg_t *x = NULL;
    uCov = uLen = 0;
    if(rev) 
    {
        yPosBeg = yReads->len - yPosBeg - 1;
        yPosEnd = yReads->len - yPosEnd - 1;
        offset = yPosBeg; yPosBeg = yPosEnd; yPosEnd = offset;
    }
    

    IdxBeg = IdxEnd = -1; 
    x = xReads; i_beg = xPosBeg; i_end = xPosEnd;
    for (i = 0, offset = 0; i < x->n; i++)
    {
        rId = x->a[i]>>33;
        r_beg = offset; r_end = offset + (long long)(read_g->seq[rId].len) - 1;
        offset += (uint32_t)x->a[i];
        ovlp = (long long)(MIN(r_end, i_end)) - (long long)(MAX(r_beg, i_beg)) + 1;
        if(ovlp <= 0 || ovlp < read_g->seq[rId].len * 0.8)
        {
            if(IdxBeg != -1 && IdxEnd != -1) break;
            continue;
        }

        if(IdxBeg == -1) IdxBeg = i;
        IdxEnd = i;
    }
    if(IdxBeg != -1 && IdxEnd != -1)
    {
        for (i = IdxBeg; (long long)i <= IdxEnd; i++)
        {
            rId = x->a[i]>>33;
            uCov += cov->cov[rId];
            uLen += cov->read_g->seq[rId].len;
        }
    }    



    IdxBeg = IdxEnd = -1; 
    x = yReads; i_beg = yPosBeg; i_end = yPosEnd;
    for (i = 0, offset = 0; i < x->n; i++)
    {
        rId = x->a[i]>>33;
        r_beg = offset; r_end = offset + (long long)(read_g->seq[rId].len) - 1;
        offset += (uint32_t)x->a[i];
        ovlp = (long long)(MIN(r_end, i_end)) - (long long)(MAX(r_beg, i_beg)) + 1;
        if(ovlp <= 0 || ovlp < read_g->seq[rId].len * 0.8)
        {
            if(IdxBeg != -1 && IdxEnd != -1) break;
            continue;
        }

        if(IdxBeg == -1) IdxBeg = i;
        IdxEnd = i;
    }
    if(IdxBeg != -1 && IdxEnd != -1)
    {
        for (i = IdxBeg; (long long)i <= IdxEnd; i++)
        {
            rId = x->a[i]>>33;
            uCov += cov->cov[rId];
            uLen += cov->read_g->seq[rId].len;
        }
    }

    return (uLen == 0? 0 : uCov / uLen);
}


void get_pair_hap_similarity_by_base(ma_utg_t *xReads, asg_t *read_g, uint32_t target_uId, 
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, long long xBegPos, long long xEndPos, 
double* Match, double* Total)
{
    uint32_t i, j, qn, tn, is_Unitig, uId, min_count = 0, max_count = 0;
    long long offset, r_beg, r_end, ovlp;
    
    for (i = 0, offset = 0; i < xReads->n; i++)
    {
        qn = xReads->a[i]>>33;
        r_beg = offset; r_end = offset + (long long)(read_g->seq[qn].len) - 1;
        offset += (uint32_t)xReads->a[i];
        
        ovlp = (long long)(MIN(r_end, xEndPos)) - (long long)(MAX(r_beg, xBegPos)) + 1;
        if(ovlp <= 0) continue;

        if(reverse_sources[qn].length > 0) min_count++;
        if(reverse_sources[qn].length == 0) continue;
        for (j = 0; j < reverse_sources[qn].length; j++)
        {
            tn = Get_tn(reverse_sources[qn].buffer[j]);
            if(read_g->seq[tn].del == 1)
            {
                get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
                if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue;
            }


            get_R_to_U(ruIndex, tn, &uId, &is_Unitig);
            if(uId!=(uint32_t)-1 && is_Unitig == 1 && uId == target_uId)
            {
                max_count++;
                break;
            }
        }
    }

    (*Match) = max_count;
    (*Total) = min_count;
}

void get_pair_hap_similarity(uint64_t* readIDs, uint32_t Len, uint32_t target_uId, 
ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match, double* Total)
{
    #define CUTOFF_THRES 1000
    uint32_t i, j, qn, tn, is_Unitig, uId, min_count = 0, max_count = 0, cutoff = 0;;
    for (i = 0; i < Len; i++)
    {
        if(cutoff > CUTOFF_THRES)
        {
            max_count = 0;
            min_count = Len;
            break;
        }
        qn = readIDs[i]>>33;
        if(reverse_sources[qn].length > 0) min_count++;
        if(reverse_sources[qn].length == 0) continue;
        for (j = 0; j < reverse_sources[qn].length; j++)
        {
            tn = Get_tn(reverse_sources[qn].buffer[j]);
            if(read_g->seq[tn].del == 1)
            {
                get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
                if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue;
            }


            get_R_to_U(ruIndex, tn, &uId, &is_Unitig);
            if(uId!=(uint32_t)-1 && is_Unitig == 1 && uId == target_uId)
            {
                max_count++;
                break;
            }
        }

        //means no match
        if(j == reverse_sources[qn].length)
        {
            cutoff++;
        }
        else
        {
            cutoff = 0;
        }
    }

    (*Match) = max_count;
    (*Total) = min_count;
}
/**
void get_pair_hap_similarity_deduplicate(uint64_t* readIDs, uint32_t Len, uint32_t target_uId, 
ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match, double* Total)
{
    get_pair_hap_similarity(readIDs, Len, target_uId, reverse_sources, read_g, ruIndex, Match, Total);
    return;

    #define CUTOFF_THRES 100
    uint32_t i, j, qn, tn, is_Unitig, uId, min_count = 0, max_count = 0, cutoff = 0, is_found;
    for (i = 0; i < Len; i++)
    {
        if(cutoff > CUTOFF_THRES)
        {
            max_count = 0;
            min_count = Len;
            break;
        }
        qn = readIDs[i]>>33;
        is_found = 0;

        for (j = 0; j < reverse_sources[qn].length; j++)
        {
            tn = Get_tn(reverse_sources[qn].buffer[j]);
            if(read_g->seq[tn].del == 1)
            {
                get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
                if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue;
            }


            get_R_to_U(ruIndex, tn, &uId, &is_Unitig);
            if(uId!=(uint32_t)-1 && is_Unitig == 1 && uId == target_uId)
            {
                max_count++;
            }
            min_count++;
            is_found = 1;
        }

        //means there is a match
        if(is_found)
        {
            cutoff = 0;
        }
        else
        {
            cutoff++;
        }
    }

    (*Match) = max_count;
    (*Total) = min_count;
}
**/

inline void check_hap_match(uint32_t qn, uint32_t targetBeg, uint32_t targetEnd, uint32_t targetID, 
uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, uint32_t* is_found, uint32_t* is_match)
{
    uint32_t j, tn, uId, is_Unitig, offset;
    (*is_found) = (*is_match) = 0;
    if(reverse_sources[qn].length > 0) (*is_found) = 1;
    for (j = 0; j < reverse_sources[qn].length; j++)
    {
        tn = Get_tn(reverse_sources[qn].buffer[j]);
        if(read_g->seq[tn].del == 1)
        {
            get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
            if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue;
        }


        get_R_to_U(ruIndex, tn, &uId, &is_Unitig);
        if(uId!=(uint32_t)-1 && is_Unitig == 1 && uId == targetID)
        {
            offset = (uint32_t)(position_index[tn]);
            if(offset >= targetBeg && offset <= targetEnd)
            {
                (*is_match) = 1;
                break;
            }                    
        }
    }

}

/**
inline void check_hap_match_deduplicate(uint32_t qn, uint32_t targetBeg, uint32_t targetEnd, uint32_t targetID, 
uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, uint32_t* is_found, uint32_t* is_match)
{
    check_hap_match(qn, targetBeg, targetEnd, targetID, position_index, reverse_sources, read_g, 
    ruIndex, is_found, is_match);
    return;


    uint32_t j, tn, uId, is_Unitig, offset;
    (*is_found) = (*is_match) = 0;

    for (j = 0; j < reverse_sources[qn].length; j++)
    {
        tn = Get_tn(reverse_sources[qn].buffer[j]);
        if(read_g->seq[tn].del == 1)
        {
            get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
            if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue;
        }


        get_R_to_U(ruIndex, tn, &uId, &is_Unitig);
        if(uId!=(uint32_t)-1 && is_Unitig == 1 && uId == targetID)
        {
            offset = (uint32_t)(position_index[tn]);
            if(offset >= targetBeg && offset <= targetEnd)
            {
                (*is_match)++;
            }                    
        }
        (*is_found)++;
    }

}
**/

void determin_hap_alignment_boundary_single_side(uint64_t* readIDs, long long queryLen, long long targetBeg, 
long long targetEnd, long long targetID, long long eMatch, long long eTotal, long long dir, 
float H_rate, int is_local, uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, 
R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_count)
{
    if(queryLen == 0)
    {
        (*n_matchLen) = (*n_min_count) = (*n_max_count) = 0;
        return;
    }
    long long i, maxId, min_count = eTotal, max_count = eMatch, matchLen = 0;
    long long rLen, score = 0, max_score = 0;
    uint32_t is_found, is_match;
    if(dir == 0)
    {
        for (i = 0, maxId = 0; i < queryLen; i++)
        {
            
            check_hap_match(readIDs[i]>>33, targetBeg, targetEnd, targetID, position_index, reverse_sources, 
            read_g, ruIndex, &is_found, &is_match);

            min_count += is_found;
            max_count += is_match;
            if(max_count > min_count*H_rate) maxId = i;

            if(is_local && is_found)
            {
                rLen = read_g->seq[readIDs[i]>>33].len;
                score += (is_match? rLen : (rLen*(-1)));
                if(score >= max_score) max_score = score, maxId = i;
            }
        }

        for (i = maxId; i >= 0; i--)
        {
            check_hap_match(readIDs[i]>>33, targetBeg, targetEnd, targetID, position_index, reverse_sources, 
            read_g, ruIndex, &is_found, &is_match);

            ///if(is_found > 0 && is_match > 0 && is_match > is_found*Hap_rate)
            if(is_found == 1 && is_match == 1)
            {
                break;
            }
            min_count -= is_found;
            max_count -= is_match;
        }
        
        matchLen = i+1;
    }
    else
    {
        for (i = queryLen - 1, maxId = queryLen - 1; i >= 0; i--)
        {
            check_hap_match(readIDs[i]>>33, targetBeg, targetEnd, targetID, position_index, reverse_sources, 
            read_g, ruIndex, &is_found, &is_match);

            min_count += is_found;
            max_count += is_match;
            if(max_count > min_count*H_rate) maxId = i;

            if(is_local && is_found)
            {
                rLen = read_g->seq[readIDs[i]>>33].len;
                score += (is_match? rLen : (rLen*(-1)));
                if(score >= max_score) max_score = score, maxId = i;
            }
        }

        for (i = maxId; i < queryLen; i++)
        {
            check_hap_match(readIDs[i]>>33, targetBeg, targetEnd, targetID, position_index, reverse_sources, 
            read_g, ruIndex, &is_found, &is_match);

            ///if(is_found > 0 && is_match > 0 && is_match > is_found*Hap_rate)
            if(is_found == 1 && is_match == 1)
            {
                break;
            }
            min_count -= is_found;
            max_count -= is_match;
        }

        matchLen = queryLen - i;
    }

    ///need to check if min_count == 0
    if(min_count == 0)
    {
        (*n_matchLen) = (*n_min_count) = (*n_max_count) = 0;
        return;
    }

    (*n_matchLen) = matchLen;
    (*n_min_count) = min_count;
    (*n_max_count) = max_count;
}

inline void modify_target_interval(long long beg, long long end, long long len,
long long* target_beg, long long* target_end)
{
    #define TARGET_SGIFT 3
    beg -= TARGET_SGIFT;
    end += TARGET_SGIFT;
    if(beg < 0) beg = 0;
    if(end >= len) end = len - 1;
    (*target_beg) = beg;
    (*target_end) = end;
}


void bi_direction_hap_alignment_extention(ma_utg_t* xReads, uint32_t xLeftBeg, uint32_t xLeftLen,
uint32_t xRightBeg, uint32_t xRightLen, uint32_t targetUid, uint32_t target_beg, uint32_t target_end,
float Hap_rate, int is_local, uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex,
uint32_t rev, long long* x_interval_beg, long long* x_interval_end)
{
    if(rev)
    {
        uint32_t k;
        k = xLeftBeg; xLeftBeg = xRightBeg; xRightBeg = k;
        k = xLeftLen; xLeftLen = xRightLen; xRightLen = k;
    }
    uint32_t n_matchLenLeft, x_max_countLeft, x_min_countLeft;
    uint32_t n_matchLenRight, x_max_countRight, x_min_countRight;
    n_matchLenLeft = x_max_countLeft = x_min_countLeft = 0;
    determin_hap_alignment_boundary_single_side(xReads->a+xLeftBeg, xLeftLen,
    target_beg, target_end, targetUid, x_max_countLeft, x_min_countLeft, 1, Hap_rate, is_local,
    position_index, reverse_sources, read_g, ruIndex, &n_matchLenLeft, &x_max_countLeft, 
    &x_min_countLeft);

    n_matchLenRight = x_max_countRight = x_min_countRight = 0;
    determin_hap_alignment_boundary_single_side(xReads->a+xRightBeg, xRightLen,
    target_beg, target_end, targetUid, x_max_countRight, x_min_countRight, 0, Hap_rate, is_local,
    position_index, reverse_sources, read_g, ruIndex, &n_matchLenRight, &x_max_countRight,
    &x_min_countRight);

    if(x_max_countLeft >= x_max_countRight)
    {
        determin_hap_alignment_boundary_single_side(xReads->a+xRightBeg, xRightLen,
        target_beg, target_end, targetUid, x_max_countLeft, x_min_countLeft, 0, Hap_rate, is_local,
        position_index, reverse_sources, read_g, ruIndex, &n_matchLenRight, &x_max_countRight,
        &x_min_countRight);
    }
    else
    {
        determin_hap_alignment_boundary_single_side(xReads->a+xLeftBeg, xLeftLen,
        target_beg, target_end, targetUid, x_max_countRight, x_min_countRight, 1, Hap_rate, is_local,
        position_index, reverse_sources, read_g, ruIndex, &n_matchLenLeft, &x_max_countLeft, 
        &x_min_countLeft);
    }

    (*x_interval_beg) = xLeftBeg + xLeftLen; (*x_interval_beg) -= n_matchLenLeft;
    (*x_interval_end) = xRightBeg + n_matchLenRight; (*x_interval_end) -= 1;
}

void get_hap_alignment_boundary(ma_utg_t* xReads, ma_utg_t* yReads, uint32_t type, 
uint32_t xLeftMatch, uint32_t xLeftTotal, uint32_t yLeftMatch, uint32_t yLeftTotal, 
uint32_t xRightMatch, uint32_t xRightTotal, uint32_t yRightMatch, uint32_t yRightTotal,
uint32_t xLeftBeg, uint32_t xLeftLen, uint32_t yLeftBeg, uint32_t yLeftLen,
uint32_t xRightBeg, uint32_t xRightLen, uint32_t yRightBeg, uint32_t yRightLen,
uint32_t xUid, uint32_t yUid, float Hap_rate, int is_local, uint64_t* position_index,
ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, uint32_t rev,
long long* r_x_interval_beg, long long* r_x_interval_end, 
long long* r_y_interval_beg, long long* r_y_interval_end)
{
    
    uint32_t x_max_count, x_min_count, y_max_count, y_min_count, n_matchLen;
    long long x_interval_beg, x_interval_end, y_interval_beg, y_interval_end;
    long long target_beg, target_end;
    x_max_count = x_min_count = y_max_count = y_min_count = 0;
    if(type == X2Y)
    {
        /********************x*********************/
        x_max_count = xRightMatch;
        x_min_count = xRightTotal;

        modify_target_interval(yLeftBeg, yLeftBeg+yLeftLen-1, yReads->n, &target_beg, &target_end);
        determin_hap_alignment_boundary_single_side(xReads->a+xLeftBeg, xLeftLen, 
        /**yLeftBeg, yLeftBeg+yLeftLen-1,**/ target_beg, target_end, yUid, 
        x_max_count, x_min_count, 1, Hap_rate, is_local, position_index, reverse_sources, 
        read_g, ruIndex, &n_matchLen, &x_max_count, &x_min_count);

        x_interval_beg = xLeftBeg + xLeftLen; x_interval_beg -= n_matchLen;
        x_interval_end = xRightBeg + xRightLen; x_interval_end -= 1;
        /********************x*********************/

        /********************y*********************/
        y_max_count = yLeftMatch;
        y_min_count = yLeftTotal;
        modify_target_interval(xRightBeg, xRightBeg+xRightLen-1, xReads->n, &target_beg, &target_end);
        determin_hap_alignment_boundary_single_side(yReads->a+yRightBeg, yRightLen,
        /**xRightBeg, xRightBeg+xRightLen-1,**/ target_beg, target_end, xUid, 
        y_max_count, y_min_count, rev, Hap_rate, is_local, position_index, reverse_sources, 
        read_g, ruIndex, &n_matchLen, &y_max_count, &y_min_count);
        if(rev == 0)
        {
            y_interval_beg = yLeftBeg;
            y_interval_end = yRightBeg + n_matchLen; y_interval_end -= 1;
        }
        else
        {
            y_interval_beg = yRightBeg + yRightLen; y_interval_beg -= n_matchLen;
            y_interval_end = yLeftBeg + yLeftLen; y_interval_end -= 1;
        }
        
        /********************y*********************/
    }
    else if(type == Y2X)
    {
        /********************x*********************/
        x_max_count = xLeftMatch;
        x_min_count = xLeftTotal;
        modify_target_interval(yRightBeg, yRightBeg+yRightLen-1, yReads->n, &target_beg, &target_end);
        determin_hap_alignment_boundary_single_side(xReads->a+xRightBeg, xRightLen,
        /**yRightBeg, yRightBeg+yRightLen-1,**/ target_beg, target_end, yUid, 
        x_max_count, x_min_count, 0, Hap_rate, is_local, position_index, reverse_sources, 
        read_g, ruIndex, &n_matchLen, &x_max_count, &x_min_count);

        x_interval_beg = xLeftBeg;
        x_interval_end = xRightBeg + n_matchLen; x_interval_end -= 1;
        /********************x*********************/

        /********************y*********************/
        y_max_count = yRightMatch;
        y_min_count = yRightTotal;
        modify_target_interval(xLeftBeg, xLeftBeg+xLeftLen-1, xReads->n, &target_beg, &target_end);
        determin_hap_alignment_boundary_single_side(yReads->a+yLeftBeg, yLeftLen,
        /**xLeftBeg, xLeftBeg+xLeftLen-1,**/ target_beg, target_end, xUid, 
        y_max_count, y_min_count, 1-rev, Hap_rate, is_local, position_index, reverse_sources, 
        read_g, ruIndex, &n_matchLen, &y_max_count, &y_min_count);
        if(rev == 0)
        {
            y_interval_beg = yLeftBeg + yLeftLen; y_interval_beg -= n_matchLen;
            y_interval_end = yRightBeg + yRightLen; y_interval_end -= 1;
        }
        else
        {
            y_interval_beg = yRightBeg;
            y_interval_end = yLeftBeg + n_matchLen; y_interval_end -= 1;
        }
        /********************y*********************/
    }
    else if(type == XCY)
    {
        /********************x*********************/   
        bi_direction_hap_alignment_extention(xReads, xLeftBeg, xLeftLen, xRightBeg, xRightLen, 
        yUid, 0, yReads->n - 1, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, 0,
        &x_interval_beg, &x_interval_end);
        /********************x*********************/

        /********************y*********************/
        y_interval_beg = 0;
        y_interval_end = yReads->n; y_interval_end -= 1;
        /********************y*********************/
    }
    else if(type == YCX)
    {
        /********************x*********************/  
        x_interval_beg = 0;
        x_interval_end = xReads->n; x_interval_end -= 1;
        /********************x*********************/  

        /********************y*********************/
        bi_direction_hap_alignment_extention(yReads, yLeftBeg, yLeftLen, yRightBeg, yRightLen,
        xUid, 0, xReads->n - 1, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, rev,
        &y_interval_beg, &y_interval_end);
        /********************y*********************/
    } else abort();

    (*r_x_interval_beg) = x_interval_beg;
    (*r_x_interval_end) = x_interval_end;
    (*r_y_interval_beg) = y_interval_beg;
    (*r_y_interval_end) = y_interval_end;    
}

uint32_t vote_overlap_type(kvec_asg_arc_t_offset* u_buffer, hap_candidates* hap_can, 
uint64_t* position_index, ma_utg_t* xReads, ma_utg_t* yReads)
{
    uint32_t i, xBasePos, yBasePos;
    asg_arc_t_offset* arch = NULL;
    uint32_t flag[4];
    flag[X2Y] = flag[Y2X] = flag[XCY] = flag[YCX] = 0;


    for (i = hap_can->index_beg; i <= hap_can->index_end; i++)
    {
        arch = &(u_buffer->a.a[i]);
        xBasePos = (uint32_t)(arch->Off>>32);
        yBasePos = (uint32_t)(arch->Off);
        flag[classify_hap_overlap(xBasePos, xBasePos, xReads->len, yBasePos, yBasePos, yReads->len,
        NULL, NULL, NULL, NULL)]++;
    }
    
    uint32_t max_flag_i = 0;
    for (i = 0; i < 4; i++)
    {
        if(i == max_flag_i) continue;
        if(flag[i] > flag[max_flag_i])
        {
            max_flag_i = i;
        }
    }
    
    return max_flag_i;
}


void get_base_boundary(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, 
asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads,
uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, long long yEndIndex, 
uint32_t dir, uint32_t rev, uint32_t* x_off, uint32_t* y_off)
{
    long long k, j, offset;
    ma_hit_t_alloc *xR = NULL;
    ma_hit_t *h = NULL;
    ma_sub_t *sq = NULL, *st = NULL;
    int32_t r;
    asg_arc_t t;
    uint32_t rId, Hap_uId, is_Unitig, v, w, v_dir, w_dir, is_found = 0, oLen = 0;
    uint64_t tmp;
    (*x_off) = (*y_off) = (uint32_t)-1;
    if(dir == 1)
    {
        for (k = xEndIndex; k >= xBegIndex; k--)
        {
            xR = &(reverse_sources[xReads->a[k]>>33]);
            is_found = 0; oLen = 0;
            for (j = 0; j < xR->length; j++)
            {
                h = &(xR->buffer[j]);
                sq = &(coverage_cut[Get_qn(*h)]);
                st = &(coverage_cut[Get_tn(*h)]);
                if(st->del || read_g->seq[Get_tn(*h)].del) continue;

                r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, 
                                asm_opt.max_hang_rate, min_ovlp, &t);
                ///if it is a contained overlap, skip
                if(r < 0) continue;

                rId = t.v>>1;
                if(read_g->seq[rId].del == 1) continue;
                ///there are two cases: 
                ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id  
                ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1
                get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig);
                if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
                if(Hap_uId != yUid) continue;

                v = xReads->a[k]>>32;
                get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig);
                if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
                if(Hap_uId != xUid) continue;
                if((uint32_t)(position_index[v>>1]) != k) continue;

                w = (yReads->a[(uint32_t)(position_index[rId])])>>32;

                v_dir = ((t.ul>>32)==v)?1:0;
                w_dir = (t.v == w)?1:0;
                if(rev == 0 && v_dir != w_dir) continue;
                if(rev == 1 && v_dir == w_dir) continue;

                /****************************may have bugs********************************/
                offset = (uint32_t)(position_index[rId]);
                if(offset < yBegIndex || offset > yEndIndex) continue;
                /****************************may have bugs********************************/

                tmp = get_xy_pos(read_g, &t, v, w, xReads->len, yReads->len, position_index, &(t.el));
                if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue;

                ///if(is_found == 0 || ((uint32_t)(tmp>>32) > (*x_off) && ((uint32_t)tmp) > (*y_off)))
                if(is_found == 0 || t.ol > oLen)
                {
                    (*x_off) = tmp>>32;
                    (*y_off) = (uint32_t)tmp;
                    oLen = t.ol;
                }

                is_found = 1; 
            }
            if(is_found) return;
        }   
    }
    else
    {
        for (k = xBegIndex; k <= xEndIndex; k++)
        {
            xR = &(reverse_sources[xReads->a[k]>>33]);
            is_found = 0; oLen = 0;
            for (j = 0; j < xR->length; j++)
            {
                h = &(xR->buffer[j]);
                sq = &(coverage_cut[Get_qn(*h)]);
                st = &(coverage_cut[Get_tn(*h)]);
                if(st->del || read_g->seq[Get_tn(*h)].del) continue;

                r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, 
                                asm_opt.max_hang_rate, min_ovlp, &t);
                ///if it is a contained overlap, skip
                if(r < 0) continue;

                rId = t.v>>1;
                if(read_g->seq[rId].del == 1) continue;
                ///there are two cases: 
                ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id  
                ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1
                get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig);
                if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
                if(Hap_uId != yUid) continue;

                v = xReads->a[k]>>32;
                get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig);
                if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
                if(Hap_uId != xUid) continue;
                if((uint32_t)(position_index[v>>1]) != k) continue;

                w = (yReads->a[(uint32_t)(position_index[rId])])>>32;

                v_dir = ((t.ul>>32)==v)?1:0;
                w_dir = (t.v == w)?1:0;
                if(rev == 0 && v_dir != w_dir) continue;
                if(rev == 1 && v_dir == w_dir) continue;
            
                /****************************may have bugs********************************/
                offset = (uint32_t)(position_index[rId]);
                if(offset < yBegIndex || offset > yEndIndex) continue;
                /****************************may have bugs********************************/

                tmp = get_xy_pos(read_g, &t, v, w, xReads->len, yReads->len, position_index, &(t.el));
                if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue;

                ///if(is_found == 0 || ((uint32_t)(tmp>>32) < (*x_off) && ((uint32_t)tmp) < (*y_off)))
                if(is_found == 0 || t.ol > oLen)
                {
                    (*x_off) = tmp>>32;
                    (*y_off) = (uint32_t)tmp;
                    oLen = t.ol;
                }  

                is_found = 1;         
            }
            if(is_found) return;
        }
    }


    (*x_off) = (*y_off) = (uint32_t)-1;
    
}

void print_asg_arc_t_offset(asg_arc_t_offset* x, long long n, const char* info)
{
    fprintf(stderr,"\n\n(%s)n: %lld\n", info, n);
    long long i, x_off, y_off;
    for (i = 0; i < n; i++)
    {
        x_off = (long long)(x[i].Off>>32);
        y_off = (long long)((uint32_t)x[i].Off);
        fprintf(stderr, "i: %lld, x_off: %lld, y_off: %lld, weight: %lu, rev: %u, ol: %u\n", 
        i, x_off, y_off, (unsigned long)x[i].weight, x[i].x.el, x[i].x.ol);
    }
}


// Binary search 
inline int GetCeilIndex(asg_arc_t_offset* arr, kvec_t_i32_warp* T, int l, int r, uint32_t key) 
{ 
    while (r - l > 1) { 
        int m = l + (r - l) / 2; 
        if (Get_yOff(arr[T->a.a[m]].Off) >= key) 
            r = m; 
        else
            l = m; 
    } 
  
    return r; 
} 


void quick_LIS(asg_arc_t_offset* x, uint32_t n, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex)
{
    tailIndex->a.n = prevIndex->a.n = 0;
    if(n == 0) return;

    kv_resize(int32_t, tailIndex->a, n);
    kv_resize(int32_t, prevIndex->a, n);

    long long len = 1, i, pos, m; ///the length of chain must be >=1
    tailIndex->a.a[0] = 0;
    prevIndex->a.a[0] = -1;

    ///x has already sorted by x_pos
    for(i = 1; i < (long long)n; i++) 
    { 
        if(Get_yOff(x[i].Off) < Get_yOff(x[tailIndex->a.a[0]].Off))  
        { 
            // new smallest value 
            tailIndex->a.a[0] = i; ///doesn't matter too much
        } 
        else if(Get_yOff(x[i].Off) > Get_yOff(x[tailIndex->a.a[len - 1]].Off)) 
        { 
            // arr[i] wants to extend largest subsequence 
            prevIndex->a.a[i] = tailIndex->a.a[len - 1]; 
            tailIndex->a.a[len++] = i; 
        } 
        else 
        { 
            // arr[i] wants to be a potential condidate of 
            // future subsequence 
            // It will replace ceil value in tailIndices 
            pos = GetCeilIndex(x, tailIndex, -1, len - 1, Get_yOff(x[i].Off)); 
            prevIndex->a.a[i] = pos > 0? tailIndex->a.a[pos - 1] : -1; 
            tailIndex->a.a[pos] = i; 
        } 
    }


    for (m = 0, i = tailIndex->a.a[len - 1]; m < len; i = prevIndex->a.a[i], m++)
    {
        tailIndex->a.a[len-m-1] = i;
    }

    tailIndex->a.n = len;
}

uint64_t get_xy_pos_by_pos(asg_t *read_g, asg_arc_t* t, uint32_t v_in_unitig, uint32_t w_in_unitig, 
uint32_t v_in_pos, uint32_t w_in_pos, uint32_t xUnitigLen, uint32_t yUnitigLen, uint8_t* rev)
{
    uint32_t x_pos, y_pos, x_dir = 0, y_dir = 0;
    uint64_t tmp;
    x_pos = y_pos = (uint32_t)-1;
    if((t->ul>>32)==v_in_unitig)///end pos
    {
        x_pos = v_in_pos + read_g->seq[v_in_unitig>>1].len - 1;
        x_dir = 0;
    }
    else if((t->ul>>32)==(v_in_unitig^1))///start pos
    {
        x_pos = v_in_pos;
        x_dir = 1;
    }
    else
    {
        fprintf(stderr, "ERROR\n");
    }
    
    if(t->v == w_in_unitig)
    {
        y_pos = w_in_pos + t->ol - 1;
        y_dir = 0;
    }
    else if(t->v == (w_in_unitig^1))
    {
        y_pos = w_in_pos + read_g->seq[w_in_unitig>>1].len - t->ol;
        y_dir = 1;
    }
    else
    {
        fprintf(stderr, "ERROR\n");
    }

    (*rev) = x_dir^y_dir;
    if((*rev))
    {
        if(yUnitigLen <= y_pos)
        {
            y_pos = (uint32_t)-1;
        }
        else
        {
            y_pos = yUnitigLen - y_pos - 1;
        }
    } 

    if(x_pos>=xUnitigLen) x_pos = (uint32_t)-1;
    if(y_pos>=yUnitigLen) y_pos = (uint32_t)-1;
    
    tmp = x_pos; tmp = tmp << 32; tmp = tmp | y_pos;
    return tmp;
}

void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd)
{
    ma_hit_t_alloc* reverse_sources = cov->reverse_sources;
    ma_sub_t *coverage_cut = cov->coverage_cut;
    int max_hang = cov->max_hang;
    int min_ovlp = cov->min_ovlp;
    kvec_asg_arc_t_offset* u_buffer = &(cov->u_buffer);
    kvec_t_i32_warp* tailIndex = &(cov->tailIndex);
    kvec_t_i32_warp* prevIndex = &(cov->prevIndex);
    ma_hit_t_alloc *xR = NULL;
    ma_hit_t *h = NULL;
    ma_sub_t *sq = NULL, *st = NULL;
    int32_t r;
    asg_arc_t t;
    uint32_t rId, v, w;
    uint64_t tmp;
    asg_arc_t_offset t_offset;
    u_buffer->a.n = 0;
    ///(*xEnd) = (uint32_t)-1;
    (*xEnd) = 0;
    uint32_t u_i, r_i, k, j, len, p_v, *a = xReads->b.a, uid, ori, l, m, aOcc, nv, xOcc = (uint32_t)-1;
    ma_utg_t* u = NULL;
    asg_arc_t *av  = NULL;


    for (u_i = r_i = len = aOcc = 0, xOcc = (uint32_t)-1, p_v = (uint32_t)-1; u_i < xReads->b.n; u_i++)
    {
        uid = a[u_i] >> 1;
        ori = a[u_i] & 1;
        u = &(ug->u.a[uid]);
        if(u->n == 0) continue;

        for (r_i = 0; r_i < u->n; r_i++, aOcc++)
        {
            l = 0;
            v = (ori == 1?((uint64_t)((u->a[u->n - r_i - 1])^(uint64_t)(0x100000000)))>>32:((uint64_t)(u->a[r_i]))>>32);

            if(p_v != (uint32_t)-1)
            {
                av = asg_arc_a(read_sg, p_v);
                nv = asg_arc_n(read_sg, p_v);
                
                for (k = 0; k < nv; k++)
                {
                    if(av[k].del) continue;
                    if(av[k].v == v) 
                    {
                        l = asg_arc_len(av[k]);
                        break;
                    }
                }
                if(k == nv) fprintf(stderr, "ERROR\n");
            }

            p_v = v; len += l;
            if(len >= targetBaseLen)
            {
                xOcc = aOcc;
                break;
            }
        }

        if(xOcc != (uint32_t)-1) break;
    }
    if(xOcc == (uint32_t)-1) xOcc = aOcc;
    if(xOcc == 0) xOcc = 1;



    for (u_i = r_i = len = aOcc = 0, p_v = (uint32_t)-1; u_i < xReads->b.n; u_i++)
    {
        uid = a[u_i] >> 1;
        ori = a[u_i] & 1;
        u = &(ug->u.a[uid]);
        if(u->n == 0) continue;

        for (r_i = 0; r_i < u->n; r_i++, aOcc++)
        {
            if(aOcc >= xOcc) break;
            l = 0;
            v = (ori == 1?((uint64_t)((u->a[u->n - r_i - 1])^(uint64_t)(0x100000000)))>>32:((uint64_t)(u->a[r_i]))>>32);

            if(p_v != (uint32_t)-1)
            {
                av = asg_arc_a(read_sg, p_v);
                nv = asg_arc_n(read_sg, p_v);
                
                for (k = 0; k < nv; k++)
                {
                    if(av[k].del) continue;
                    if(av[k].v == v) 
                    {
                        l = asg_arc_len(av[k]);
                        break;
                    }
                }
                if(k == nv) fprintf(stderr, "ERROR\n");
            }

            p_v = v; len += l;

            xR = &(reverse_sources[v>>1]);
            for (j = 0; j < xR->length; j++)
            {
                h = &(xR->buffer[j]);
                sq = &(coverage_cut[Get_qn(*h)]);
                st = &(coverage_cut[Get_tn(*h)]);
                if(st->del || read_sg->seq[Get_tn(*h)].del) continue;
                r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, 
                            asm_opt.max_hang_rate, min_ovlp, &t);
                ///if it is a contained overlap, skip
                if(r < 0) continue;

                rId = t.v>>1;
                if(read_sg->seq[rId].del == 1) continue;
                if(cov->pos_idx[rId] == (uint64_t)-1) continue;
                w = (uint32_t)(cov->pos_idx[rId]);
                if(rId != (w>>1)) continue;

                tmp = get_xy_pos_by_pos(read_sg, &t, v, w, len, cov->pos_idx[w>>1]>>32,
                                                    (uint32_t)-1, targetBaseLen, &(t.el));
                if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue;
                if(t.el) continue; ///must

                t_offset.Off = tmp;
                t_offset.x = t;
                t_offset.weight = 1;
                kv_push(asg_arc_t_offset, u_buffer->a, t_offset); 
            }
        }

        if(aOcc >= xOcc) break;
    }

    if(u_buffer->a.n == 0) return;

    qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment_chaining);

    ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "before");


    for (k = 1, l = 0, m = 0; k <= u_buffer->a.n; ++k)
    {
        if (k == u_buffer->a.n || u_buffer->a.a[k].x.el != u_buffer->a.a[l].x.el || 
                                            u_buffer->a.a[k].Off != u_buffer->a.a[l].Off) 
        {
            u_buffer->a.a[m] = u_buffer->a.a[l];
            for (l += 1; l < k; l++)
            {
                u_buffer->a.a[m].weight += u_buffer->a.a[l].weight;
                if(u_buffer->a.a[l].x.ol > u_buffer->a.a[m].x.ol)
                {
                    u_buffer->a.a[m].x = u_buffer->a.a[l].x;
                }
            }
            l = k;
            m++;
        }
    } 
    u_buffer->a.n = m;

    ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "after");
    quick_LIS(u_buffer->a.a, u_buffer->a.n, tailIndex, prevIndex);
    if(tailIndex->a.n == 0) return;


    uint32_t xLen_thres = (uint32_t)-1;
    asg_arc_t_offset* best = &(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]]);
    for (u_i = r_i = len = aOcc = 0, p_v = (uint32_t)-1; u_i < xReads->b.n; u_i++)
    {
        uid = a[u_i] >> 1;
        ori = a[u_i] & 1;
        u = &(ug->u.a[uid]);
        if(u->n == 0) continue;

        for (r_i = 0; r_i < u->n; r_i++, aOcc++)
        {
            l = 0;
            v = (ori == 1?((uint64_t)((u->a[u->n - r_i - 1])^(uint64_t)(0x100000000)))>>32:((uint64_t)(u->a[r_i]))>>32);

            if(p_v != (uint32_t)-1)
            {
                av = asg_arc_a(read_sg, p_v);
                nv = asg_arc_n(read_sg, p_v);
                
                for (k = 0; k < nv; k++)
                {
                    if(av[k].del) continue;
                    if(av[k].v == v) 
                    {
                        l = asg_arc_len(av[k]);
                        break;
                    }
                }
                if(k == nv) fprintf(stderr, "ERROR\n");
            }

            p_v = v; len += l;

            if((v>>1) == (best->x.ul>>33) && xLen_thres == (uint32_t)-1)
            {
                ///cov->pos_idx[v>>1] = len; 
                xR = &(reverse_sources[v>>1]);
                for (j = 0; j < xR->length; j++)
                {
                    h = &(xR->buffer[j]);
                    sq = &(coverage_cut[Get_qn(*h)]);
                    st = &(coverage_cut[Get_tn(*h)]);
                    if(st->del || read_sg->seq[Get_tn(*h)].del) continue;
                    r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, 
                                asm_opt.max_hang_rate, min_ovlp, &t);
                    ///if it is a contained overlap, skip
                    if(r < 0) continue;

                    rId = t.v>>1;
                    if(read_sg->seq[rId].del == 1) continue;
                    if(cov->pos_idx[rId] == (uint64_t)-1) continue;
                    w = (uint32_t)(cov->pos_idx[rId]);
                    if(rId != (w>>1)) continue;

                    tmp = get_xy_pos_by_pos(read_sg, &t, v, w, len, cov->pos_idx[w>>1]>>32,
                                                        (uint32_t)-1, targetBaseLen, &(t.el));
                    if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue;
                    if(t.el) continue; ///must

                    t_offset.Off = tmp;
                    t_offset.x = t;
                    t_offset.weight = 1;
                    if(t_offset.Off == best->Off && t_offset.x.v == best->x.v && t_offset.x.ul == best->x.ul)
                    {
                        xLen_thres = Get_xOff(best->Off) + targetBaseLen - Get_yOff(best->Off);
                    }
                }
            }

            if(len >= xLen_thres)
            {
                (*xEnd) = aOcc;
                return;
            }
        }
    }

    (*xEnd) = aOcc;
}


void get_base_boundary_advance_back(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, 
asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads,
uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, long long yEndIndex, 
uint32_t rev, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex,
uint32_t* xBeg, uint32_t* xEnd, uint32_t* yBeg, uint32_t* yEnd)
{
    long long k, j, offset, m;
    ma_hit_t_alloc *xR = NULL;
    ma_hit_t *h = NULL;
    ma_sub_t *sq = NULL, *st = NULL;
    int32_t r;
    asg_arc_t t;
    uint32_t rId, Hap_uId, is_Unitig, v, w, v_dir, w_dir;
    uint64_t tmp;
    asg_arc_t_offset t_offset;
    u_buffer->a.n = 0;
    (*xBeg) = (*xEnd) = (*yBeg) = (*yEnd) = (uint32_t)-1;
    for (k = xBegIndex; k <= xEndIndex; k++)
    {
        xR = &(reverse_sources[xReads->a[k]>>33]);
        for (j = 0; j < xR->length; j++)
        {
            h = &(xR->buffer[j]);
            sq = &(coverage_cut[Get_qn(*h)]);
            st = &(coverage_cut[Get_tn(*h)]);
            if(st->del || read_g->seq[Get_tn(*h)].del) continue;

            r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, 
                            asm_opt.max_hang_rate, min_ovlp, &t);
            ///if it is a contained overlap, skip
            if(r < 0) continue;

            rId = t.v>>1;
            if(read_g->seq[rId].del == 1) continue;
            ///there are two cases: 
            ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id  
            ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1
            get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig);
            if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
            if(Hap_uId != yUid) continue;

            v = xReads->a[k]>>32;
            get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig);
            if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
            if(Hap_uId != xUid) continue;
            if((uint32_t)(position_index[v>>1]) != k) continue;

            w = (yReads->a[(uint32_t)(position_index[rId])])>>32;

            v_dir = ((t.ul>>32)==v)?1:0;
            w_dir = (t.v == w)?1:0;
            if(rev == 0 && v_dir != w_dir) continue;
            if(rev == 1 && v_dir == w_dir) continue;
        
            /****************************may have bugs********************************/
            offset = (uint32_t)(position_index[rId]);
            if(offset < yBegIndex || offset > yEndIndex) continue;
            /****************************may have bugs********************************/

            tmp = get_xy_pos(read_g, &t, v, w, xReads->len, yReads->len, position_index, &(t.el));
            if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue;

            t_offset.Off = tmp;
            t_offset.x = t;
            t_offset.weight = 1;
            kv_push(asg_arc_t_offset, u_buffer->a, t_offset);      
        }
    }
    if(u_buffer->a.n == 0) return;

    qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment_chaining);

    ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "before");

    for (k = 1, m = 1; k < (long long)u_buffer->a.n; k++)
    {
        if(u_buffer->a.a[m-1].Off == u_buffer->a.a[k].Off)
        {
            u_buffer->a.a[m-1].weight += u_buffer->a.a[k].weight;
            if(u_buffer->a.a[k].x.ol > u_buffer->a.a[m-1].x.ol)
            {
                u_buffer->a.a[m-1].x = u_buffer->a.a[k].x;
            }
            continue;
        }
        u_buffer->a.a[m] = u_buffer->a.a[k];
        m++;
    }
    u_buffer->a.n = m;

    ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "after");
    quick_LIS(u_buffer->a.a, u_buffer->a.n, tailIndex, prevIndex);

    if(tailIndex->a.n == 0) return;

    (*xBeg) = Get_xOff(u_buffer->a.a[tailIndex->a.a[0]].Off);
    (*yBeg) = Get_yOff(u_buffer->a.a[tailIndex->a.a[0]].Off);
    (*xEnd) = Get_xOff(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].Off);
    (*yEnd) = Get_yOff(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].Off);
}

uint32_t determine_hap_overlap_type_advance_back(hap_candidates* hap_can, ma_utg_t *xReads, ma_utg_t *yReads,
R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, uint64_t* position_index, 
int max_hang, int min_ovlp, uint32_t xUid, uint32_t yUid, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, 
kvec_t_i32_warp* prevIndex, long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end)
{
    uint32_t x_pos_beg, y_pos_beg, x_pos_end, y_pos_end;
    /*************************x***************************/
    get_base_boundary_advance_back(ruIndex, reverse_sources, coverage_cut, read_g, position_index, 
    max_hang, min_ovlp, xReads, yReads, xUid, yUid, Get_x_beg(*hap_can), Get_x_end(*hap_can), 
    Get_y_beg(*hap_can), Get_y_end(*hap_can), Get_rev(*hap_can), u_buffer, tailIndex, prevIndex, 
    &x_pos_beg, &x_pos_end, &y_pos_beg, &y_pos_end);
    /*************************x***************************/

    if(x_pos_beg == (uint32_t)-1 || y_pos_beg == (uint32_t)-1 
                    || x_pos_end == (uint32_t)-1 || y_pos_end == (uint32_t)-1)
    {
        return (uint32_t)-1;
    }
    if(x_pos_beg > x_pos_end || y_pos_beg > y_pos_end) return (uint32_t)-1;
    
    /**
    #define X2Y 0
    #define Y2X 1
    #define XCY 2
    #define YCX 3
    **/
    return classify_hap_overlap(x_pos_beg, x_pos_end, xReads->len, y_pos_beg, y_pos_end, yReads->len, 
    r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end);
}

void get_idx_by_base(ma_utg_t *x, asg_t *read_g, long long beg_base, long long end_base, 
                                            long long* beg_idx, long long* end_idx)
{
    long long offset, r_beg, r_end;
    uint64_t i, rId;
    (*beg_idx) = (*end_idx) = -1;
    for (i = 0, offset = 0; i < x->n; i++)
    {
        rId = x->a[i]>>33;
        r_beg = offset; r_end = offset + (long long)(read_g->seq[rId].len) - 1;
        offset += (uint32_t)x->a[i];
        if(beg_base > r_end || r_beg > end_base)
        {
            if((*beg_idx) != -1 && (*end_idx) != -1) break;
            continue;
        } 
        if((*beg_idx) == -1) (*beg_idx) = i;
        (*end_idx) = i;
    }
}

int get_base_boundary_chain(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, 
asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads,
uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, long long yEndIndex, 
uint32_t rev, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex)
{
    long long k, j, l, offset, m;
    ma_hit_t_alloc *xR = NULL;
    ma_hit_t *h = NULL;
    ma_sub_t *sq = NULL, *st = NULL;
    int32_t r;
    asg_arc_t t;
    uint32_t rId, Hap_uId, is_Unitig, v, w, v_dir, w_dir;
    uint64_t tmp;
    asg_arc_t_offset t_offset;
    u_buffer->a.n = 0;
    for (k = xBegIndex; k <= xEndIndex; k++)
    {
        xR = &(reverse_sources[xReads->a[k]>>33]);
        for (j = 0; j < xR->length; j++)
        {
            h = &(xR->buffer[j]);
            sq = &(coverage_cut[Get_qn(*h)]);
            st = &(coverage_cut[Get_tn(*h)]);
            if(st->del || read_g->seq[Get_tn(*h)].del) continue;

            r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, 
                            asm_opt.max_hang_rate, min_ovlp, &t);
            ///if it is a contained overlap, skip
            if(r < 0) continue;

            rId = t.v>>1;
            if(read_g->seq[rId].del == 1) continue;
            ///there are two cases: 
            ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id  
            ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1
            get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig);
            if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
            if(Hap_uId != yUid) continue;

            v = xReads->a[k]>>32;
            get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig);
            if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
            if(Hap_uId != xUid) continue;
            if((uint32_t)(position_index[v>>1]) != k) continue;

            w = (yReads->a[(uint32_t)(position_index[rId])])>>32;

            v_dir = ((t.ul>>32)==v)?1:0;
            w_dir = (t.v == w)?1:0;
            if(rev == 0 && v_dir != w_dir) continue;
            if(rev == 1 && v_dir == w_dir) continue;
        
            /****************************may have bugs********************************/
            offset = (uint32_t)(position_index[rId]);
            if(offset < yBegIndex || offset > yEndIndex) continue;
            /****************************may have bugs********************************/

            tmp = get_xy_pos(read_g, &t, v, w, xReads->len, yReads->len, position_index, &(t.el));
            if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue;

            t_offset.Off = tmp;
            t_offset.x = t;
            t_offset.weight = 1;
            kv_push(asg_arc_t_offset, u_buffer->a, t_offset);      
        }
    }
    if(u_buffer->a.n == 0) return 0;

    qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment_chaining);

    ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "before");
    for (k = 1, l = 0, m = 0; k <= (long long)u_buffer->a.n; ++k)
    {
        if (k == (long long)u_buffer->a.n || u_buffer->a.a[k].Off != u_buffer->a.a[l].Off) 
        {
            u_buffer->a.a[m] = u_buffer->a.a[l];
            for (l += 1; l < k; l++)
            {
                u_buffer->a.a[m].weight += u_buffer->a.a[l].weight;
                if(u_buffer->a.a[l].x.ol > u_buffer->a.a[m].x.ol)
                {
                    u_buffer->a.a[m].x = u_buffer->a.a[l].x;
                }
            }
            l = k;
            m++;
        }
    } 
    u_buffer->a.n = m;

    ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "after");
    quick_LIS(u_buffer->a.a, u_buffer->a.n, tailIndex, prevIndex);

    if(tailIndex->a.n == 0) return 0;
    return 1;
}


void get_base_boundary_advance(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, 
asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads,
uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, long long yEndIndex, 
uint32_t rev, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex,
uint32_t* xBeg, uint32_t* xEnd, uint32_t* yBeg, uint32_t* yEnd)
{
    long long offset;
    long long new_xBeg, new_yBeg, new_xEnd, new_yEnd;
    long long new_xIdxBeg, new_yIdxBeg, new_xIdxEnd, new_yIdxEnd;

    (*xBeg) = (*xEnd) = (*yBeg) = (*yEnd) = (uint32_t)-1;
    if(!get_base_boundary_chain(ruIndex, reverse_sources, coverage_cut, read_g, position_index, 
    max_hang, min_ovlp, xReads, yReads, xUid, yUid, xBegIndex, xEndIndex, yBegIndex, yEndIndex, 
    rev, u_buffer, tailIndex, prevIndex))
    {
        return;
    }
    ///base
    new_xBeg = Get_xOff(u_buffer->a.a[tailIndex->a.a[0]].Off);
    new_yBeg = Get_yOff(u_buffer->a.a[tailIndex->a.a[0]].Off);
    new_xEnd = Get_xOff(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].Off);
    new_yEnd = Get_yOff(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].Off);

    if(new_xBeg > new_xEnd || new_yBeg > new_yEnd) return;
    
    classify_hap_overlap(new_xBeg, new_xEnd, xReads->len, new_yBeg, new_yEnd, yReads->len, 
                                                        &new_xBeg, &new_xEnd, &new_yBeg, &new_yEnd);

    if(rev) 
    {
        new_yBeg = yReads->len - new_yBeg - 1;
        new_yEnd = yReads->len - new_yEnd - 1;
        offset = new_yBeg; new_yBeg = new_yEnd; new_yEnd = offset;
    }
    ///idx
    get_idx_by_base(xReads, read_g, new_xBeg, new_xEnd, &new_xIdxBeg, &new_xIdxEnd);
    get_idx_by_base(yReads, read_g, new_yBeg, new_yEnd, &new_yIdxBeg, &new_yIdxEnd);
    if(new_xIdxBeg == -1 || new_xIdxEnd == -1 || new_yIdxBeg == -1 || new_yIdxEnd == -1) return;
    
    if(!get_base_boundary_chain(ruIndex, reverse_sources, coverage_cut, read_g, position_index, 
    max_hang, min_ovlp, xReads, yReads, xUid, yUid, new_xIdxBeg, new_xIdxEnd, new_yIdxBeg, 
    new_yIdxEnd, rev, u_buffer, tailIndex, prevIndex))
    {
        return;
    }

    (*xBeg) = Get_xOff(u_buffer->a.a[tailIndex->a.a[0]].Off);
    (*yBeg) = Get_yOff(u_buffer->a.a[tailIndex->a.a[0]].Off);
    (*xEnd) = Get_xOff(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].Off);
    (*yEnd) = Get_yOff(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].Off);
}

#define generic_key(x) (x)
KRADIX_SORT_INIT(i32, int32_t, generic_key, sizeof(int32_t))

long long get_chain_score(ma_utg_t *xReads, asg_t *read_g, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* idx, 
ma_hit_t_alloc* reverse_sources, long long xBegPos, long long xEndPos)
{
    long long offset, r_beg, r_end, inp_beg, inp_end, hap_beg, hap_end, inp_match, hap_match, ovlp;
    uint64_t i, k, rId;
    idx->a.n = 0;
    
    for (i = k = 0; i < tailIndex->a.n; i++)
    {
        rId = u_buffer->a.a[tailIndex->a.a[i]].x.ul>>33;
        
        
        for (; k < xReads->n; k++)
        {
            if(rId == (xReads->a[k]>>33)) break;
        }

        if(k >= xReads->n)
        {
            for (k = 0; k < xReads->n; k++)
            {
                if(rId == (xReads->a[k]>>33)) break;
            }
        }

        if(k < xReads->n) kv_push(int32_t, idx->a, k);
        else
        {
            fprintf(stderr, "\nERROR-get_chain_score: tailIndex->a.n: %lu, xReads->n: %lu\n", (uint64_t)tailIndex->a.n, (uint64_t)xReads->n);
        } 
    }
    

    radix_sort_i32(idx->a.a, idx->a.a + idx->a.n);
    inp_beg = -1; inp_end = -2;
    hap_beg = -1; hap_end = -2;
    for (i = k = 0, offset = 0, inp_match = hap_match = 0; i < xReads->n; i++)
    {
        rId = xReads->a[i]>>33;
        r_beg = offset; r_end = offset + (long long)(read_g->seq[rId].len) - 1;
        offset += (uint32_t)xReads->a[i];

        if(reverse_sources[rId].length > 0)
        {
            if(r_beg <= hap_end)
            {
                hap_end = MAX(hap_end, r_end);
            } 
            else
            {
                ///match += (hap_end - hap_beg + 1);
                ovlp = (long long)(MIN(hap_end, xEndPos)) - (long long)(MAX(hap_beg, xBegPos)) + 1;
                hap_match += (ovlp >= 0? ovlp : 0);
                hap_beg = r_beg; hap_end = r_end;
            }
        }

        for (; k < idx->a.n; k++)
        {
            if(i <= (uint64_t)idx->a.a[k]) break;
        }

        if(k >= idx->a.n) continue;
        
        if(i == (uint64_t)idx->a.a[k])
        {
            if(r_beg <= inp_end)
            {
                inp_end = MAX(inp_end, r_end);
            } 
            else
            {
                ovlp = (long long)(MIN(inp_end, xEndPos)) - (long long)(MAX(inp_beg, xBegPos)) + 1;
                inp_match += (ovlp >= 0? ovlp : 0);
                inp_beg = r_beg; inp_end = r_end;
            }
        }
    }
    
    ovlp = (long long)(MIN(inp_end, xEndPos)) - (long long)(MAX(inp_beg, xBegPos)) + 1;
    inp_match += (ovlp >= 0? ovlp : 0);
    
    ovlp = (long long)(MIN(hap_end, xEndPos)) - (long long)(MAX(hap_beg, xBegPos)) + 1;
    hap_match += (ovlp >= 0? ovlp : 0);

    // if(inp_match > (xEndPos - xBegPos + 1)) fprintf(stderr, "ERROR1\n");
    // if(hap_match > (xEndPos - xBegPos + 1)) fprintf(stderr, "ERRO2\n");
    // if(inp_match > hap_match) fprintf(stderr, "ERROR3\n");
    // fprintf(stderr, "tailIndex->a.n: %u, xReads->n: %u, total_match: %lld, hap_match: %lld, inp_match: %lld\n", 
    //                     tailIndex->a.n, xReads->n, (xEndPos - xBegPos + 1), hap_match, inp_match);
    
    return ((double)(inp_match)*CHAIN_MATCH) - ((double)(hap_match-inp_match)*CHAIN_UNMATCH);
}


uint32_t determine_hap_overlap_type_advance(hap_candidates* hap_can, ma_utg_t *xReads, ma_utg_t *yReads,
R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, uint64_t* position_index, 
int max_hang, int min_ovlp, uint32_t xUid, uint32_t yUid, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, 
kvec_t_i32_warp* prevIndex, long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end)
{
    uint32_t x_pos_beg, y_pos_beg, x_pos_end, y_pos_end;
    /*************************x***************************/
    get_base_boundary_advance(ruIndex, reverse_sources, coverage_cut, read_g, position_index, 
    max_hang, min_ovlp, xReads, yReads, xUid, yUid, Get_x_beg(*hap_can), Get_x_end(*hap_can), 
    Get_y_beg(*hap_can), Get_y_end(*hap_can), Get_rev(*hap_can), u_buffer, tailIndex, prevIndex, 
    &x_pos_beg, &x_pos_end, &y_pos_beg, &y_pos_end);
    /*************************x***************************/
    if(x_pos_beg == (uint32_t)-1 || y_pos_beg == (uint32_t)-1 
                    || x_pos_end == (uint32_t)-1 || y_pos_end == (uint32_t)-1)
    {
        return (uint32_t)-1;
    }
    if(x_pos_beg > x_pos_end || y_pos_beg > y_pos_end) return (uint32_t)-1;

    /**
    #define X2Y 0
    #define Y2X 1
    #define XCY 2
    #define YCX 3
    **/
    hap_can->index_end = classify_hap_overlap(x_pos_beg, x_pos_end, xReads->len, y_pos_beg, y_pos_end, yReads->len, 
    r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end);
    hap_can->x_beg_pos = MIN((uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[0]].x.ul>>33]), 
                                (uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].x.ul>>33]));
    hap_can->x_end_pos = MAX((uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[0]].x.ul>>33]), 
                                (uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].x.ul>>33]));
    hap_can->y_beg_pos = MIN((uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[0]].x.v>>1]), 
                                (uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].x.v>>1]));
    hap_can->y_end_pos = MAX((uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[0]].x.v>>1]), 
                                (uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].x.v>>1]));
    double xLeftMatch, xLeftTotal;
    get_pair_hap_similarity(xReads->a + hap_can->x_beg_pos, hap_can->x_end_pos + 1 - hap_can->x_beg_pos, 
        yUid, reverse_sources, read_g, ruIndex, &xLeftMatch, &xLeftTotal);
    if(xLeftMatch == 0 || xLeftTotal == 0) return (uint32_t)-1;
    hap_can->weight = xLeftMatch;
    hap_can->index_beg = xLeftTotal;
    hap_can->score = get_chain_score(xReads, read_g, u_buffer, tailIndex, prevIndex, reverse_sources, 
                                                                        (*r_x_pos_beg), (*r_x_pos_end));
    if(hap_can->score <= 0) return (uint32_t)-1;
    return hap_can->index_end;
}



uint32_t determine_hap_overlap_type(hap_candidates* hap_can, ma_utg_t *xReads, ma_utg_t *yReads,
R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, 
uint64_t* position_index, int max_hang, int min_ovlp, uint32_t xUid, uint32_t yUid,
long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end)
{
    uint32_t x_pos_beg, y_pos_beg, x_pos_end, y_pos_end;
    /*************************x***************************/
    get_base_boundary(ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, 
    min_ovlp, xReads, yReads, xUid, yUid, Get_x_beg(*hap_can), Get_x_end(*hap_can), 
    Get_y_beg(*hap_can), Get_y_end(*hap_can), 0, Get_rev(*hap_can), &x_pos_beg, &y_pos_beg);

    get_base_boundary(ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, 
    min_ovlp, xReads, yReads, xUid, yUid, Get_x_beg(*hap_can), Get_x_end(*hap_can), 
    Get_y_beg(*hap_can), Get_y_end(*hap_can), 1, Get_rev(*hap_can), &x_pos_end, &y_pos_end);
    /*************************x***************************/

    if(x_pos_beg == (uint32_t)-1 || y_pos_beg == (uint32_t)-1 
                    || x_pos_end == (uint32_t)-1 || y_pos_end == (uint32_t)-1)
    {
        return (uint32_t)-1;
    }
    if(x_pos_beg > x_pos_end || y_pos_beg > y_pos_end) return (uint32_t)-1;
    
    /**
    #define X2Y 0
    #define Y2X 1
    #define XCY 2
    #define YCX 3
    **/
    return classify_hap_overlap(x_pos_beg, x_pos_end, xReads->len, y_pos_beg, y_pos_end, yReads->len, 
    r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end);
}


uint32_t calculate_pair_hap_similarity_advance(hap_candidates* hap_can, 
uint64_t* position_index, uint32_t xUid, uint32_t yUid, ma_utg_t* xReads, ma_utg_t* yReads,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, ma_sub_t *coverage_cut, 
float Hap_rate, int is_local, int max_hang, int min_ovlp, uint64_t cov_threshold, kvec_asg_arc_t_offset* u_buffer, 
kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, hap_cov_t *cov, long long* r_x_pos_beg, long long* r_x_pos_end, 
long long* r_y_pos_beg, long long* r_y_pos_end)
{
    uint32_t max_count = 0, min_count = 0, flag;
    uint32_t xLen = xReads->n, xIndex;
    uint32_t yLen = yReads->n, yIndex;
    uint32_t xLeftBeg, xLeftLen, yLeftBeg, yLeftLen;
    uint32_t xRightBeg, xRightLen, yRightBeg, yRightLen;
    double xLeftMatch = 0, xLeftTotal = 0, yLeftMatch = 0, yLeftTotal = 0;
    double xRightMatch = 0, xRightTotal = 0, yRightMatch = 0, yRightTotal = 0;
    asg_arc_t* arch = NULL;
    

    arch = &(hap_can->t);
    xIndex = (uint32_t)(position_index[arch->ul>>33]);
    yIndex = (uint32_t)(position_index[arch->v>>1]);
    
    if(hap_can->rev == 0)
    {
        xLeftBeg = 0; xLeftLen = xIndex; xRightBeg = xIndex; xRightLen = xLen - xRightBeg;
        yLeftBeg = 0; yLeftLen = yIndex; yRightBeg = yIndex; yRightLen = yLen - yRightBeg;
    }
    else
    {
        xLeftBeg = 0; xLeftLen = xIndex; xRightBeg = xIndex; xRightLen = xLen - xRightBeg;

        yLeftBeg = yIndex + 1; yLeftLen = yLen - yLeftBeg;
        yRightBeg = 0; yRightLen = yIndex + 1;
    }

    flag = Get_type(*hap_can);


    if(flag == XCY)
    {
        get_pair_hap_similarity(yReads->a, yLen, xUid, reverse_sources, read_g, ruIndex, 
        &yLeftMatch, &yLeftTotal);
        max_count = yLeftMatch;
        min_count = yLeftTotal;
    }
    else if(flag == YCX)
    {
        get_pair_hap_similarity(xReads->a, xLen, yUid, reverse_sources, read_g, ruIndex, 
        &xLeftMatch, &xLeftTotal);
        max_count = xLeftMatch;
        min_count = xLeftTotal;
    }
    else if(flag == X2Y)
    {
        get_pair_hap_similarity(yReads->a+yLeftBeg, yLeftLen, xUid, reverse_sources, read_g, ruIndex,
        &yLeftMatch, &yLeftTotal);
        get_pair_hap_similarity(xReads->a+xRightBeg, xRightLen, yUid, reverse_sources, read_g, ruIndex, 
        &xRightMatch, &xRightTotal);
        max_count = yLeftMatch + xRightMatch;
        min_count = yLeftTotal + xRightTotal;
    }
    else if(flag == Y2X)
    {
        get_pair_hap_similarity(xReads->a+xLeftBeg, xLeftLen, yUid, reverse_sources, read_g, ruIndex, 
        &xLeftMatch, &xLeftTotal);
        get_pair_hap_similarity(yReads->a+yRightBeg, yRightLen, xUid, reverse_sources, read_g, ruIndex,
        &yRightMatch, &yRightTotal);
        max_count = xLeftMatch + yRightMatch;
        min_count = xLeftTotal + yRightTotal;
    } else abort();

    hap_can->weight = hap_can->index_beg = 0;
    if(min_count == 0) return NON_PLOID;
    if((max_count > min_count*Hap_rate) || is_local)
    {
        long long r_x_interval_beg, r_x_interval_end, r_y_interval_beg, r_y_interval_end;
        uint64_t ploid_coverage = 0;

        ///for containment, don't need to do anything
        get_hap_alignment_boundary(xReads, yReads, flag, xLeftMatch, xLeftTotal, 
        yLeftMatch, yLeftTotal, xRightMatch, xRightTotal, yRightMatch, yRightTotal,
        xLeftBeg, xLeftLen, yLeftBeg, yLeftLen, xRightBeg, xRightLen, yRightBeg, yRightLen,
        xUid, yUid, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, 
        hap_can->rev, &r_x_interval_beg, &r_x_interval_end, &r_y_interval_beg, &r_y_interval_end);

        if(r_x_interval_beg < 0 || r_x_interval_end < 0 || r_y_interval_beg < 0 || r_y_interval_end < 0)
        {
            return NON_PLOID;
        }

        get_pair_hap_similarity(xReads->a + r_x_interval_beg, r_x_interval_end + 1 - r_x_interval_beg, 
        yUid, reverse_sources, read_g, ruIndex, &xLeftMatch, &xLeftTotal);
        if(xLeftMatch == 0 || xLeftTotal == 0 || (is_local == 0 && xLeftMatch <= xLeftTotal*Hap_rate))
        {
            return NON_PLOID;
        }
                
        
        hap_can->weight = xLeftMatch;
        hap_can->index_beg = xLeftTotal;
        hap_can->index_end = flag;
        hap_can->x_beg_pos = r_x_interval_beg;
        hap_can->x_end_pos = r_x_interval_end;
        hap_can->y_beg_pos = r_y_interval_beg;
        hap_can->y_end_pos = r_y_interval_end;
        hap_can->index_end = determine_hap_overlap_type_advance(hap_can, xReads, yReads,
        ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, min_ovlp, 
        xUid, yUid, u_buffer, tailIndex, prevIndex, r_x_pos_beg, r_x_pos_end, r_y_pos_beg, 
        r_y_pos_end);
        if(hap_can->index_end == XCY && yReads->len > (xReads->len*2)) return NON_PLOID;
        if(hap_can->index_end == YCX && xReads->len > (yReads->len*2)) return NON_PLOID;
        if(hap_can->index_end == (uint32_t)-1) return NON_PLOID;

        ploid_coverage = get_pair_purge_coverage(xReads, *r_x_pos_beg, *r_x_pos_end, 
        yReads, *r_y_pos_beg, *r_y_pos_end, hap_can->rev, read_g, cov);

        if(cov_threshold > 0 && ploid_coverage >= cov_threshold) return NON_PLOID;

        get_pair_hap_similarity_by_base(xReads, read_g, yUid, reverse_sources, ruIndex, 
        *r_x_pos_beg, *r_x_pos_end, &xLeftMatch, &xLeftTotal);
        if(xLeftMatch == 0 || xLeftTotal == 0 || xLeftMatch <= xLeftTotal*Hap_rate)
        {
            return NON_PLOID;
        }

        return PLOID;
    } 
	return NON_PLOID;
}


void print_hap_paf(ma_ug_t *ug, hap_overlaps* ovlp)
{
    fprintf(stderr, "utg%.6d%c\t%u(%u)\t%u(%u)\t%u(%u)\t%c\tutg%.6d%c\t%u(%u)\t%u(%u)\t%u(%u)\t%u\t%u\t%lld(%u)\n", 
    ovlp->xUid+1, "lc"[ug->u.a[ovlp->xUid].circ], ug->u.a[ovlp->xUid].len, ug->u.a[ovlp->xUid].n,
    ovlp->x_beg_pos, ovlp->x_beg_id, ovlp->x_end_pos, ovlp->x_end_id, "+-"[ovlp->rev], 
    ovlp->yUid+1, "lc"[ug->u.a[ovlp->yUid].circ], ug->u.a[ovlp->yUid].len, ug->u.a[ovlp->yUid].n,
    ovlp->y_beg_pos, ovlp->y_beg_id, ovlp->y_end_pos, ovlp->y_end_id, ovlp->type, ovlp->weight, 
    ovlp->score, ovlp->status);
}

inline long long get_max_index(asg_arc_t_offset* x, int32_t* Scores, uint8_t* Flag, long long n,
long long x_readLen, long long y_readLen)
{
    long long i = 0, max_result = -1, max_i = -1, min_xLen = x_readLen * 2 + 2, x_off, y_off, tmp_xLen;
    for (i = 0; i < n; i++)
    {
        if(Flag[i] != 0) continue;
        x_off = (long long)(x[i].Off>>32);
        y_off = (long long)((uint32_t)x[i].Off);
        if(Scores[i] > max_result)
        {
            max_result = Scores[i];
            max_i = i;
            min_xLen = get_hap_overlapLen(x_off, x_off, x_readLen, y_off, y_off, y_readLen, 
            NULL, NULL, NULL, NULL);
        }
        else if(Scores[i] == max_result)
        {
            tmp_xLen = get_hap_overlapLen(x_off, x_off, x_readLen, y_off, y_off, y_readLen, 
            NULL, NULL, NULL, NULL);

            if(tmp_xLen < min_xLen)
            {
                max_result = Scores[i];
                max_i = i;
                min_xLen = tmp_xLen;
            }
        }
    }

    return max_i;
}


inline void get_chain_details(int32_t* Pres, int32_t* Results, uint8_t* Flag, long long max_i,
long long* chainLen, long long* dup)
{
    long long i = max_i;
    (*chainLen) = 0;
    (*dup) = 0;

    while (i >= 0)
    {
        if(Flag[i] == 1) (*dup)++;
        Results[(*chainLen)] = i;
        i = Pres[i];
        (*chainLen)++;
    }
}

inline void push_hap_can(asg_arc_t_offset* x, kvec_hap_candidates* u_can, int32_t* Results, 
long long chainLen, long long x_readLen, long long y_readLen)
{
    if(chainLen <= 0) return;
    hap_candidates hap_can;
    hap_can.rev = x[Results[0]].x.el;
    hap_can.x_beg_pos = hap_can.x_end_pos = (uint32_t)(x[Results[0]].Off>>32);
    hap_can.y_beg_pos = hap_can.y_end_pos = (uint32_t)(x[Results[0]].Off);
    hap_can.weight = 0;
    long long i = 0;
    uint64_t totalWeigth = 0;
    ///fprintf(stderr, "^^^chainLen: %lld\n", chainLen);
    for (i = 0; i < chainLen; i++)
    {
        ///fprintf(stderr, "i: %lld, Results[i]: %d\n", i, Results[i]);
        hap_can.x_beg_pos = (uint32_t)(x[Results[i]].Off>>32);
        hap_can.y_beg_pos = (uint32_t)(x[Results[i]].Off);
        hap_can.weight += x[Results[i]].weight;
    }

    for (i = 0; i < chainLen; i++)
    {
        totalWeigth += x[Results[i]].weight;
        if(totalWeigth >= (hap_can.weight/2)) break;
    }

    if(i >= chainLen) i = chainLen-1;
    ///Get_total(hap_can) = Results[i];
    hap_can.t = x[Results[i]].x;
    Get_type(hap_can) = classify_hap_overlap(hap_can.x_beg_pos, hap_can.x_end_pos, 
    x_readLen, hap_can.y_beg_pos, hap_can.y_end_pos, y_readLen, NULL, NULL, NULL, NULL);
    kv_push(hap_candidates, u_can->a, hap_can);
    if(hap_can.x_beg_pos > hap_can.x_end_pos || hap_can.y_beg_pos > hap_can.y_end_pos)
    {
        fprintf(stderr, "ERROR\n");
    }
}



void print_chain_data(int32_t* Scores, int32_t* Pres, int32_t* Begs, long long n)
{
    fprintf(stderr,"*****\nn_chain: %lld\n", n);
    long long i;
    for (i = 0; i < n; i++)
    {
        fprintf(stderr, "i: %lld, Scores: %d, Pres: %d, Begs: %d\n", 
        i, Scores[i], Pres[i], Begs[i]);
    }
}



void hap_chaining(asg_arc_t_offset* x, uint32_t n, kvec_t_i32_warp* score_vc, kvec_t_i32_warp* prevIndex_vec, 
kvec_t_i32_warp* begIndex_vec, kvec_t_u8_warp* flag_vec, float band_width_threshold, long long max_skip, 
long long x_readLen, long long y_readLen, kvec_hap_candidates* u_can)
{
    #define DUP_OVLP_RATE 0.75
    score_vc->a.n = prevIndex_vec->a.n = begIndex_vec->a.n = flag_vec->a.n = 0;
    if(n == 0) return;
    kv_resize(int32_t, score_vc->a, n);
    kv_resize(int32_t, prevIndex_vec->a, n);
    kv_resize(int32_t, begIndex_vec->a, n);
    kv_resize(uint8_t, flag_vec->a, n);
    
    int32_t* Scores = score_vc->a.a;
    int32_t* Pres = prevIndex_vec->a.a;
    int32_t* Begs = begIndex_vec->a.a;
    uint8_t* Flag = flag_vec->a.a;
    long long i, j, n_max_skip, x_off, y_off, max_beg, max_j = -1, max_score, score;
    long long distance_x, distance_y, total_distance_x, total_distance_y, distance_gap;
    float gap_rate, band_width_penalty = 1 / band_width_threshold;
    long long max_result, max_i, min_xLen, tmp_xLen, chainLen = 0, dup = 0;
    max_result = max_i = -1; min_xLen = x_readLen * 2 + 2;
    for (i = 0; i < n; i++)
    {
		n_max_skip = 0;

        x_off = (long long)(x[i].Off>>32);
        y_off = (long long)((uint32_t)x[i].Off);
        max_j = -1;
        max_score = x[i].weight;
        max_beg = i; //i itself
        ///may have a pre-cut condition for j
        for (j = i - 1; j >= 0; --j) 
        {
            distance_x = x_off - (long long)(x[j].Off>>32);
            distance_y = y_off - (long long)((uint32_t)x[j].Off);
            ///x has been sorted by x_off
            if(distance_x <= 0 || distance_y <= 0) continue;

            total_distance_x = x_off - (long long)(x[Begs[j]].Off>>32);
            total_distance_y = y_off - (long long)((uint32_t)x[Begs[j]].Off);

            distance_gap = total_distance_x - total_distance_y;
            if(distance_gap < 0) distance_gap = -distance_gap;
            if(distance_gap > band_width_threshold * total_distance_x)
            {
                continue;
            }

            score = x[i].weight;
            gap_rate = (float)((float)(distance_gap)/(float)(total_distance_x));
            score -= (long long)(score * gap_rate * band_width_penalty);
            score += Scores[j];

            ///find a new max score
			if (score > max_score) {
				max_score = score;
				max_j = j;
				max_beg = Begs[j];
				n_max_skip = 0;
			} 
            else 
            {
				if (++n_max_skip > max_skip) break;
			}
        }

        Scores[i] = max_score;
        Pres[i] = max_j;
        Begs[i] = max_beg;

        if(Scores[i] > max_result)
        {
            max_result = Scores[i];
            max_i = i;
            min_xLen = get_hap_overlapLen(x_off, x_off, x_readLen, y_off, y_off, y_readLen, 
            NULL, NULL, NULL, NULL);
        }
        else if(Scores[i] == max_result)
        {
            tmp_xLen = get_hap_overlapLen(x_off, x_off, x_readLen, y_off, y_off, y_readLen, 
            NULL, NULL, NULL, NULL);

            if(tmp_xLen < min_xLen)
            {
                max_result = Scores[i];
                max_i = i;
                min_xLen = tmp_xLen;
            }
        }
        Flag[i] = 0;
    }

    // print_asg_arc_t_offset(x, n);
    // print_chain_data(Scores, Pres, Begs, n);
    while (max_i != -1)
    {
        get_chain_details(Pres, Begs, Flag, max_i, &chainLen, &dup);
        if(chainLen == 0) break;
        if(dup > chainLen*DUP_OVLP_RATE)
        {
            for (i = 0; i < chainLen; i++)
            {
                if(Flag[Begs[i]] == 1) continue;
                Flag[Begs[i]] = 2;
            }
            
        }
        else
        {
            push_hap_can(x, u_can, Begs, chainLen, x_readLen, y_readLen);

            for (i = 0; i < chainLen; i++)
            {
                Flag[Begs[i]] = 1;
            }
        }

        max_i = get_max_index(x, Scores, Flag, n, x_readLen, y_readLen);
    }
}

void get_candidate_hap_alignment(kvec_hap_candidates* u_can, kvec_asg_arc_t_offset* u_buffer,
kvec_t_i32_warp* score_vc, kvec_t_i32_warp* prevIndex_vec, kvec_t_i32_warp* begIndex_vec,
kvec_t_u8_warp* flag_vec, float band_width_threshold, long long max_skip, long long x_readLen, 
long long y_readLen)
{
    u_can->a.n = 0;
    if(u_buffer->a.n == 0) return;
    uint32_t i = 0, anchor_i = 0, m = 1, break_point = (uint32_t)-1, is_merge;

    qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment_chaining);

    ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n);

    for (i = 1; i < u_buffer->a.n; i++)
    {
        is_merge = 0;
        if(u_buffer->a.a[m-1].x.el == u_buffer->a.a[i].x.el)
        {
            if(u_buffer->a.a[m-1].Off == u_buffer->a.a[i].Off) is_merge = 1;
            if(is_merge == 0 && (Get_xOff(u_buffer->a.a[m-1].Off)==Get_xOff(u_buffer->a.a[i].Off)))
            {
                if((Get_yOff(u_buffer->a.a[i].Off)-(Get_yOff(u_buffer->a.a[m-1].Off))) ==
                  (i-anchor_i))///not sure why, does it use for tolerate indels in overlaps? 
                {
                    is_merge = 1;
                }
            }

            if(is_merge)
            {
                u_buffer->a.a[m-1].weight += u_buffer->a.a[i].weight;
                continue;
            }
        } 
        u_buffer->a.a[m] = u_buffer->a.a[i];
        anchor_i = i;
        if(u_buffer->a.a[m].x.el != u_buffer->a.a[m-1].x.el) break_point = m;
        m++;
    }
    u_buffer->a.n = m;
    if(break_point > u_buffer->a.n) break_point = u_buffer->a.n;

    ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n);

    hap_chaining(u_buffer->a.a, break_point, score_vc, prevIndex_vec, begIndex_vec, flag_vec,
    band_width_threshold, max_skip, x_readLen, y_readLen, u_can);
    
    hap_chaining(u_buffer->a.a + break_point, u_buffer->a.n - break_point, score_vc, prevIndex_vec, 
    begIndex_vec, flag_vec, band_width_threshold, max_skip, x_readLen, y_readLen, u_can);
}

int filter_secondary_chain(long long max_score, long long cur_score, double rate)
{
    if(cur_score >= max_score) return 1;
    long long diff = max_score - cur_score;
    if(max_score < 0) max_score *= -1;
    if(diff >= max_score*(1-rate)) return 0;
    return 1;
}

static void hap_alignment_advance_worker(void *_data, long eid, int tid)
{
    hap_alignment_struct_pip* hap_buf = (hap_alignment_struct_pip*)_data;
    ma_ug_t *ug = hap_buf->ug;
    asg_t *read_g = hap_buf->read_g;
    ma_hit_t_alloc* sources = hap_buf->sources;
    ma_hit_t_alloc* reverse_sources = hap_buf->reverse_sources;
    R_to_U* ruIndex = hap_buf->ruIndex;
    ma_sub_t *coverage_cut = hap_buf->coverage_cut;
    uint64_t* position_index = hap_buf->position_index;
    float Hap_rate = hap_buf->Hap_rate;
    int max_hang = hap_buf->max_hang;
    int min_ovlp = hap_buf->min_ovlp;
    float chain_rate = hap_buf->chain_rate;
    hap_overlaps_list* all_ovlp = hap_buf->all_ovlp;
    uint32_t Input_uId = eid;
    uint64_t* vote_counting = hap_buf->buf[tid].vote_counting;
    uint8_t* visit = hap_buf->buf[tid].visit;
    kvec_t_u64_warp* u_vecs = &(hap_buf->buf[tid].u_vecs);
    kvec_asg_arc_t_offset* u_buffer = &(hap_buf->buf[tid].u_buffer);
    kvec_hap_candidates* u_can = &(hap_buf->buf[tid].u_can);
    kvec_t_i32_warp* score_vc = &(hap_buf->buf[tid].u_buffer_tailIndex);
    kvec_t_i32_warp* prevIndex_vec = &(hap_buf->buf[tid].u_buffer_prevIndex);
    kvec_t_i32_warp* begIndex_vec = &(hap_buf->buf[tid].u_buffer_beg);
    kvec_t_u8_warp* flag_vec = &(hap_buf->buf[tid].u_buffer_flag);
    uint64_t cov_threshold = hap_buf->cov_threshold;
    hap_cov_t *cov = hap_buf->cov;
    if(hap_buf->cov_threshold < 0) cov_threshold = (uint64_t)-1;
    ma_utg_t *xReads = NULL, *yReads = NULL;
    ma_hit_t_alloc *xR = NULL;
    ma_hit_t *h = NULL;
    ma_sub_t *sq = NULL, *st = NULL;
    asg_t* nsg = ug->g;
    uint32_t i, j, v, rId, k, is_Unitig, Hap_uId, xUid, yUid, seedOcc;
    uint64_t tmp, max_weight, m;
    long long r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end, max_score;
    int32_t r;
    asg_arc_t t;
    asg_arc_t_offset t_offset;
    hap_overlaps hap_align;
    hap_overlaps *hap_align_x = NULL;
    xUid = Input_uId;
    if(nsg->seq[xUid].del || nsg->seq[xUid].c == ALTER_LABLE) return;
    memset(vote_counting, 0, sizeof(uint64_t)*nsg->n_seq);
    memset(visit, 0, nsg->n_seq);
    u_vecs->a.n = 0;
    u_can->a.n = 0;

    xReads = &(ug->u.a[xUid]);
    for (i = 0; i < xReads->n; i++)
    {
        xR = &(reverse_sources[xReads->a[i]>>33]);
        
        for (k = 0; k < xR->length; k++)
        {
            rId = Get_tn(xR->buffer[k]);
        
            if(read_g->seq[rId].del == 1)
            {
                ///get the id of read that contains it 
                get_R_to_U(ruIndex, rId, &rId, &is_Unitig);
                if(rId == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[rId].del == 1) continue;
            }

            ///there are two cases: 
            ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id  
            ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1
            get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig);
            if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
            ///here rId is the id of the read coming from the different haplotype
            ///Hap_cId is the id of the corresponding contig (note here is the contig, instead of untig)
            if(visit[Hap_uId]!=0) continue; ///one read only has one vote for one hap unitig
            visit[Hap_uId] = 1;
            if(vote_counting[Hap_uId] < UINT64_MAX) vote_counting[Hap_uId]++;
        }

        clean_visit_flag(visit, read_g, ruIndex, nsg->n_seq, xR);
    }



    u_vecs->a.n = 0;
    for (i = 0; i < nsg->n_seq; i++)
    {
        if(i == xUid) continue;
        if(vote_counting[i] == 0) continue;
        tmp = vote_counting[i]; tmp = tmp << 32; tmp = tmp | (uint64_t)i;
        kv_push(uint64_t, u_vecs->a, tmp);
    }

    if(u_vecs->a.n == 0) return;
    sort_kvec_t_u64_warp(u_vecs, 1);


    ///scan each candidate unitig
    for (i = 0; i < u_vecs->a.n; i++)
    {
        yUid = (uint32_t)u_vecs->a.a[i];
        seedOcc = u_vecs->a.a[i]>>32;
        xReads = &(ug->u.a[xUid]);
        yReads = &(ug->u.a[yUid]);
        u_buffer->a.n = 0;

        for (k = 0; k < xReads->n; k++)
        {
            xR = &(reverse_sources[xReads->a[k]>>33]);

            for (j = 0; j < xR->length; j++)
            {
                h = &(xR->buffer[j]);
                sq = &(coverage_cut[Get_qn(*h)]);
                st = &(coverage_cut[Get_tn(*h)]);
                if(st->del || read_g->seq[Get_tn(*h)].del) continue;

                r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, 
                                asm_opt.max_hang_rate, min_ovlp, &t);
                ///if it is a contained overlap, skip
                if(r < 0) continue;

                rId = t.v>>1;
                if(read_g->seq[rId].del == 1) continue;
                ///there are two cases: 
                ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id  
                ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1
                get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig);
                if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
                if(Hap_uId != yUid) continue;

                v = xReads->a[k]>>32;
                get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig);
                if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
                if(Hap_uId != xUid) continue;
                if((uint32_t)(position_index[v>>1]) != k) continue;

                if(asm_opt.purge_level_primary <= 2 &&
                    (prefilter((uint32_t)(position_index[v>>1]), (uint32_t)(position_index[rId]), 
                    xReads->n, yReads->n, 0, Hap_rate, seedOcc)==NON_PLOID) && 
                   (prefilter((uint32_t)(position_index[v>>1]), (uint32_t)(position_index[rId]), 
                    xReads->n, yReads->n, 1, Hap_rate, seedOcc)==NON_PLOID))
                {
                    continue;
                }
            
                t_offset.Off = get_xy_pos(read_g, &t, v, (yReads->a[(uint32_t)(position_index[rId])])>>32, 
                xReads->len, yReads->len, position_index, &(t.el));
                if(((t_offset.Off>>32) == (uint32_t)-1) || (((uint32_t)t_offset.Off) == (uint32_t)-1)) continue;

                t_offset.x = t;
                t_offset.weight = 1;

                kv_push(asg_arc_t_offset, u_buffer->a, t_offset);
            }

            deduplicate_edge(u_buffer);
        }

        if(u_buffer->a.n == 0) continue;

        
        get_candidate_hap_alignment(u_can, u_buffer, score_vc, prevIndex_vec, begIndex_vec,
        flag_vec, chain_rate, 50, xReads->len, yReads->len);

        if(u_can->a.n == 0) continue;
        
        qsort(u_can->a.a, u_can->a.n, sizeof(hap_candidates), cmp_hap_candidates);

        memset(&hap_align, 0, sizeof(hap_overlaps));    
        m = all_ovlp->x[xUid].a.n;     
        max_weight = 0; max_score = 0;
        for (k = 0; k < u_can->a.n; k++)
        {
            if(u_can->a.a[k].weight < max_weight*0.33) continue;
            if(calculate_pair_hap_similarity_advance(&(u_can->a.a[k]), position_index, xUid, yUid, 
            xReads, yReads, sources, reverse_sources, read_g, ruIndex, coverage_cut, Hap_rate,
            (asm_opt.purge_level_primary<=2? 0:1), max_hang, min_ovlp, cov_threshold, u_buffer, 
            score_vc, prevIndex_vec, cov, &r_x_pos_beg, &r_x_pos_end, &r_y_pos_beg, &r_y_pos_end)!=PLOID)
            {
                continue;
            }
            ///max_weight == 0 means the first matched chain
            if(max_weight == 0 || max_score < u_can->a.a[k].score) max_score = u_can->a.a[k].score;
            if(max_weight < u_can->a.a[k].weight) max_weight = u_can->a.a[k].weight;
            ///if one is positive and another one is negative, it is wrong
            if(!filter_secondary_chain(max_score, u_can->a.a[k].score, CHAIN_FILTER_RATE)) continue;
                            
            hap_align.rev = Get_rev(u_can->a.a[k]);
            hap_align.type = Get_type(u_can->a.a[k]);
            hap_align.x_beg_id = Get_x_beg(u_can->a.a[k]);
            hap_align.x_end_id = Get_x_end(u_can->a.a[k]) + 1;
            hap_align.y_beg_id = Get_y_beg(u_can->a.a[k]);
            hap_align.y_end_id = Get_y_end(u_can->a.a[k]) + 1;
            hap_align.weight = Get_match(u_can->a.a[k]);
            hap_align.score = u_can->a.a[k].score;
            hap_align.x_beg_pos = r_x_pos_beg;
            hap_align.x_end_pos = r_x_pos_end + 1;
            if(hap_align.rev == 0)
            {
                hap_align.y_beg_pos = r_y_pos_beg;
                hap_align.y_end_pos = r_y_pos_end + 1;
            }
            else
            {
                hap_align.y_beg_pos = yReads->len - r_y_pos_end - 1;
                hap_align.y_end_pos = yReads->len - r_y_pos_beg - 1 + 1;
            }
            hap_align.xUid = xUid;
            hap_align.yUid = yUid;
            hap_align.status = SELF_EXIST;
            kv_push(hap_overlaps, all_ovlp->x[hap_align.xUid].a, hap_align);
        }
        /**
        ///chains with same xUid && yUid
        for (k = m; k < all_ovlp->x[xUid].a.n; k++)
        {
            if(!filter_secondary_chain(max_score, 
                        all_ovlp->x[xUid].a.a[k].score, CHAIN_FILTER_RATE))
            {
                continue;
            } 
            all_ovlp->x[xUid].a.a[m] = all_ovlp->x[xUid].a.a[k];
            m++;
        }
        all_ovlp->x[xUid].a.n = m;
        **/
        
        hap_align_x = NULL;
        for (k = m; k < all_ovlp->x[xUid].a.n; k++)
        {
            if(all_ovlp->x[xUid].a.a[k].score != max_score) continue;
            if(hap_align_x == NULL || all_ovlp->x[xUid].a.a[k].weight > hap_align_x->weight)
            {
                hap_align_x = &(all_ovlp->x[xUid].a.a[k]);
            }
            else if(all_ovlp->x[xUid].a.a[k].weight == hap_align_x->weight)
            {
                if((all_ovlp->x[xUid].a.a[k].x_end_pos 
                            - all_ovlp->x[xUid].a.a[k].x_beg_pos) < 
                    (hap_align_x->x_end_pos - hap_align_x->x_beg_pos))
                {
                    hap_align_x = &(all_ovlp->x[xUid].a.a[k]);
                }
            }
        }
        if(hap_align_x)
        {
            all_ovlp->x[xUid].a.a[m] = (*hap_align_x);
            all_ovlp->x[xUid].a.n = m + 1;
        }
    }

}

int get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn)
{
    uint32_t i;
    for (i = 0; i < x->a.n; i++)
    {
        if(x->a.a[i].xUid == qn && x->a.a[i].yUid == tn)
        {
            return i;
        }        
    }

    return -1;
}

void set_reverse_hap_overlap(hap_overlaps* dest, hap_overlaps* source, uint32_t* types)
{
    dest->status = REVE_EXIST;
    dest->rev = source->rev;
    dest->type = types[source->type];
    dest->weight = source->weight;
    dest->xUid = source->yUid;
    dest->yUid = source->xUid;
    dest->x_beg_pos = source->y_beg_pos;
    dest->x_end_pos = source->y_end_pos;
    dest->y_beg_pos = source->x_beg_pos;
    dest->y_end_pos = source->x_end_pos;
    dest->x_beg_id = source->y_beg_id;
    dest->x_end_id = source->y_end_id;
    dest->y_beg_id = source->x_beg_id;
    dest->y_end_id = source->x_end_id;
    dest->score = source->score;
}

/**
#define X2Y 0
#define Y2X 1
#define XCY 2
#define YCX 3
**/
void normalize_hap_overlaps(hap_overlaps_list* all_ovlp, hap_overlaps_list* back_all_ovlp)
{
    hap_overlaps *x = NULL, *y = NULL;
    uint32_t v, i, uId, qn, tn;
    uint32_t types[4];
    types[X2Y] = Y2X; types[Y2X] = X2Y; types[XCY] = YCX; types[YCX] = XCY;
    int index;
    for (v = 0; v < all_ovlp->num; v++)
    {
        uId = v;
        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
        {
            qn = all_ovlp->x[uId].a.a[i].xUid;
            tn = all_ovlp->x[uId].a.a[i].yUid;
            x = &(all_ovlp->x[uId].a.a[i]);
            index = get_specific_hap_overlap(&(all_ovlp->x[tn]), tn, qn);
            if(index != -1)
            {
                y = &(all_ovlp->x[tn].a.a[index]);
                if(x->rev == y->rev && types[x->type]==y->type) continue;
                if(x->weight >= y->weight)
                {
                    kv_push(hap_overlaps, back_all_ovlp->x[tn].a, (*y));
                    set_reverse_hap_overlap(y, x, types);
                }
                else
                {
                    kv_push(hap_overlaps, back_all_ovlp->x[qn].a, (*x));
                    set_reverse_hap_overlap(x, y, types);
                }
            }
            else
            {
                kv_pushp(hap_overlaps, all_ovlp->x[tn].a, &y);
                set_reverse_hap_overlap(y, x, types);
            }
        }
    }
}

inline uint64_t calculate_bi_weight(hap_overlaps *x, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, 
R_to_U* ruIndex)
{
    double xMatch, xTotal;
    uint64_t weight = x->weight;
    ma_utg_t *yReads = &(ug->u.a[x->yUid]);
    get_pair_hap_similarity(yReads->a + x->y_beg_id, x->y_end_id - x->y_beg_id, 
    x->xUid, reverse_sources, read_g, ruIndex, &xMatch, &xTotal);
    weight += xMatch;
    return weight;
}

void normalize_hap_overlaps_advance(hap_overlaps_list* all_ovlp, hap_overlaps_list* back_all_ovlp,
ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex)
{
    hap_overlaps *x = NULL, *y = NULL;
    uint32_t v, i, uId, qn, tn;
    uint32_t types[4];
    
    types[X2Y] = Y2X; types[Y2X] = X2Y; types[XCY] = YCX; types[YCX] = XCY;
    int index;
    for (v = 0; v < all_ovlp->num; v++)
    {
        uId = v;
        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
        {
            qn = all_ovlp->x[uId].a.a[i].xUid;
            tn = all_ovlp->x[uId].a.a[i].yUid;
            x = &(all_ovlp->x[uId].a.a[i]);
            index = get_specific_hap_overlap(&(all_ovlp->x[tn]), tn, qn);
            if(index != -1)
            {
                y = &(all_ovlp->x[tn].a.a[index]);
                if(x->rev == y->rev && types[x->type]==y->type) continue;
                ///if(x->weight >= y->weight)
                // if((calculate_bi_weight(x, ug, read_g, reverse_sources, ruIndex)) >= 
                //    (calculate_bi_weight(y, ug, read_g, reverse_sources, ruIndex)))
                if(x->score >= y->score)
                {
                    kv_push(hap_overlaps, back_all_ovlp->x[tn].a, (*y));
                    set_reverse_hap_overlap(y, x, types);
                }
                else
                {
                    kv_push(hap_overlaps, back_all_ovlp->x[qn].a, (*x));
                    set_reverse_hap_overlap(x, y, types);
                }
            }
            else
            {
                kv_pushp(hap_overlaps, all_ovlp->x[tn].a, &y);
                set_reverse_hap_overlap(y, x, types);
            }
        }
    }
}

void get_p_nodes(p_g_t *pg, p_node_t **x, uint32_t *x_occ, uint32_t id)
{
    if(x) (*x) = pg->pg_het_node.a + pg->pg_h_lev_idx.a[id].beg;
    if(x_occ) (*x_occ) = pg->pg_h_lev_idx.a[id].occ;
}

void normalize_hap_overlaps_advance_by_p_g_t(hap_overlaps_list* all_ovlp, hap_overlaps_list* back_all_ovlp,
ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, p_g_t *pg, hap_cov_t *cov, 
double filter_rate)
{
    hap_overlaps *x = NULL, *y = NULL;
    uint32_t v, i, uId, qn, tn;
    uint32_t types[4];
    

    types[X2Y] = Y2X; types[Y2X] = X2Y; types[XCY] = YCX; types[YCX] = XCY;
    int index;
    uint32_t k, qs, qe, ts, te, occ, as, ae, ovlp, hetLen, homLen;
    p_node_t *a = NULL;


    for (v = 0; v < all_ovlp->num; v++)
    {
        uId = v;
        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
        {
            /*****************qn*****************/
            qn = all_ovlp->x[uId].a.a[i].xUid;
            qs = all_ovlp->x[uId].a.a[i].x_beg_pos;
            qe = all_ovlp->x[uId].a.a[i].x_end_pos - 1;
            get_p_nodes(pg, &a, &occ, qn);
            for (k = 0, hetLen = 0, homLen = 0; k < occ; k++)
            {
                as = a[k].baseBeg;
                ae = a[k].baseEnd;
                ovlp = ((MIN(qe, ae) >= MAX(qs, as))? MIN(qe, ae) - MAX(qs, as) + 1 : 0);
                if(homLen + hetLen > 0 && ovlp == 0) break;
                if(ovlp == 0) continue;
                if(a[k].h_status == N_HET)
                {
                    homLen += ovlp;
                } 
                else if(asm_opt.polyploidy <= 2 && (a[k].h_status&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].h_status&S_HET))
                {
                    homLen += ovlp;
                }
                else
                {
                    hetLen += ovlp;
                }                
            }

            if(hetLen <= ((hetLen + homLen) * filter_rate))
            {
                all_ovlp->x[uId].a.a[i].status = DELETE;
                continue;
            }
            /*****************qn*****************/


            /*****************tn*****************/
            tn = all_ovlp->x[uId].a.a[i].yUid;
            ts = all_ovlp->x[uId].a.a[i].y_beg_pos;
            te = all_ovlp->x[uId].a.a[i].y_end_pos - 1;
            get_p_nodes(pg, &a, &occ, tn);
            for (k = 0, hetLen = 0, homLen = 0; k < occ; k++)
            {
                as = a[k].baseBeg;
                ae = a[k].baseEnd;
                ovlp = ((MIN(te, ae) >= MAX(ts, as))? MIN(te, ae) - MAX(ts, as) + 1 : 0);
                if(homLen + hetLen > 0 && ovlp == 0) break;
                if(ovlp == 0) continue;
                if(a[k].h_status == N_HET)
                {
                    homLen += ovlp;
                } 
                else if(asm_opt.polyploidy <= 2 && (a[k].h_status&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].h_status&S_HET))
                {
                    homLen += ovlp;
                }
                else
                {
                    hetLen += ovlp;
                }                 
            }

            if(hetLen <= ((hetLen + homLen) * filter_rate))
            {
                all_ovlp->x[uId].a.a[i].status = DELETE;
                continue;
            }
            /*****************tn*****************/
        }
    }

    for (v = 0; v < all_ovlp->num; v++)
    {
        uId = v;
        k = 0;
        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
        {
            if(all_ovlp->x[uId].a.a[i].status == DELETE) continue;
            all_ovlp->x[uId].a.a[k] = all_ovlp->x[uId].a.a[i];
            k++;
        }
        all_ovlp->x[uId].a.n = k;
    }


    for (v = 0; v < all_ovlp->num; v++)
    {
        uId = v;
        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
        {
            qn = all_ovlp->x[uId].a.a[i].xUid;
            tn = all_ovlp->x[uId].a.a[i].yUid;
            x = &(all_ovlp->x[uId].a.a[i]);
            index = get_specific_hap_overlap(&(all_ovlp->x[tn]), tn, qn);
            if(index != -1)
            {
                y = &(all_ovlp->x[tn].a.a[index]);
                if(x->rev == y->rev && types[x->type]==y->type) continue;
                if(x->score >= y->score)
                {
                    kv_push(hap_overlaps, back_all_ovlp->x[tn].a, (*y));
                    set_reverse_hap_overlap(y, x, types);
                }
                else
                {
                    kv_push(hap_overlaps, back_all_ovlp->x[qn].a, (*x));
                    set_reverse_hap_overlap(x, y, types);
                }
            }
            else
            {
                kv_pushp(hap_overlaps, all_ovlp->x[tn].a, &y);
                set_reverse_hap_overlap(y, x, types);
            }
        }
    }
}

void filter_hap_overlaps_by_length(hap_overlaps_list* all_ovlp, uint32_t minLen)
{
    if(minLen == 0) return;
    hap_overlaps *x = NULL;
    uint32_t v, i, m, uId;    

    for (v = 0; v < all_ovlp->num; v++)
    {
        uId = v;
        m = 0;
        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
        {
            x = &(all_ovlp->x[uId].a.a[i]);
            if(x->x_end_id - x->x_beg_id < minLen) continue;
            all_ovlp->x[uId].a.a[m] = (*x);
            m++;
        }
        all_ovlp->x[uId].a.n = m;
    }
}

void debug_hap_overlaps(hap_overlaps_list* all_ovlp, hap_overlaps_list* back_all_ovlp)
{
    hap_overlaps *x = NULL, *y = NULL;
    uint32_t v, i, uId, qn, tn;
    uint32_t types[4];
    types[X2Y] = Y2X; types[Y2X] = X2Y; types[XCY] = YCX; types[YCX] = XCY;
    int index;
    for (v = 0; v < all_ovlp->num; v++)
    {
        uId = v;
        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
        {
            qn = all_ovlp->x[uId].a.a[i].xUid;
            tn = all_ovlp->x[uId].a.a[i].yUid;
            x = &(all_ovlp->x[uId].a.a[i]);
            index = get_specific_hap_overlap(&(all_ovlp->x[tn]), tn, qn);
            if(index == -1)
            {
                fprintf(stderr, "ERROR 0\n");
                continue;
            }

            y = &(all_ovlp->x[tn].a.a[index]);
            if(x->rev != y->rev || types[x->type] != y->type)
            {
                fprintf(stderr, "ERROR 1\n");
                continue;
            }

            if(x->status == REVE_EXIST && y->status != SELF_EXIST)
            {
                fprintf(stderr, "ERROR 2\n");
                continue;
            }

            if(x->status == REVE_EXIST)
            {
                if(x->weight != y->weight) fprintf(stderr, "ERROR 3\n"); 
                if(x->xUid != y->yUid) fprintf(stderr, "ERROR 4\n"); 
                if(x->yUid != y->xUid) fprintf(stderr, "ERROR 5\n"); 
                if(x->x_beg_pos != y->y_beg_pos) fprintf(stderr, "ERROR 6\n"); 
                if(x->x_end_pos != y->y_end_pos) fprintf(stderr, "ERROR 7\n"); 
                if(x->y_beg_pos != y->x_beg_pos) fprintf(stderr, "ERROR 8\n");
                if(x->y_end_pos != y->x_end_pos) fprintf(stderr, "ERROR 9\n");
                if(x->x_beg_id != y->y_beg_id) fprintf(stderr, "ERROR 10\n");
                if(x->x_end_id != y->y_end_id) fprintf(stderr, "ERROR 11\n");
                if(x->y_beg_id != y->x_beg_id) fprintf(stderr, "ERROR 12\n");
                if(x->y_beg_id != y->x_beg_id) fprintf(stderr, "ERROR 13\n");
                if(x->y_end_id != y->x_end_id) fprintf(stderr, "ERROR 14\n");

                index = get_specific_hap_overlap(&(back_all_ovlp->x[qn]), qn, tn);
                if(index != -1)
                {
                    if(back_all_ovlp->x[qn].a.a[index].weight > x->weight) fprintf(stderr, "ERROR 15\n");
                }
            }
           
        }
    }
}

void print_purge_gfa(ma_ug_t *ug, asg_t *purge_g)
{
    uint32_t v, i, n_vtx = purge_g->n_seq * 2;
    for (v = 0; v < n_vtx; v++)
    {
        if(v%2==0) fprintf(stderr, "\n");
        if(purge_g->seq[v>>1].del)
        {
            fprintf(stderr, "(D) v>>1: %u, v&1: %u, utg%.6d%c\n", v>>1, v&1, (v>>1)+1,
            "lc"[ug->u.a[v>>1].circ]);
            continue;
        } 

        fprintf(stderr, "(E) v>>1: %u, v&1: %u, utg%.6dl%c\n", v>>1, v&1, (v>>1)+1,
        "lc"[ug->u.a[v>>1].circ]);
        
        uint32_t nv = asg_arc_n(purge_g, v);
		asg_arc_t *av = asg_arc_a(purge_g, v);
        for (i = 0; i < nv; i++)
        {
            if(av[i].del) continue;
            fprintf(stderr, "av[i].ul: %u (utg%.6d%c, dir: %u, len: %u), av[i].v: %u (utg%.6d%c, dir: %u, len: %u), ol: %u\n", 
                (uint32_t)(av[i].ul>>33), (uint32_t)(av[i].ul>>33)+1, "lc"[ug->u.a[av[i].ul>>33].circ], (uint32_t)(av[i].ul>>32)&1, ug->u.a[av[i].ul>>33].len, 
                av[i].v>>1, (av[i].v>>1)+1, "lc"[ug->u.a[av[i].v>>1].circ], av[i].v&1, ug->u.a[av[i].v>>1].len, av[i].ol);
        }

    }
    
}

long long decode_score(uint32_t h_bits, uint32_t l_bits)
{
    uint64_t x; 
    x = h_bits; x <<= 32; x += l_bits;
    long long score = ((uint64_t)((uint64_t)x<<1)>>1);
    if((x>>63) == 0) score *= -1;
    return score;
}

void encode_score(long long i_s, uint32_t *h_bits, uint32_t *l_bits)
{
    uint64_t score = (i_s >= 0? (i_s) : (i_s*(-1)));
    if(i_s >= 0) score += (((uint64_t)1)<<63);
    (*l_bits) = (uint32_t)score; (*h_bits)= (score>>32);
}

uint64_t asg_bub_pop1_purge_graph(asg_t *g, uint32_t v0, int max_dist, buf_t *b)
{   
	uint32_t i, n_pending = 0, n_tips, tip_end;
	uint64_t n_pop = 0;
	///if this node has been deleted
	if (g->seq[v0>>1].del || g->seq[v0>>1].c == ALTER_LABLE) return 0; // already deleted
	///if ((uint32_t)g->idx[v0] < 2) return 0; // no bubbles
    if(get_real_length(g, v0, NULL)<2) return 0;
	///S saves nodes with all incoming edges visited
	b->S.n = b->T.n = b->b.n = b->e.n = 0;
	///for each node, b->a saves all related information
	b->a[v0].c = b->a[v0].d = b->a[v0].m = b->a[v0].nc = b->a[v0].np = 0;
	///b->S is the nodes with all incoming edges visited
	kv_push(uint32_t, b->S, v0);
    n_tips = 0;
    tip_end = (uint32_t)-1;

	do {
		///v is a node that all incoming edges have been visited
		///d is the distance from v0 to v
		uint32_t v = kv_pop(b->S), d = b->a[v].d;
		uint32_t nv = asg_arc_n(g, v);
		asg_arc_t *av = asg_arc_a(g, v);
        long long t_s = decode_score(b->a[v].c, b->a[v].m), c_s;
		///why we have this assert?
		///assert(nv > 0);
		///all out-edges of v
		for (i = 0; i < nv; ++i) { // loop through v's neighbors
            ///if this edge has been deleted
			if (av[i].del) continue;
			uint32_t w = av[i].v; // v->w with length l
			binfo_t *t = &b->a[w];
			///that means there is a circle, directly terminate the whole bubble poping
			///if (w == v0) goto pop_reset;
            if ((w>>1) == (v0>>1)) goto pop_reset;
            c_s = decode_score((uint32_t)av[i].ul, av[i].ol);
			///push the edge
            ///high 32-bit of g->idx[v] is the start point of v's edges
            //so here is the point of this specfic edge
			kv_push(uint32_t, b->e, (g->idx[v]>>32) + i);
			///find a too far path? directly terminate the whole bubble poping
			if (d + 1 > (uint32_t)max_dist) break; // too far

            ///if this node
			if (t->s == 0) { // this vertex has never been visited
				kv_push(uint32_t, b->b, w); // save it for revert
				///t->p is the parent node of 
				///t->s = 1 means w has been visited
				///d is len(v0->v), l is len(v->w), so t->d is len(v0->w)
				t->p = v, t->s = 1, t->d = d + 1;
                encode_score(t_s + c_s, &(t->c), &(t->m));
				///incoming edges of w
				///t->r = count_out(g, w^1);
                t->r = get_real_length(g, w^1, NULL);
				++n_pending;
			} else { // visited before
                if((t_s + c_s)> decode_score(t->c, t->m))
                {
                    t->p = v;
                    encode_score(t_s + c_s, &(t->c), &(t->m));
                }
                ///it is the shortest edge
				if (d + 1 < t->d) t->d = d + 1; // update dist
			}
			///assert(t->r > 0);
			//if all incoming edges of w have visited
			//push it to b->S
			if (--(t->r) == 0) {
                uint32_t x = get_real_length(g, w, NULL);
                /****************************may have bugs for bubble********************************/
                if(x > 0)
                {
                    kv_push(uint32_t, b->S, w);
                }
                else
                {
                    ///at most one tip
                    if(n_tips != 0) goto pop_reset;
                    n_tips++;
                    tip_end = w;
                }
                /****************************may have bugs for bubble********************************/
				--n_pending;
			}
		}
        //if found a tip
        /****************************may have bugs for bubble********************************/
        if(n_tips == 1)
        {
            if(tip_end != (uint32_t)-1 && n_pending == 0 && b->S.n == 0)
            {
                kv_push(uint32_t, b->S, tip_end);
                break;
            }
            else
            {
                goto pop_reset;
            }
        }
        /****************************may have bugs for bubble********************************/
		///if i < nv, that means (d + l > max_dist)
		if (i < nv || b->S.n == 0) goto pop_reset;
	} while (b->S.n > 1 || n_pending);

	asg_bub_backtrack_primary(g, v0, b);

    n_pop = 1;
pop_reset:
	for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices
		binfo_t *t = &b->a[b->b.a[i]];
		t->s = t->c = t->d = t->m = t->nc = t->np = 0;
	}
	return n_pop;
}



// pop bubbles
int asg_pop_bubble_purge_graph(asg_t *purge_g)
{
	uint32_t v, n_vtx = purge_g->n_seq * 2;
	uint64_t n_pop = 0;
	buf_t b;
	if (!purge_g->is_symm) asg_symm(purge_g);
	memset(&b, 0, sizeof(buf_t));
	///set information for each node
	b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t));
	//traverse all node with two directions 
	for (v = 0; v < n_vtx; ++v) {
		uint32_t i, n_arc = 0, nv = asg_arc_n(purge_g, v);
		asg_arc_t *av = asg_arc_a(purge_g, v);
		///some node could be deleted
		if (nv < 2 || purge_g->seq[v>>1].del || purge_g->seq[v>>1].c == ALTER_LABLE) continue;
		///some edges could be deleted
		for (i = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs
			if (!av[i].del) ++n_arc;
		if (n_arc > 1)
			n_pop += asg_bub_pop1_purge_graph(purge_g, v, purge_g->n_seq, &b);            
	}
	free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a);
	if (n_pop) asg_cleanup(purge_g);

	if(VERBOSE >= 1)
    {
        fprintf(stderr, "[M::%s] popped %lu bubbles\n", __func__, (unsigned long)n_pop);
    }
    return n_pop;
}

int get_hap_arch(hap_overlaps* hap, uint32_t qLen, uint32_t tLen, int max_hang, float max_hang_rate, 
int min_ovlp, asg_arc_t* t)
{
    int r;
    ma_hit_t h;
    h.qns = hap->xUid;
    h.qns = h.qns << 32;
    h.qns = h.qns | hap->x_beg_pos;
    h.qe = hap->x_end_pos;
    h.tn = hap->yUid;
    h.ts = hap->y_beg_pos;
    h.te = hap->y_end_pos;
    h.rev = hap->rev;
    h.del = 0;
    h.bl = h.el = h.ml = h.no_l_indel = 0;

    r = ma_hit2arc(&h, qLen, tLen, max_hang, max_hang_rate, min_ovlp, t);
    if(r < 0) return r;
    uint64_t score = (hap->score >= 0? (hap->score) : (hap->score*(-1)));
    if(hap->score >= 0) score += (((uint64_t)1)<<63);
    t->ol = (uint32_t)score;
    t->ul >>= 32; t->ul <<= 32; t->ul |= (score>>32);
    return r;
}

typedef struct {
    uint64_t eid;
    uint64_t score;
}e_score;

typedef struct {
    size_t n, m;
    e_score* a;
}e_score_warp;

#define e_score_key(a) ((a).score)
KRADIX_SORT_INIT(e_score, e_score, e_score_key, member_size(e_score, score))

int purge_g_arc_del_short_diploid_by_score(asg_t *g, float drop_ratio)
{
    e_score_warp b;
    kv_init(b);
    e_score *p = NULL;

	uint32_t v, n_vtx = g->n_seq * 2;
    long long n_cut = 0;
    
    for (v = 0; v < n_vtx; ++v) 
    {
        if(g->seq[v>>1].c == ALTER_LABLE || g->seq[v>>1].del) continue;
        asg_arc_t *av = asg_arc_a(g, v);
        uint32_t nv = asg_arc_n(g, v);
        if (nv < 2) continue;
        uint64_t i;
        for (i = 0; i < nv; ++i)
        {
            kv_pushp(e_score, b, &p);
            p->eid = av - g->arc + i;
            p->score = (uint32_t)av[i].ul; 
            p->score <<= 32;
            p->score |= av[i].ol;
        }
	}

    radix_sort_e_score(b.a, b.a + b.n);

    uint64_t k;
    for (k = 0; k < b.n; k++)
    {
        asg_arc_t *a = &g->arc[b.a[k].eid];
		///v is self id, w is the id of another end
		uint32_t i, v = (a->ul)>>32;
		uint32_t nv = asg_arc_n(g, v), kv;
		long long ovlp_max = 0, ovlp;
		asg_arc_t *av = NULL;
		///nv must be >= 2
		if (nv <= 1) continue;
        av = asg_arc_a(g, v);

        ///calculate the longest edge for v and w
		for (i = 0, kv = 0; i < nv; ++i) {
			if (av[i].del) continue;
            ovlp = decode_score((uint32_t)av[i].ul, av[i].ol);
            if (kv == 0 || ovlp_max < ovlp) ovlp_max = ovlp;
			++kv;
		}

        if (kv <= 1) continue;
        ovlp = decode_score((uint32_t)a->ul, a->ol);
		if (kv >= 2)
        {   
            if(ovlp >= 0 && ovlp_max >= 0 && ovlp > ovlp_max * drop_ratio) continue;
        }
         
        a->del = 1;
        asg_arc_del(g, a->v^1, av->ul>>32^1, 1);
        ++n_cut; 
    }

    kv_destroy(b); 
	if (n_cut) 
    {
		asg_cleanup(g);
		asg_symm(g);
	}

    
	return n_cut;
}


void clean_purge_graph(asg_t *purge_g, float drop_ratio, uint32_t is_force_break)
{
    uint64_t operation = 1;
    while (operation > 0)
    {
        operation = 0;
        operation += asg_pop_bubble_purge_graph(purge_g);
        operation += purge_g_arc_del_short_diploid_by_score(purge_g, drop_ratio);
    }
    
    if(is_force_break) purge_g_arc_del_short_diploid_by_score(purge_g, 1);
}


void get_node_boundary_advance(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, 
asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads,
uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, 
long long yEndIndex, uint32_t dir, uint32_t rev, kvec_asg_arc_t_offset* u_buffer, 
kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, asg_arc_t* reture_t_f, asg_arc_t* reture_t_r)
{
    long long k, j, offset, m;
    ma_hit_t_alloc *xR = NULL;
    ma_hit_t *h = NULL;
    ma_sub_t *sq = NULL, *st = NULL;
    int r, index;
    asg_arc_t t_f, t_r;
    uint32_t rId, Hap_uId, is_Unitig, v, w, v_dir, w_dir;
    uint64_t tmp;
    asg_arc_t_offset t_offset;
    reture_t_f->del = reture_t_r->del = 1;
    u_buffer->a.n = 0;


    for (k = xBegIndex; k <= xEndIndex; k++)
    {
        xR = &(reverse_sources[xReads->a[k]>>33]);
        for (j = 0; j < xR->length; j++)
        {
            h = &(xR->buffer[j]);
            sq = &(coverage_cut[Get_qn(*h)]);
            st = &(coverage_cut[Get_tn(*h)]);
            if(st->del || read_g->seq[Get_tn(*h)].del) continue;

            r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, 
                            asm_opt.max_hang_rate, min_ovlp, &t_f);
            ///if it is a contained overlap, skip
            if(r < 0) continue;

            rId = t_f.v>>1;
            if(read_g->seq[rId].del == 1) continue;
            ///there are two cases: 
            ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id  
            ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1
            get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig);
            if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
            if(Hap_uId != yUid) continue;

            v = xReads->a[k]>>32;
            get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig);
            if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
            if(Hap_uId != xUid) continue;
            if((uint32_t)(position_index[v>>1]) != k) continue;

            w = (yReads->a[(uint32_t)(position_index[rId])])>>32;

            v_dir = ((t_f.ul>>32)==v)?1:0;
            w_dir = (t_f.v == w)?1:0;
            if(rev == 0 && v_dir != w_dir) continue;
            if(rev == 1 && v_dir == w_dir) continue;
            if(dir == v_dir) continue;

            /****************************may have bugs********************************/
            offset = (uint32_t)(position_index[rId]);
            if(offset < yBegIndex || offset > yEndIndex) continue;
            /****************************may have bugs********************************/
        
            /************************get reverse edge*************************/
            index = get_specific_overlap(&(reverse_sources[Get_tn(*h)]), Get_tn(*h), Get_qn(*h));
            if(index == -1) continue;
            h = &(reverse_sources[Get_tn(*h)].buffer[index]);
            sq = &(coverage_cut[Get_qn(*h)]);
            st = &(coverage_cut[Get_tn(*h)]);
            if(st->del || read_g->seq[Get_tn(*h)].del) continue;
            r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, 
                            asm_opt.max_hang_rate, min_ovlp, &t_r);
            if(r < 0) continue;
            /************************get reverse edge*************************/

            tmp = get_xy_pos(read_g, &t_f, v, w, xReads->len, yReads->len, position_index, &(t_f.el));
            if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue;

            t_offset.Off = tmp;
            t_offset.x = t_f;
            t_offset.weight = 1;
            kv_push(asg_arc_t_offset, u_buffer->a, t_offset);

        }
    }

    if(u_buffer->a.n == 0) return;

    qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment_chaining);

    for (k = 1, m = 1; k < (long long)u_buffer->a.n; k++)
    {
        if(u_buffer->a.a[m-1].Off == u_buffer->a.a[k].Off)
        {
            u_buffer->a.a[m-1].weight += u_buffer->a.a[k].weight;
            if(u_buffer->a.a[k].x.ol > u_buffer->a.a[m-1].x.ol)
            {
                u_buffer->a.a[m-1].x = u_buffer->a.a[k].x;
            }
            continue;
        }
        u_buffer->a.a[m] = u_buffer->a.a[k];
        m++;
    }
    u_buffer->a.n = m;

    quick_LIS(u_buffer->a.a, u_buffer->a.n, tailIndex, prevIndex);

    if(tailIndex->a.n == 0) return;

    if(dir == 0)
    {
        for (k = 0; k < (long long)tailIndex->a.n; k++)
        {
            v = u_buffer->a.a[tailIndex->a.a[k]].x.v>>1;
            w = u_buffer->a.a[tailIndex->a.a[k]].x.ul>>33;
            index = get_specific_overlap(&(reverse_sources[v]), v, w);
            if(index == -1) continue;
            h = &(reverse_sources[v].buffer[index]);
            sq = &(coverage_cut[Get_qn(*h)]);
            st = &(coverage_cut[Get_tn(*h)]);
            if(st->del || read_g->seq[Get_tn(*h)].del) continue;
            r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, 
                            asm_opt.max_hang_rate, min_ovlp, &t_r);
            if(r < 0) continue;

            (*reture_t_f) = u_buffer->a.a[tailIndex->a.a[k]].x;
            (*reture_t_r) = t_r;
            return;
        }
    }
    else
    {
        for (k = tailIndex->a.n-1; k >= 0; k--)
        {
            v = u_buffer->a.a[tailIndex->a.a[k]].x.v>>1;
            w = u_buffer->a.a[tailIndex->a.a[k]].x.ul>>33;
            index = get_specific_overlap(&(reverse_sources[v]), v, w);
            if(index == -1) continue;
            h = &(reverse_sources[v].buffer[index]);
            sq = &(coverage_cut[Get_qn(*h)]);
            st = &(coverage_cut[Get_tn(*h)]);
            if(st->del || read_g->seq[Get_tn(*h)].del) continue;
            r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, 
                            asm_opt.max_hang_rate, min_ovlp, &t_r);
            if(r < 0) continue;

            (*reture_t_f) = u_buffer->a.a[tailIndex->a.a[k]].x;
            (*reture_t_r) = t_r;
            return;
        }
    }
}


void fill_unitig(uint64_t* buffer, uint32_t bufferLen, asg_t* read_g, kvec_asg_arc_t_warp* edge,
uint32_t is_circle, uint64_t* rLen)
{
    uint32_t i, k, totalLen, v, w, nv, l;
    asg_arc_t *av  = NULL;
    (*rLen) = totalLen = 0;
    for (i = 0; i < bufferLen - 1; i++)
    {
        v = (uint64_t)(buffer[i])>>32;
        w = (uint64_t)(buffer[i + 1])>>32; 
        av = asg_arc_a(read_g, v);
        nv = asg_arc_n(read_g, v);
        l = 0;
        
        for (k = 0; k < nv; k++)
        {
            if(av[k].del) continue;
            if(av[k].v == w) 
            {
                l = asg_arc_len(av[k]);
                break;
            }
        }

        if(k == nv)
        {
            for (k = 0; k < edge->a.n; k++)
            {
                if(edge->a.a[k].del) continue;
                if((edge->a.a[k].ul>>32) == v && edge->a.a[k].v == w)
                {
                    l = asg_arc_len(edge->a.a[k]);
                    break;
                }
            }

            if(k == edge->a.n)
            {
                fprintf(stderr, "####ERROR1-fill: i: %u, v>>1: %u, v&1: %u, w>>1: %u, w&1: %u\n",
                i, v>>1, v&1, w>>1, w&1);
            }
        } 


        
        buffer[i] = v; buffer[i] = buffer[i]<<32; buffer[i] = buffer[i] | (uint64_t)(l);
        totalLen += l;
    }
    
    if(i < bufferLen)
    {
        if(is_circle)
        {
            v = (uint64_t)(buffer[i])>>32;
            w = (uint64_t)(buffer[0])>>32; 
            av = asg_arc_a(read_g, v);
            nv = asg_arc_n(read_g, v);
            l = 0;
            
            for (k = 0; k < nv; k++)
            {
                if(av[k].del) continue;
                if(av[k].v == w) 
                {
                    l = asg_arc_len(av[k]);
                    break;
                }
            }

            if(k == nv)
            {
                for (k = 0; k < edge->a.n; k++)
                {
                    if(edge->a.a[k].del) continue;
                    if((edge->a.a[k].ul>>32) == v && edge->a.a[k].v == w)
                    {
                        l = asg_arc_len(edge->a.a[k]);
                        break;
                    }
                }

                if(k == edge->a.n)
                {
                    fprintf(stderr, "####ERROR2-fill: i: %u, v>>1: %u, v&1: %u, w>>1: %u, w&1: %u\n",
                    i, v>>1, v&1, w>>1, w&1);
                }
            } 

            buffer[i] = v; buffer[i] = buffer[i]<<32; buffer[i] = buffer[i] | (uint64_t)(l);
            totalLen += l;
        }
        else
        {
            v = (uint64_t)(buffer[i])>>32;
            l = read_g->seq[v>>1].len;
            buffer[i] = v;
            buffer[i] = buffer[i]<<32;
            buffer[i] = buffer[i] | (uint64_t)(l);
            totalLen += l;
        }
    }

    (*rLen) = totalLen;

}

void collect_trans_purge_cov(hap_cov_t *cov, ma_ug_t *ug, hap_overlaps* x, uint32_t is_keep_X)
{
    if(ug->u.a[x->xUid].n == 0 || ug->u.a[x->yUid].n == 0) return;
    uint64_t *pri = NULL, pri_n, *aux = NULL, aux_n, i, rId, uCov = 0, uLen = 0;

    if(is_keep_X)
    {
        pri = ug->u.a[x->xUid].a + x->x_beg_id;
        pri_n = x->x_end_id - x->x_beg_id;

        aux = ug->u.a[x->yUid].a + x->y_beg_id;
        aux_n = x->y_end_id - x->y_beg_id;
    }
    else
    {
        pri = ug->u.a[x->yUid].a + x->y_beg_id;
        pri_n = x->y_end_id - x->y_beg_id;

        aux = ug->u.a[x->xUid].a + x->x_beg_id;
        aux_n = x->x_end_id - x->x_beg_id;
    }


    uCov = uLen = 0;
    for (i = 0; i < aux_n; i++)
    {
        rId = aux[i]>>33;
        uCov += cov->cov[rId];
    }

    for (i = 0; i < pri_n; i++)
    {
        rId = pri[i]>>33;
        uLen += cov->read_g->seq[rId].len;
    }
    
    uCov = (uLen == 0? 0 : uCov / uLen);


    for (i = 0; i < pri_n; i++)
    {
        rId = pri[i]>>33;
        cov->cov[rId] += (uCov * cov->read_g->seq[rId].len);
    }
}


void collect_trans_purge_joint_cov(hap_cov_t *cov, ma_ug_t *ug, hap_overlaps* x)
{
    if(ug->u.a[x->xUid].n == 0 || ug->u.a[x->yUid].n == 0) return;
    uint64_t *a[2], a_n[2], uCov[2], uLen[2], uDepth[2], i, rId;
    
    a[0] = ug->u.a[x->xUid].a + x->x_beg_id;
    a_n[0] = x->x_end_id - x->x_beg_id;
    uCov[0] = uLen[0] = 0;
    for (i = 0; i < a_n[0]; i++)
    {
        rId = a[0][i]>>33;
        uCov[0] += cov->cov[rId];
        uLen[0] += cov->read_g->seq[rId].len;
    }
    
    a[1] = ug->u.a[x->yUid].a + x->y_beg_id;
    a_n[1] = x->y_end_id - x->y_beg_id;
    uCov[1] = uLen[1] = 0;
    for (i = 0; i < a_n[1]; i++)
    {
        rId = a[1][i]>>33;
        uCov[1] += cov->cov[rId];
        uLen[1] += cov->read_g->seq[rId].len;
    }
    
    uDepth[0] = (uLen[0] == 0? 0 : uCov[1] / uLen[0]);
    uDepth[1] = (uLen[1] == 0? 0 : uCov[0] / uLen[1]);

    for (i = 0; i < a_n[0]; i++)
    {
        rId = a[0][i]>>33;
        cov->cov[rId] += (uDepth[0] * cov->read_g->seq[rId].len);
    }

    for (i = 0; i < a_n[1]; i++)
    {
        rId = a[1][i]>>33;
        cov->cov[rId] += (uDepth[1] * cov->read_g->seq[rId].len);
    }
}



void purge_merge(asg_t *purge_g, ma_ug_t *ug, hap_overlaps_list* all_ovlp, buf_t* b_0,
R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, 
uint64_t* position_index, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, 
kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit,
hap_cov_t *cov)
{
    uint32_t i, nv, k, v, w, x_beg_index, x_end_index, y_beg_index, y_end_index, cut_beg, cut_end, begIndex, endIndex, keepUid;
    hap_overlaps *x = NULL/**, *y = NULL**/;
    ma_utg_t *xReads = NULL, *yReads = NULL;
    asg_arc_t t_forward, t_backward;
    asg_arc_t *av = NULL;
    kvec_t(uint64_t) buffer;
    uint64_t totalLen;
    int index = 0;
    i = 0;
    while (i < b_0->b.n)
    {
        cut_beg = 0; cut_end = (uint32_t)-1;
        kv_init(buffer);
        /********************for the first node********************/
        v = b_0->b.a[i];
        keepUid = v>>1;
        xReads = &(ug->u.a[v>>1]);
        if(v&1)
        {
            for (k = 0; k < xReads->n; k++)
            {
                kv_push(uint64_t, buffer, (xReads->a[xReads->n - k - 1])^(uint64_t)(0x100000000));
            }
        }
        else
        {
            for (k = 0; k < xReads->n; k++)
            {
                kv_push(uint64_t, buffer, xReads->a[k]);
            }
        }
        cut_beg = 0; cut_end = xReads->n - 1;
        i++;
        /********************for the first node********************/

        
        for (; i < b_0->b.n; i++)
        {
            ///x = y = NULL;
            x = NULL;

            v = b_0->b.a[i-1];
            w = b_0->b.a[i];
            

            index = get_specific_hap_overlap(&(all_ovlp->x[v>>1]), v>>1, w>>1);
            x = &(all_ovlp->x[v>>1].a.a[index]);


        
            xReads = &(ug->u.a[v>>1]);
            yReads = &(ug->u.a[w>>1]);

            begIndex = x->x_beg_id;
            if(cut_beg > begIndex) begIndex = cut_beg;

            endIndex = x->x_end_id-1;
            if(cut_end < endIndex) endIndex = cut_end;

            get_node_boundary_advance(ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, 
            min_ovlp, xReads, yReads, v>>1, w>>1, begIndex, endIndex, x->y_beg_id, x->y_end_id-1, v&1, 
            x->rev, u_buffer, tailIndex, prevIndex, &t_forward, &t_backward);

            if(t_forward.del || t_backward.del) break;
            
            kv_push(asg_arc_t, edge->a, t_forward);
            kv_push(asg_arc_t, edge->a, t_backward);

            x_beg_index = 0; x_end_index = xReads->n - 1;
            y_beg_index = 0; y_end_index = yReads->n - 1;
            
            if((v&1) == 0)
            {
                x_end_index = (uint32_t)position_index[t_forward.ul>>33];
                buffer.n = buffer.n - (cut_end - x_end_index);
            }
            else
            {
                x_beg_index = (uint32_t)position_index[t_forward.ul>>33];
                buffer.n = buffer.n - (x_beg_index - cut_beg);
            }

            if((w&1) == 1)
            {
                y_end_index = (uint32_t)position_index[t_forward.v>>1];
            }
            else
            {
                y_beg_index = (uint32_t)position_index[t_forward.v>>1];
            }

            cut_beg = y_beg_index;
            cut_end = y_end_index;

            if((w&1) == 1)
            {
                for (k = y_end_index; k >= y_beg_index; k--)
                {
                    kv_push(uint64_t, buffer, (yReads->a[k])^(uint64_t)(0x100000000));
                    if(k==0) break;
                }
            }
            else
            {
                for (k = y_beg_index; k <= y_end_index; k++)
                {
                    kv_push(uint64_t, buffer, yReads->a[k]);
                }
            }

            purge_g->seq[w>>1].c = ALTER_LABLE;
            if(cov) collect_trans_purge_joint_cov(cov, ug, x);

            // if(buffer.n > 1)
            // {
            //     for (k = 0; k < buffer.n - 1; k++)
            //     {
            //         if((buffer.a[k]>>32) == 854769 && (buffer.a[k+1]>>32) == 64486)
            //         {
            //             fprintf(stderr, "+++++++v: %u, w: %u, xReads->n: %u, yReads->n: %u\n", 
            //                         v, w, xReads->n, yReads->n);
            //             fprintf(stderr, "x->rev: %u, x->x_beg_id: %u, x->x_end_id: %u, x->y_beg_id: %u, x->y_end_id: %u\n", 
            //                         x->rev, x->x_beg_id, x->x_end_id, x->y_beg_id, x->y_end_id);
            //             fprintf(stderr, "t_forward.ul>>32: %u, t_forward.v: %u, y_beg_index: %u, y_end_index: %u\n", 
            //                         t_forward.ul>>32, t_forward.v, y_beg_index, y_end_index);
            //             fprintf(stderr, "type: %u, x->x_beg_pos: %u, x->x_end_pos: %u, xReads->len: %u\n", 
            //                         x->type, x->x_beg_pos, x->x_end_pos, xReads->len);
            //             fprintf(stderr, "x->y_beg_pos: %u, x->y_end_pos: %u, yReads->len: %u\n", 
            //                         x->y_beg_pos, x->y_end_pos, yReads->len);
            //         }
            //     }
            // }
        }

        // fprintf(stderr, "+keepUid: %u, i: %u, b_0->b.n: %u, buffer.n: %u\n", 
        //                             keepUid, i, (uint32_t)b_0->b.n, (uint32_t)buffer.n);
        fill_unitig(buffer.a, buffer.n, read_g, edge, 0, &totalLen);
        ///fprintf(stderr, "-keepUid: %u\n", keepUid);

        xReads = &(ug->u.a[keepUid]);
        free(xReads->a);
        xReads->a = buffer.a;
        xReads->n = buffer.n;
        xReads->m = buffer.m;
        xReads->len = totalLen;
        xReads->circ = 0;
        if(xReads->start != (xReads->a[0]>>32))
        {
            xReads->start = xReads->a[0]>>32;
            v = (keepUid<<1)+1;
            av = asg_arc_a(ug->g, v);
            nv = asg_arc_n(ug->g, v);
            for (k = 0; k < nv; k++)
            {
                if(av[k].del) continue;
                asg_arc_del(ug->g, av[k].ul>>32, av[k].v, 1);
                asg_arc_del(ug->g, av[k].v^1, av[k].ul>>32^1, 1);
            }
        }
        
        if(xReads->end != ((xReads->a[xReads->n-1]>>32)^1))
        {
            xReads->end = ((xReads->a[xReads->n-1]>>32)^1);
            v = (keepUid<<1);
            av = asg_arc_a(ug->g, v);
            nv = asg_arc_n(ug->g, v);
            for (k = 0; k < nv; k++)
            {
                if(av[k].del) continue;
                asg_arc_del(ug->g, av[k].ul>>32, av[k].v, 1);
                asg_arc_del(ug->g, av[k].v^1, av[k].ul>>32^1, 1);
            }
        }
    }

    for (i = 0; i < b_0->b.n; i++)
    {
        v = b_0->b.a[i];
        visit[v>>1] = 1;
        if(purge_g->seq[v>1].c != ALTER_LABLE) continue;
        asg_seq_drop(purge_g, v>1);
    }
}

void print_het_ovlp(p_g_t *pg, ma_ug_t *ug, hap_overlaps_list* ha, double filter_rate)
{
    uint32_t v, i, k, n_vtx = pg->pg_h_lev->n_seq * 2, nv, qn, qs, qe, tn, ts, te, as, ae, occ, ovlp, hetLen, homLen;
    asg_arc_t *av = NULL;
    hap_overlaps *x = NULL;
    p_node_t *a = NULL;
    int index;
    for (v = 0; v < n_vtx; v++)
    {
        av = asg_arc_a(pg->pg_h_lev, v);
        nv = asg_arc_n(pg->pg_h_lev, v);
        for (i = 0; i < nv; i++)
        {
            if(av[i].del) continue;
            index = get_specific_hap_overlap(&(ha->x[av[i].ul>>33]), av[i].ul>>33, av[i].v>>1);
            if(index == -1) fprintf(stderr, "ERROR\n");
            x = &(ha->x[av[i].ul>>33].a.a[index]);

            qn = x->xUid;
            qs = x->x_beg_pos;
            qe = x->x_end_pos - 1;
            get_p_nodes(pg, &a, &occ, qn);
            for (k = 0, hetLen = 0, homLen = 0; k < occ; k++)
            {
                as = a[k].baseBeg;
                ae = a[k].baseEnd;
                ovlp = ((MIN(qe, ae) >= MAX(qs, as))? MIN(qe, ae) - MAX(qs, as) + 1 : 0);
                if(homLen + hetLen > 0 && ovlp == 0) break;
                if(ovlp == 0) continue;
                if(a[k].h_status == N_HET)
                {
                    homLen += ovlp;
                } 
                else if(asm_opt.polyploidy <= 2 && (a[k].h_status&P_HET))
                {
                    homLen += ovlp;
                }
                else
                {
                    hetLen += ovlp;
                }  
            }

            if(hetLen <= ((hetLen + homLen) * filter_rate))
            {
                ///all_ovlp->x[uId].a.a[i].status = DELETE;
                fprintf(stderr, "********XY********\n");
                print_hap_paf(ug, x);
            }


            tn = x->yUid;
            ts = x->y_beg_pos;
            te = x->y_end_pos - 1;
            get_p_nodes(pg, &a, &occ, tn);
            for (k = 0, hetLen = 0, homLen = 0; k < occ; k++)
            {
                as = a[k].baseBeg;
                ae = a[k].baseEnd;
                ovlp = ((MIN(te, ae) >= MAX(ts, as))? MIN(te, ae) - MAX(ts, as) + 1 : 0);
                if(homLen + hetLen > 0 && ovlp == 0) break;
                if(ovlp == 0) continue;
                if(a[k].h_status == N_HET)
                {
                    homLen += ovlp;
                } 
                else if(asm_opt.polyploidy <= 2 && (a[k].h_status&P_HET))
                {
                    homLen += ovlp;
                }
                else
                {
                    hetLen += ovlp;
                }                 
            }

            if(hetLen <= ((hetLen + homLen) * filter_rate))
            {
                ///all_ovlp->x[uId].a.a[i].status = DELETE;
                fprintf(stderr, "********YX********\n");
                print_hap_paf(ug, x);
            }
        }
    }

    for (v = 0; v < ha->num; v++)
    {
        for (i = 0; i < ha->x[v].a.n; i++)
        {
            if(ha->x[v].a.a[i].status == DELETE)
            {
                x = &(ha->x[v].a.a[i]);

                qn = x->xUid;
                qs = x->x_beg_pos;
                qe = x->x_end_pos - 1;
                get_p_nodes(pg, &a, &occ, qn);
                for (k = 0, hetLen = 0, homLen = 0; k < occ; k++)
                {
                    as = a[k].baseBeg;
                    ae = a[k].baseEnd;
                    ovlp = ((MIN(qe, ae) >= MAX(qs, as))? MIN(qe, ae) - MAX(qs, as) + 1 : 0);
                    if(homLen + hetLen > 0 && ovlp == 0) break;
                    if(ovlp == 0) continue;
                    if(a[k].h_status == N_HET)
                    {
                        homLen += ovlp;
                    } 
                    else if(asm_opt.polyploidy <= 2 && (a[k].h_status&P_HET))
                    {
                        homLen += ovlp;
                    }
                    else
                    {
                        hetLen += ovlp;
                    }  
                }

                if(hetLen <= ((hetLen + homLen) * filter_rate))
                {
                    fprintf(stderr, "********C(X)********hetLen-%u, homLen-%u\n", hetLen, homLen);
                    print_hap_paf(ug, x);
                }


                tn = x->yUid;
                ts = x->y_beg_pos;
                te = x->y_end_pos - 1;
                get_p_nodes(pg, &a, &occ, tn);
                for (k = 0, hetLen = 0, homLen = 0; k < occ; k++)
                {
                    as = a[k].baseBeg;
                    ae = a[k].baseEnd;
                    ovlp = ((MIN(te, ae) >= MAX(ts, as))? MIN(te, ae) - MAX(ts, as) + 1 : 0);
                    if(homLen + hetLen > 0 && ovlp == 0) break;
                    if(ovlp == 0) continue;
                    if(a[k].h_status == N_HET)
                    {
                        homLen += ovlp;
                    } 
                    else if(asm_opt.polyploidy <= 2 && (a[k].h_status&P_HET))
                    {
                        homLen += ovlp;
                    }
                    else
                    {
                        hetLen += ovlp;
                    }                 
                }

                if(hetLen <= ((hetLen + homLen) * filter_rate))
                {
                    fprintf(stderr, "********C(Y)********hetLen-%u, homLen-%u\n", hetLen, homLen);
                    print_hap_paf(ug, x);
                }
            }
        }
    }
}

void link_unitigs(asg_t *purge_g, ma_ug_t *ug, hap_overlaps_list* all_ovlp,
R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, 
uint64_t* position_index, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, 
kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit,
hap_cov_t *cov)
{
    uint32_t v, n_vtx = purge_g->n_seq * 2, beg, end;
    long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen;
    buf_t b_0;
    memset(&b_0, 0, sizeof(buf_t));
    memset(visit, 0, purge_g->n_seq);
    for (v = 0; v < n_vtx; ++v) 
    {
        if(purge_g->seq[v>>1].c == ALTER_LABLE || purge_g->seq[v>>1].del || visit[v>>1]) continue;
        if(get_real_length(purge_g, v, NULL) != 1) continue;
        if(get_real_length(purge_g, v^1, NULL) != 0) continue;

        beg = v;
        b_0.b.n = 0;
        if(get_unitig(purge_g, NULL, beg, &end, &nodeLen, &baseLen, &max_stop_nodeLen, 
                    &max_stop_baseLen, 1, &b_0) == LOOP)
        {
            continue;
        }

        ///if(cov->link) collect_reverse_unitigs_purge(&b_0, cov->link, ug, all_ovlp);
        purge_merge(purge_g, ug, all_ovlp, &b_0, ruIndex, reverse_sources, coverage_cut, 
        read_g, position_index, u_buffer, tailIndex, prevIndex,max_hang, min_ovlp, edge, visit, cov);
    }
    free(b_0.b.a);
}

void print_all_purge_ovlp(ma_ug_t *ug, hap_overlaps_list* all_ovlp, const char* cmd)
{
    fprintf(stderr, "\n%s--->ug->u.n: %u\n", cmd, (uint32_t)ug->u.n);
    uint32_t v, uId, i;
    for (v = 0; v < all_ovlp->num; v++)
    {
        uId = v;
        ///if(uId != 96 && uId != 272) continue;
        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
        {
            print_hap_paf(ug, &(all_ovlp->x[uId].a.a[i]));
        }
    }

}

inline int get_available_cnt(asg_t *g, uint32_t v, uint8_t* del, asg_arc_t* v_s)
{
    //v has direction
    if(del && del[v>>1]) return 0;
    uint32_t i, kv = 0;
    asg_arc_t *av = asg_arc_a(g, v);
    uint32_t nv = asg_arc_n(g, v);

    for (i = 0, kv = 0; i < nv; i++)
    {
        if(!av[i].del)
        {
            if(del && del[av[i].v>>1]) continue;
            if(v_s) v_s[kv] = av[i];
            kv++;
        }
    }

    return kv;
}

long long get_specific_contig_length(asg_t *g, uint8_t *del)
{
    asg_cleanup(g);
    uint32_t v, n_vtx = g->n_seq * 2, q_occ;
    uint8_t *mark = NULL;
    ///is a queue
	//kdq_t(uint64_t) *q;
    ///each node has two directions
    //q = kdq_init(uint64_t);


    mark = (uint8_t*)calloc(n_vtx, 1);
    
    long long totalLen = 0;
    for (v = 0; v < n_vtx; ++v) 
    {
        uint32_t w, x, l, start, end, len;
        asg_arc_t arc;
        if (g->seq[v>>1].del || mark[v]) continue;
        if (get_available_cnt(g, v, del, NULL) == 0 && get_available_cnt(g, (v^1), del, NULL) != 0) continue;
        if (del[v>>1]) continue;

        mark[v] = 1;
		//q->count = 0, start = v, end = v^1, len = 0;
        q_occ =0, start = v, end = v^1, len = 0;
		// forward
		w = v;


        while (1) 
        {
			/**
			 * w----->x
			 * w<-----x
			 * that means the only suffix of w is x, and the only prefix of x is w
			 **/
			if (get_available_cnt(g, w, del, NULL) != 1) break;
            get_available_cnt(g, w, del, &arc);
			x = arc.v; // w->x
			if (get_available_cnt(g, x^1, del, NULL) != 1) break;

			/**
			 * another direction of w would be marked as used (since w has been used)
			**/
			mark[x] = mark[w^1] = 1;
			///l is the edge length, instead of overlap length
            ///note: edge length is different with overlap length
			///l = asg_arc_len(arc_first(g, w));
            get_available_cnt(g, w, del, &arc);
            l = ((uint32_t)((arc).ul));
			//kdq_push(uint64_t, q, (uint64_t)w<<32 | l);
            q_occ++;
			end = x^1, len += l;
			w = x;
			if (x == v) break;
		}


        //if (start != (end^1) || kdq_size(q) == 0) { // linear unitig
        if (start != (end^1) || q_occ == 0) { // linear unitig
			///length of seq, instead of edge
			l = g->seq[end>>1].len;
			//kdq_push(uint64_t, q, (uint64_t)(end^1)<<32 | l);
            q_occ++;
			len += l;
		} else { // circular unitig
			start = end = UINT32_MAX;
			goto add_unitig; // then it is not necessary to do the backward
		}

        // backward
		x = v;
		while (1) { // similar to forward but not the same
			if (get_available_cnt(g, x^1, del, NULL) != 1) break;
            get_available_cnt(g, x^1, del, &arc);
            w = arc.v ^ 1;
			if (get_available_cnt(g, w, del, NULL) != 1) break;
			mark[x] = mark[w^1] = 1;
			///l = asg_arc_len(arc_first(g, w));
            get_available_cnt(g, w, del, &arc);
            l = ((uint32_t)((arc).ul));
			///w is the seq id + direction, l is the length of edge
			///push element to the front of a queue
			//kdq_unshift(uint64_t, q, (uint64_t)w<<32 | l);
            q_occ++;

			start = w, len += l;
			x = w;
		}


        add_unitig:
        if (start != UINT32_MAX) mark[start] = mark[end] = 1;
        totalLen += len;
    }
    //kdq_destroy(uint64_t, q);
    free(mark);
    return totalLen;
}


void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen)
{
    uint8_t *del = (uint8_t *)malloc(sizeof(uint8_t)*g->n_seq);
    uint32_t v, k;
    ma_utg_t* u = NULL;
    memset(del, 1, g->n_seq);
    (*primaryLen) = (*alterLen) = 0;

    for (v = 0; v < ug->g->n_seq; ++v) 
    {
        if(ug->g->seq[v].del) continue;
        if(ug->g->seq[v].c == ALTER_LABLE) continue;
        u = &(ug->u.a[v]);
        if(u->m == 0) continue;
        for (k = 0; k < u->n; k++)
        {
            del[u->a[k]>>33] = 0;
        }
    }
    (*primaryLen) = get_specific_contig_length(g, del);


    for (v = 0; v < g->n_seq; ++v) 
    {
        del[v] = 1 - del[v];
    }

    (*alterLen) = get_specific_contig_length(g, del);

    free(del);
}


int if_ploid_sample(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, 
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, 
hap_alignment_struct_pip* hap_buf, hap_overlaps_list* all_ovlp, hap_overlaps_list* back_all_ovlp, 
uint32_t minLen, double purge_threshold)
{
    asg_t* nsg = ug->g;
    uint64_t v, k, total_bases = 0, alter_bases = 0, primary_bases = 0, purge_bases = 0;
    kt_for(asm_opt.thread_num, hap_alignment_advance_worker, hap_buf, nsg->n_seq);
    
    filter_hap_overlaps_by_length(all_ovlp, minLen);
    normalize_hap_overlaps_advance(all_ovlp, back_all_ovlp, ug, read_g, reverse_sources, ruIndex);

    get_contig_length(ug, read_g, &primary_bases, &alter_bases);
    total_bases = primary_bases + alter_bases;
    // fprintf(stderr, "primary_bases: %lu\n", primary_bases);
    // fprintf(stderr, "alter_bases: %lu\n", alter_bases);
    // fprintf(stderr, "total_bases: %lu\n", total_bases);


    for (v = 0; v < all_ovlp->num; v++)
    {
        for (k = 0; k < all_ovlp->x[v].a.n; k++)
        {
            purge_bases += all_ovlp->x[v].a.a[k].x_end_pos - all_ovlp->x[v].a.a[k].x_beg_pos;
        }
    }
    purge_bases = purge_bases/2;
    ///fprintf(stderr, "purge_bases: %lu\n", purge_bases);
    alter_bases = alter_bases + purge_bases;
    ///fprintf(stderr, "new alter_bases: %lu\n", alter_bases);


    for (v = 0; v < all_ovlp->num; v++)
    {
        all_ovlp->x[v].a.n = 0;
    }

    for (v = 0; v < back_all_ovlp->num; v++)
    {
        back_all_ovlp->x[v].a.n = 0;
    }

    if(alter_bases > total_bases * purge_threshold) return 1;
    return 0;
}

int cmp_chain_score(const void * a, const void * b)
{
    if((*(hap_overlaps*)a).score < (*(hap_overlaps*)b).score) return 1;
    if((*(hap_overlaps*)a).score > (*(hap_overlaps*)b).score) return -1;
    
    return 0;
}
long long get_ovlp_len(long long a_beg, long long a_end, long long b_beg, long long b_end)
{
    long long ovlp = (long long)(MIN(a_end, b_end)) - (long long)(MAX(a_beg, b_beg)) + 1;
    return ovlp <= 0? 0 : ovlp;
}
void sort_hap_chain(hap_overlaps_list* all_ovlp)
{
    hap_overlaps *x = NULL, *p = NULL;
    uint32_t v, i, k, uId;    
    long long ovlp, xLen, pLen;
    kvec_t(hap_overlaps) pri; kv_init(pri);
    kvec_t(hap_overlaps) alt; kv_init(alt);

    for (v = 0; v < all_ovlp->num; v++)
    {
        uId = v;
        qsort(all_ovlp->x[uId].a.a, all_ovlp->x[uId].a.n, sizeof(hap_overlaps), cmp_chain_score);
        pri.n = alt.n = 0;        
        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
        {
            x = &(all_ovlp->x[uId].a.a[i]);
            xLen = x->x_end_pos - x->x_beg_pos;
            for (k = 0; k < pri.n; k++)
            {
                p = &(pri.a[k]);
                pLen = p->x_end_pos - p->x_beg_pos;
                ovlp = get_ovlp_len(x->x_beg_pos, x->x_end_pos-1, p->x_beg_pos, p->x_end_pos-1);
                if(ovlp == 0) continue;
                if(ovlp >= (MIN(xLen, pLen))*0.5) break;
            }

            if(k < pri.n)
            {
                x->xUid = k;
                kv_push(hap_overlaps, alt, *x);
            }
            else
            {
                kv_push(hap_overlaps, pri, *x);
            } 
        }

    }

    kv_destroy(pri); kv_destroy(alt);
}

void remove_contained_haplotig(hap_overlaps_list* all_ovlp, ma_ug_t *ug, asg_t* nsg, asg_t *purge_g, hap_cov_t *cov)
{
    uint32_t v, i, uId, xUid;
    hap_overlaps *p = NULL;
    for (v = 0; v < all_ovlp->num; v++)
    {
        uId = v; p = NULL;
        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
        {
            if(p == NULL || p->score < all_ovlp->x[uId].a.a[i].score)
            {
                p = &(all_ovlp->x[uId].a.a[i]);
            }
        }

        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
        {
            if(all_ovlp->x[uId].a.a[i].type == YCX)
            {
                if(!filter_secondary_chain(p->score, all_ovlp->x[uId].a.a[i].score, 0.95))
                {
                    continue;
                } 
                
                xUid = all_ovlp->x[uId].a.a[i].xUid;

                nsg->seq[xUid].c = ALTER_LABLE;
                purge_g->seq[xUid].c = ALTER_LABLE;
                purge_g->seq[xUid].del = 1;

                all_ovlp->x[uId].a.a[i].status = DELETE;
                ///if(cov->link) collect_reverse_unitig_pair(cov->link, ug, &(all_ovlp->x[uId].a.a[i]));
                collect_trans_purge_cov(cov, ug, &(all_ovlp->x[uId].a.a[i]), 0);
            }

            ///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
        }
    }

    // for (v = 0; v < all_ovlp.num; v++)
    // {
    //     uId = v;
    //     for (i = 0; i < all_ovlp.x[uId].a.n; i++)
    //     {
    //         if(all_ovlp.x[uId].a.a[i].type == YCX)
    //         {
    //             nsg->seq[all_ovlp.x[uId].a.a[i].xUid].c = ALTER_LABLE;
    //             purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].c = ALTER_LABLE;
    //             purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].del = 1;
    //             all_ovlp.x[uId].a.a[i].status = DELETE;
    //             if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i]));
    //             collect_trans_purge_cov(cov, ug, &(all_ovlp.x[uId].a.a[i]), 0);
    //         }

    //         if(all_ovlp.x[uId].a.a[i].type == XCY)
    //         {
    //             nsg->seq[all_ovlp.x[uId].a.a[i].yUid].c = ALTER_LABLE;
    //             purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].c = ALTER_LABLE;
    //             purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].del = 1;
    //             all_ovlp.x[uId].a.a[i].status = DELETE;
    //             if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i]));
    //             collect_trans_purge_cov(cov, ug, &(all_ovlp.x[uId].a.a[i]), 1);
    //         }
    //         ///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
    //     }
    // }
}

void debug_p_g_t(p_g_t* pg, hap_cov_t *cov, asg_t *read_g)
{
    fprintf(stderr, "----------[M::%s]----------\n", __func__);
    uint32_t i, offset, v, sid, eid, spos, epos, p_status, p_uid, occ;
    p_node_t *t = NULL;
    ma_utg_t *u = NULL;
    p_node_t *a = NULL;

    for (v = 0; v < pg->ug->u.n; v++)
    {
        ///fprintf(stderr, "\nu->n: %u, uid: %u\n", (uint32_t)(pg->ug->u.a[v].n), v);
        get_p_nodes(pg, &a, &occ, v);
        for (i = 0; i < occ; i++)
        {
            if(a[i].c_ug_id != v) fprintf(stderr, "sbsbsbsbsb\n");
            ///fprintf(stderr, "sid: %u, eid: %u\n", a[i].nodeBeg, a[i].nodeEnd);
        }
    }
    

    for (v = 0, p_status = (uint32_t)-1, p_uid = (uint32_t)-1; v < pg->pg_het_node.n; v++)
    {
        t = &(pg->pg_het_node.a[v]);
        sid = t->nodeBeg;
        eid = t->nodeEnd;
        spos = t->baseBeg;
        epos = t->baseEnd;
        // fprintf(stderr, "sid: %u, eid: %u, spos: %u, epos: %u, t->b_ug_id: %u\n", 
        //                                                 sid, eid, spos, epos, (uint32_t)t->b_ug_id);
        // fprintf(stderr, "pg->pg_het_node.n: %u\n", (uint32_t)pg->pg_het_node.n);
        if(p_uid == t->c_ug_id && p_status == t->h_status)
        {
            fprintf(stderr, "ERROR-(-1)\n");
        } 
        p_status = t->h_status;
        p_uid = t->c_ug_id;

        u = &(pg->ug->u.a[t->c_ug_id]);
        ///fprintf(stderr, "u->n: %u, sid: %u, eid: %u\n", (uint32_t)u->n, sid, eid);
        for (i = offset = 0; i < u->n; i++)
        {
            if(i == sid)
            {
                if(spos != offset)
                {
                    fprintf(stderr, "ERROR-1\n");
                }
            }

            if(i == eid)
            {
                if(epos != (offset+read_g->seq[u->a[i]>>33].len - 1))
                {
                    fprintf(stderr, "ERROR-2, real end: %u\n", 
                                    (uint32_t)(offset+read_g->seq[u->a[i]>>33].len - 1));
                }
            }
            offset += (uint32_t)u->a[i];
            if(i >= sid && i <= eid)
            {
                if(cov->t_ch->is_r_het[u->a[i]>>33] != t->h_status)
                {
                    fprintf(stderr, "ERROR-(-3): is_r_het: %u, h_status: %u\n", cov->t_ch->is_r_het[u->a[i]>>33], t->h_status);
                }
            }
        }
    }
}


void print_p_g_t_interval(p_g_t* pg, hap_cov_t *cov)
{
    fprintf(stderr, "----------[M::%s]----------\n", __func__);
    uint32_t i, v, sid, eid;
    p_node_t *t = NULL;
    ma_utg_t *u = NULL;

    for (v = 0; v < pg->pg_het_node.n; v++)
    {
        t = &(pg->pg_het_node.a[v]);
        sid = t->nodeBeg;
        eid = t->nodeEnd;

        u = &(pg->ug->u.a[t->c_ug_id]);
        fprintf(stderr, "\nu->n=%u, sid=%u, eid=%u, h_status=%u\n", 
                                                (uint32_t)u->n, sid, eid, t->h_status);
        for (i = sid; i <= eid; i++)
        {
            fprintf(stderr, "id:i:%u------>utg%.6ul\n", 
                        (uint32_t)(u->a[i]>>33), (get_origin_uid(u->a[i]>>32, cov->t_ch, NULL, NULL)>>1)+1);
        }
    }
    fprintf(stderr, "----------[M::%s]----------\n", __func__);
}


p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g)
{
    uint32_t v, uId, k, l, offset, l_pos, g_beg_idx, occ/**, ovlp, tLen, zLen**/;
    p_g_t *pg = NULL; CALLOC(pg, 1);
    pg->ug = ug;
    asg_t* nsg = pg->ug->g;
    ma_utg_t *u = NULL;
    p_node_t *t = NULL/**, *z = NULL**/;
    ///asg_arc_t *e = NULL;
    p_g_in_t *x = NULL;
    ///pg->pg_het = asg_init();
    pg->pg_h_lev = asg_init();
    kv_init(pg->pg_het_node);
    kv_init(pg->pg_h_lev_idx);

    for (v = 0; v < nsg->n_seq; v++)
    {
        uId = v;
        if(nsg->seq[uId].del || nsg->seq[uId].c == ALTER_LABLE)
        {
            asg_seq_set(pg->pg_h_lev, uId, 0, 1);
            pg->pg_h_lev->seq[uId].c = ALTER_LABLE;
            continue;
        } 
        
        asg_seq_set(pg->pg_h_lev, uId, ug->u.a[uId].len, 0);
        pg->pg_h_lev->seq[uId].c = PRIMARY_LABLE;
    }

    // if(asm_opt.polyploidy <= 2)
    // {
    //     for (v = 0; v < cov->t_ch->r_num; v++)
    //     {
    //         if(cov->t_ch->is_r_het[v]&P_HET) cov->t_ch->is_r_het[v] |= S_HET;
    //     }
    // } 

    for (v = 0; v < nsg->n_seq; v++)
    {
        uId = v;
        if(nsg->seq[uId].del || nsg->seq[uId].c == ALTER_LABLE) continue;

        u = &(ug->u.a[uId]);
        g_beg_idx = pg->pg_het_node.n;
        ///fprintf(stderr, "\n+v: %u, pg->pg_het_node.n: %u\n", v, (uint32_t)pg->pg_het_node.n);
        for (k = 1, l = 0, offset = 0, l_pos = 0; k <= u->n; ++k) 
        {   
            ///if (k == u->n || (!!cov->t_ch->is_r_het[u->a[k]>>33]) != (!!cov->t_ch->is_r_het[u->a[l]>>33]))
            if (k == u->n || cov->t_ch->is_r_het[u->a[k]>>33] != cov->t_ch->is_r_het[u->a[l]>>33])
            {
                kv_pushp(p_node_t, pg->pg_het_node, &t);
                t->c_ug_id = uId;
                t->h_status = cov->t_ch->is_r_het[u->a[l]>>33];
                t->baseBeg = l_pos;
                t->baseEnd = offset + read_g->seq[u->a[k-1]>>33].len - 1;
                t->nodeBeg = l;
                t->nodeEnd = k - 1;
                ///if(t->b_ug_id == (uint32_t)-1) fprintf(stderr, "xxxx\n");
                ///asg_seq_set(pg->pg_het, pg->pg_het_node.n-1, t->baseEnd+1-t->baseBeg, 0);
                ///fprintf(stderr, "l: %u, k: %u, u->n: %u, t->h_status: %u\n", l, k, u->n, t->h_status);
                l = k;
                l_pos = offset + (uint32_t)u->a[k-1];
            }
            offset += (uint32_t)u->a[k-1];
        }

        occ = pg->pg_het_node.n - g_beg_idx;
        kv_pushp(p_g_in_t, pg->pg_h_lev_idx, &x);
        x->beg = g_beg_idx; x->occ = occ;
        ///fprintf(stderr, "-v: %u, pg->pg_het_node.n: %u\n", v, (uint32_t)pg->pg_het_node.n);
        
        // if(occ > 1)
        // {
        //     for (k = g_beg_idx; (k + 1) < pg->pg_het_node.n; ++k) 
        //     {
        //         t = &(pg->pg_het_node.a[k]); tLen = t->baseEnd + 1 - t->baseBeg;
        //         z = &(pg->pg_het_node.a[k+1]); zLen = z->baseEnd + 1 - z->baseBeg;

        //         ovlp = ((MIN(t->baseEnd, z->baseEnd) >= MAX(t->baseBeg, z->baseBeg))? 
        //                         MIN(t->baseEnd, z->baseEnd) - MAX(t->baseBeg, z->baseBeg) + 1 : 0);

        //         e = asg_arc_pushp(pg->pg_het);
        //         e->ol = ovlp;
        //         e->ul = (k<<1); e->ul <<= 32; e->ul += (tLen - ovlp);
        //         e->v = ((k+1)<<1); e->del = 0; e->el = e->no_l_indel = e->strong = 1;

        //         e = asg_arc_pushp(pg->pg_het);
        //         e->ol = ovlp;
        //         e->ul = ((k+1)<<1)+1; e->ul <<= 32; e->ul += (zLen - ovlp);
        //         e->v = (k<<1)+1; e->del = 0; e->el = e->no_l_indel = e->strong = 1;
        //     }
        // }
    }
    ///asg_cleanup(pg->pg_het);
    ///debug_p_g_t(pg, cov, read_g);
    ///print_p_g_t_interval(pg, cov);

    return pg;
}

void destory_p_g_t(p_g_t **pg)
{
    if(pg && (*pg))
    {
        kv_destroy((*pg)->pg_het_node);
        kv_destroy((*pg)->pg_h_lev_idx);
        asg_destroy((*pg)->pg_het);
        asg_destroy((*pg)->pg_h_lev);
        free((*pg)); 
        (*pg) = NULL;
    }
}


void chain_origin_trans_uid_by_purge(hap_overlaps *x, ma_ug_t *ug, hap_cov_t *cov, uint64_t* position_index)
{
    uint32_t pri_uid, aux_uid, r_x, r_y;
    hap_candidates hap_for, hap_rev, *hap = NULL;
    long long x_pos_beg, x_pos_end, y_pos_beg, y_pos_end;

    Get_rev(hap_for) = x->rev;
    Get_x_beg(hap_for) = x->x_beg_id; Get_x_end(hap_for) = x->x_end_id - 1;
    Get_y_beg(hap_for) = x->y_beg_id; Get_y_end(hap_for) = x->y_end_id - 1;
    r_x = determine_hap_overlap_type_advance(&hap_for, &(ug->u.a[x->xUid]), &(ug->u.a[x->yUid]),
    cov->ruIndex, cov->reverse_sources, cov->coverage_cut, cov->read_g, position_index, 
    cov->max_hang, cov->min_ovlp, x->xUid, x->yUid, &(cov->u_buffer), &(cov->tailIndex), 
    &(cov->prevIndex), &x_pos_beg, &x_pos_end, &y_pos_beg, &y_pos_end);

    Get_rev(hap_rev) = x->rev;
    Get_x_beg(hap_rev) = x->y_beg_id; Get_x_end(hap_rev) = x->y_end_id - 1;
    Get_y_beg(hap_rev) = x->x_beg_id; Get_y_end(hap_rev) = x->x_end_id - 1;
    r_y = determine_hap_overlap_type_advance(&hap_rev, &(ug->u.a[x->yUid]), &(ug->u.a[x->xUid]),
    cov->ruIndex, cov->reverse_sources, cov->coverage_cut, cov->read_g, position_index, 
    cov->max_hang, cov->min_ovlp, x->yUid, x->xUid, &(cov->u_buffer), &(cov->tailIndex), 
    &(cov->prevIndex), &x_pos_beg, &x_pos_end, &y_pos_beg, &y_pos_end);

    if(r_x == (uint32_t)-1 && r_y == (uint32_t)-1)
    {
        fprintf(stderr, "ERROR\n");
        return;
    }

    Get_rev(hap_for) = x->rev;
    Get_x_beg(hap_for) = x->x_beg_id; Get_x_end(hap_for) = x->x_end_id - 1;
    Get_y_beg(hap_for) = x->y_beg_id; Get_y_end(hap_for) = x->y_end_id - 1;

    Get_rev(hap_rev) = x->rev;
    Get_x_beg(hap_rev) = x->y_beg_id; Get_x_end(hap_rev) = x->y_end_id - 1;
    Get_y_beg(hap_rev) = x->x_beg_id; Get_y_end(hap_rev) = x->x_end_id - 1;


    if(r_x != (uint32_t)-1 && r_y == (uint32_t)-1)
    {
        pri_uid = x->xUid; aux_uid = x->yUid; hap = &hap_for;
    }
    else if(r_x == (uint32_t)-1 && r_y != (uint32_t)-1)
    {
        aux_uid = x->xUid; pri_uid = x->yUid; hap = &hap_rev;
    }
    else
    {
        if(hap_for.score >= hap_rev.score)
        {
            pri_uid = x->xUid; aux_uid = x->yUid; hap = &hap_for;
        }
        else
        {
            aux_uid = x->xUid; pri_uid = x->yUid; hap = &hap_rev;
        }
    }

    determine_hap_overlap_type_advance(hap, &(ug->u.a[pri_uid]), &(ug->u.a[aux_uid]),
    cov->ruIndex, cov->reverse_sources, cov->coverage_cut, cov->read_g, position_index, 
    cov->max_hang, cov->min_ovlp, pri_uid, aux_uid, &(cov->u_buffer), &(cov->tailIndex), 
    &(cov->prevIndex), &x_pos_beg, &x_pos_end, &y_pos_beg, &y_pos_end);

    uint64_t pri_len = ug->u.a[pri_uid].len, aux_len = ug->u.a[aux_uid].len;
    pri_uid <<= 1; aux_uid <<= 1; aux_uid += hap->rev;

    // uint32_t i_n = cov->t_ch->k_trans.n, i;

    chain_origin_trans_uid_by_distance(cov, cov->read_g, &pri_uid, 1, x_pos_beg, &pri_len, 
                                                    &aux_uid, 1, y_pos_beg, &aux_len, ug, RC_2, hap->score, __func__);
    
    // fprintf(stderr, "\nocc: %u\n", (uint32_t)(cov->t_ch->k_trans.n - i_n));
    // fprintf(stderr, "#s-utg%.6ul\t%u\t%u\td-utg%.6ul\t%u\t%u\trev(%u)\n", 
    // x->xUid+1, x->x_beg_pos, x->x_end_pos, x->yUid+1, x->y_beg_pos, x->y_end_pos, x->rev);
    // for (i = i_n; i < cov->t_ch->k_trans.n; i++)
    // {
    //     fprintf(stderr, "s-utg%.6ul\t%u\t%u\td-utg%.6ul\t%u\t%u\trev(%u)\n", 
    //     cov->t_ch->k_trans.a[i].qn+1, cov->t_ch->k_trans.a[i].qs, cov->t_ch->k_trans.a[i].qe, 
    //     cov->t_ch->k_trans.a[i].tn+1, cov->t_ch->k_trans.a[i].ts, cov->t_ch->k_trans.a[i].te, 
    //     cov->t_ch->k_trans.a[i].rev);
    // }
}

void collect_purge_trans_cov(ma_ug_t *ug, hap_overlaps_list* ha, hap_cov_t *cov, uint64_t* position_index)
{   
    uint32_t v, i;
    hap_overlaps *x = NULL;
    for (v = 0; v < ha->num; v++)
    {
        for (i = 0; i < ha->x[v].a.n; i++)
        {
            x = &(ha->x[v].a.a[i]);
            if(x->yUid < x->xUid) continue;
            chain_origin_trans_uid_by_purge(x, ug, cov, position_index);
        }
    }
}

void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, 
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, 
uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain, 
uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
{
    p_g_t *pg = NULL;
    asg_t* nsg = ug->g;
    uint32_t v, rId, uId, i, offset;
    ma_utg_t* reads = NULL;
    uint64_t* position_index = NULL;
    if(cov) position_index = cov->pos_idx;
    else position_index = (uint64_t*)malloc(sizeof(uint64_t)*read_g->n_seq);
    memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq);
    
    hap_overlaps_list all_ovlp;
    init_hap_overlaps_list(&all_ovlp, nsg->n_seq);
    hap_overlaps_list back_all_ovlp;
    init_hap_overlaps_list(&back_all_ovlp, nsg->n_seq);
    ///uint32_t junk_cov, hap_cov, dip_cov, junk_occ, repeat_occ, single_cov;
    asg_arc_t t, *p = NULL;
    int r;
    hap_alignment_struct_pip hap_buf;
    long long k_mer_only, coverage_only;

    if(asm_opt.hom_global_coverage != -1)
    {
        hap_buf.cov_threshold = asm_opt.hom_global_coverage;
    }
    else
    {
        hap_buf.cov_threshold = get_read_coverage_thres(ug, read_g, ruIndex, position_index, 
        sources, coverage_cut, read_g->n_seq, COV_COUNT, &k_mer_only, &coverage_only);
    }
    
    
    for (v = 0; v < nsg->n_seq; v++)
    {
        uId = v;
        reads = &(ug->u.a[uId]);
        for (i = 0, offset = 0; i < reads->n; i++)
        {
            rId = reads->a[i]>>33;
            set_R_to_U(ruIndex, rId, uId, 1, &(read_g->seq[rId].c));

            position_index[rId] = offset;
            position_index[rId] = position_index[rId] << 32;
            position_index[rId] = position_index[rId] | (uint64_t)i;

            offset += (uint32_t)reads->a[i];
        }
    }


    init_hap_alignment_struct_pip(&hap_buf, asm_opt.thread_num, nsg->n_seq, ug, read_g,  
    sources, reverse_sources, ruIndex, coverage_cut, position_index, density, max_hang, min_ovlp, 
    0.1, &all_ovlp, cov);
    
    if(hap_buf.cov_threshold < 0)
    {
        if(if_ploid_sample(ug, read_g, ruIndex, sources, reverse_sources, coverage_cut, 
        &hap_buf, &all_ovlp, &back_all_ovlp, purege_minLen, 0.333))
        {
            ///if peak is het, coverage peak is more reliable 
            hap_buf.cov_threshold = coverage_only * HET_PEAK_RATE;
        }
        else
        {
            ///if peak is homo, k-mer peak is more reliable 
            hap_buf.cov_threshold = k_mer_only * HOM_PEAK_RATE;
        }
    }
    if(asm_opt.hom_global_coverage == -1) asm_opt.hom_global_coverage = hap_buf.cov_threshold;
    fprintf(stderr, "[M::%s] purge duplication coverage threshold: %lld\n", __func__, hap_buf.cov_threshold);
    if(just_coverage) goto end_coverage;

    kt_for(asm_opt.thread_num, hap_alignment_advance_worker, &hap_buf, nsg->n_seq);
    
    ///if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp);
    filter_hap_overlaps_by_length(&all_ovlp, purege_minLen);

    normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex);

    if(asm_opt.polyploidy <= 2)
    {
        mc_solve(&all_ovlp, cov->t_ch, NULL, ug, read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL);
        ///pt_solve(&all_ovlp, cov->t_ch, ug, read_g, 0.8, R_INF.trio_flag);
    } 
    
    if(collect_p_trans)
    {
        collect_purge_trans_cov(ug, &all_ovlp, cov, position_index);
        ///goto end_coverage;
    } 
    
    pg = init_p_g_t(ug, cov, read_g);
    
    normalize_hap_overlaps_advance_by_p_g_t(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, 
    ruIndex, pg, cov, 0.8);

    ///normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex);
    ///debug_hap_overlaps(&all_ovlp, &back_all_ovlp);

    remove_contained_haplotig(&all_ovlp, ug, nsg, pg->pg_h_lev, cov);
    
    if(just_contain == 0)
    {
        for (v = 0; v < all_ovlp.num; v++)
        {
            uId = v;
            if(pg->pg_h_lev->seq[uId].del || pg->pg_h_lev->seq[uId].c == ALTER_LABLE) continue;
            for (i = 0; i < all_ovlp.x[uId].a.n; i++)
            {
                if(all_ovlp.x[uId].a.a[i].status == DELETE) continue;
                if(pg->pg_h_lev->seq[all_ovlp.x[uId].a.a[i].xUid].c == ALTER_LABLE||
                   pg->pg_h_lev->seq[all_ovlp.x[uId].a.a[i].xUid].del||
                   pg->pg_h_lev->seq[all_ovlp.x[uId].a.a[i].yUid].c == ALTER_LABLE||
                   pg->pg_h_lev->seq[all_ovlp.x[uId].a.a[i].yUid].del)
                {
                    continue;
                }

                
                ///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
                
                r = get_hap_arch(&(all_ovlp.x[uId].a.a[i]), ug->u.a[all_ovlp.x[uId].a.a[i].xUid].len, 
                ug->u.a[all_ovlp.x[uId].a.a[i].yUid].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
                
                // if(all_ovlp.x[uId].a.a[i].xUid == 118 && all_ovlp.x[uId].a.a[i].yUid == 82)
                // {
                //     fprintf(stderr, "r: %d\n", r);
                //     print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
                // }
                
                if(r < 0) continue;
                p = asg_arc_pushp(pg->pg_h_lev);
                *p = t;
            }
        }

        asg_cleanup(pg->pg_h_lev);
		asg_symm(pg->pg_h_lev);
        ///may need to do transitive reduction
        clean_purge_graph(pg->pg_h_lev, drop_ratio, 1);

        // if(debug_enable) print_purge_gfa(ug, purge_g);
        // if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp);
        /*******************************for debug************************************/
        // print_het_ovlp(pg, ug, &all_ovlp, 0.8);
        /*******************************for debug************************************/

        link_unitigs(pg->pg_h_lev, ug, &all_ovlp, ruIndex, reverse_sources, coverage_cut, read_g, position_index, 
        &(hap_buf.buf[0].u_buffer), &(hap_buf.buf[0].u_buffer_tailIndex), &(hap_buf.buf[0].u_buffer_prevIndex),
        max_hang, min_ovlp, edge, hap_buf.buf[0].visit, cov);
    }

    for (v = 0; v < all_ovlp.num; v++)
    {
        uId = v;
        if(pg->pg_h_lev->seq[uId].c == ALTER_LABLE)
        {
            ug->g->seq[uId].c = ALTER_LABLE;
        }
    }

    end_coverage:
    uint32_t is_Unitig;
    for (v = 0; v < ruIndex->len; v++)
    {
        get_R_to_U(ruIndex, v, &uId, &is_Unitig);
        if(is_Unitig == 1) ruIndex->index[v] = (uint32_t)-1;
    }
    asg_cleanup(nsg);
    destory_hap_overlaps_list(&all_ovlp);
    destory_hap_overlaps_list(&back_all_ovlp);
    if(cov) memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq);
    else free(position_index);
    destory_hap_alignment_struct_pip(&hap_buf);
    destory_p_g_t(&pg);
}