hifiasm/Overlaps.h

#ifndef __OVERLAPS__
#define __OVERLAPS__
#include <stdio.h>
#include <stdint.h>
#include "kvec.h"
#include "kdq.h"

///#define MIN_OVERLAP_LEN 2000
///#define MIN_OVERLAP_LEN 500
///#define MIN_OVERLAP_LEN 50
///#define MIN_OVERLAP_LEN 50
///#define MIN_OVERLAP_COVERAGE 1
///#define MIN_OVERLAP_COVERAGE 0
///#define MAX_HANG_LEN 1000
///#define MAX_HANG_PRE 0.8
///#define GAP_FUZZ 1000
///#define MAX_SHORT_TIPS 3
///#define MAX_BUBBLE_DIST 10000000
#define SMALL_BUBBLE_SIZE (uint32_t)-1
//#define SMALL_BUBBLE_SIZE 1000
#define PRIMARY_LABLE 0
#define ALTER_LABLE 1
#define HAP_LABLE 2
#define FAKE_LABLE 4
#define TRIO_THRES 0.9
#define DOUBLE_CHECK_THRES 0.1
#define FINAL_DOUBLE_CHECK_THRES 0.2
#define CHIMERIC_TRIM_THRES 4
// #define PRIMARY_LABLE 1
// #define ALTER_LABLE 2
// #define HAP_LABLE 4


#define Get_qn(RECORD) ((uint32_t)((RECORD).qns>>32))
#define Get_qs(RECORD) ((uint32_t)((RECORD).qns))
#define Get_qe(RECORD) ((RECORD).qe)
#define Get_tn(RECORD) ((RECORD).tn)
#define Get_ts(RECORD) ((RECORD).ts)
#define Get_te(RECORD) ((RECORD).te)

#define LONG_TIPS 0
#define TWO_INPUT 1
#define TWO_OUTPUT 2
#define MUL_INPUT 3
#define MUL_OUTPUT 4
#define END_TIPS 5
#define LONG_TIPS_UNDER_MAX_EXT 6
#define LOOP 7

#define TRIM 10
#define CUT 11
#define CUT_DIF_HAP 12


///query is the read itself
typedef struct {
	uint64_t qns;
	uint32_t qe, tn, ts, te;
	uint32_t ml:31, rev:1;
	uint32_t bl:31, del:1;
	uint8_t el;
	uint8_t no_l_indel;
} ma_hit_t;

typedef struct {
	ma_hit_t* buffer;
    uint32_t size;
    uint32_t length;
	uint8_t is_fully_corrected;
	uint8_t is_abnormal;
} ma_hit_t_alloc;


void init_ma_hit_t_alloc(ma_hit_t_alloc* x);
void clear_ma_hit_t_alloc(ma_hit_t_alloc* x);
void resize_ma_hit_t_alloc(ma_hit_t_alloc* x, uint32_t size);
void destory_ma_hit_t_alloc(ma_hit_t_alloc* x);
void add_ma_hit_t_alloc(ma_hit_t_alloc* x, ma_hit_t* element);
void ma_hit_sort_tn(ma_hit_t *a, long long n);
void ma_hit_sort_qns(ma_hit_t *a, long long n);

int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, 
char* output_file_name);



typedef struct {
	uint32_t s:31, del:1, e;
	uint8_t c;
} ma_sub_t;

void ma_hit_sub(int min_dp, ma_hit_t_alloc* sources, long long n_read, uint64_t* readLen, 
long long mini_overlap_length, ma_sub_t** coverage_cut);
void ma_hit_cut(ma_hit_t_alloc* sources, long long n_read, uint64_t* readLen, 
long long mini_overlap_length, ma_sub_t** coverage_cut);
void ma_hit_flt(ma_hit_t_alloc* sources, long long n_read, const ma_sub_t *coverage_cut, 
int max_hang, int min_ovlp);
long long get_specific_overlap(ma_hit_t_alloc* x, uint32_t qn, uint32_t tn);


typedef struct {
	uint64_t ul;
	uint32_t v;
	uint32_t ol:31, del:1;
	uint8_t strong;
	uint8_t el;
	uint8_t no_l_indel;
} asg_arc_t;

typedef struct {
	size_t n, m;
	asg_arc_t* a;
} kv_asg_arc_t;


typedef struct {
	uint32_t len:31, circ:1; // len: length of the unitig; circ: circular if non-zero
	uint32_t start, end; // start: starting vertex in the string graph; end: ending vertex
	uint32_t m, n; // number of reads
	uint64_t *a; // list of reads
	char *s; // unitig sequence is not null
} ma_utg_t;



typedef struct {
	uint32_t len:31, del:1;
	uint8_t c;
} asg_seq_t;

typedef struct {
	uint32_t m_arc, n_arc:31, is_srt:1;
	asg_arc_t *arc;
	uint32_t m_seq, n_seq:31, is_symm:1;
	uint32_t r_seq;

	asg_seq_t *seq;
	uint64_t *idx;

	uint8_t* seq_vis;

	uint32_t n_F_seq;
	ma_utg_t* F_seq;
} asg_t;

asg_t *asg_init(void);
void asg_destroy(asg_t *g);
void asg_arc_sort(asg_t *g);
void asg_seq_set(asg_t *g, int sid, int len, int del);
void asg_arc_index(asg_t *g);
void asg_cleanup(asg_t *g);
void asg_symm(asg_t *g);
void print_gfa(asg_t *g);


typedef struct { size_t n, m; uint64_t *a; } asg64_v;


typedef struct { size_t n, m; ma_utg_t *a; } ma_utg_v;

typedef struct {
	ma_utg_v u;
	asg_t *g;
	kvec_t(uint64_t) occ;
} ma_ug_t;

typedef struct {
	uint32_t utg:31, ori:1, start, len;
} utg_intv_t;


#define MA_HT_INT        (-1)
#define MA_HT_QCONT      (-2)
#define MA_HT_TCONT      (-3)
#define MA_HT_SHORT_OVLP (-4)

///in default, max_hang = 1000, int_frac = 0.8, min_ovlp = 50
static inline int ma_hit2arc(const ma_hit_t *h, int ql, int tl, int max_hang, float int_frac, int min_ovlp, asg_arc_t *p)
{
	int32_t tl5, tl3, ext5, ext3, qs = (int32_t)h->qns;
	uint32_t u, v, l; // u: query end; v: target end; l: length from u to v

	///if query and target are in different strand
	if (h->rev) tl5 = tl - h->te, tl3 = h->ts; // tl5: 5'-end overhang (on the query strand); tl3: similar
	else tl5 = h->ts, tl3 = tl - h->te;

	///ext5 and ext3 is the hang on left side and right side, respectively
	ext5 = qs < tl5? qs : tl5;
	ext3 = ql - (int)h->qe < tl3? ql - (int)h->qe : tl3;


	/**
	if (ext5 > max_hang || ext3 > max_hang || h->qe - qs < (h->qe - qs + ext5 + ext3) * int_frac)
		return MA_HT_INT;
	**/
	///ext3 and ext5 should be always 0
	if (ext5 > max_hang || ext3 > max_hang 
	|| h->qe - qs < (h->qe - qs + ext5 + ext3) * int_frac
	|| h->te - h->ts < (h->te - h->ts + ext5 + ext3) * int_frac)
	{
		return MA_HT_INT;
	}
		
	/**
	********************************query-to-target overlap****************************
	case 1: u = 0, rev = 0                           in the view of target: direction is 1 
	query: CCCCCCCCTAATTAAAAT                        target: TAATTAAAATGGGGGG (use ex-target as query)
	               ||||||||||         <--->                  ||||||||||
	       target: TAATTAAAATGGGGGG           query: CCCCCCCCTAATTAAAAT (use ex-query as target)

	case 2: u = 0, rev = 1                           in the view of target: direction is 0
	query: CCCCCCCCTAATTAAAAT					     target: CCCCCCATTTTAATTA  (use ex-target as query)
                   ||||||||||        <--->                         ||||||||||
	       target: TAATTAAAATGGGGGG                         query: ATTTTAATTAGGGGGGGG  (use ex-query as target)
	********************************query-to-target overlap****************************

	********************************target-to-query overlap****************************
	case 3: u = 1, rev = 0                           in the view of target: direction is 0
			 query: AAATAATATCCCCCCGCG                target: GGGCCGGCAAATAATAT (use ex-target as query)
					|||||||||          <--->                          |||||||||
	target: GGGCCGGCAAATAATAT                                  query: AAATAATATCCCCCCGCG (use ex-query as target)

	case 4: u = 1, rev = 1                          in the view of target: direction is 1
	         query: AAATAATATCCCCCCGCG                        target: ATATTATTTGCCGGCCC (use ex-target as query)
                    |||||||||             <--->                       |||||||||
	target: GGGCCGGCAAATAATAT                          query: CGCGGGGGATATTATTT (use ex-query as target)
	********************************target-to-query overlap****************************
	**/

	if (qs <= tl5 && ql - (int)h->qe <= tl3) return MA_HT_QCONT; // query contained in target
	else if (qs >= tl5 && ql - (int)h->qe >= tl3) return MA_HT_TCONT; // target contained in query 
	else if (qs > tl5) u = 0, v = !!h->rev, l = qs - tl5; ///u = 0 means query-to-target overlap, l is the length of node in string graph (not the overlap length)
	else u = 1, v = !h->rev, l = (ql - h->qe) - tl3; ///u = 1 means target-to-query overlaps, l is the length of node in string graph (not the overlap length)
	if ((int)h->qe - qs + ext5 + ext3 < min_ovlp || (int)h->te - (int)h->ts + ext5 + ext3 < min_ovlp) return MA_HT_SHORT_OVLP; // short overlap
	///u = 0 / 1 means query-to-target / target-to-query overlaps, 
	///l is the length of node in string graph (not the overlap length between two reads)
	u |= h->qns>>32<<1, v |= h->tn<<1;
	/**
	p->ul: |____________31__________|__________1___________|______________32_____________|
	                    qn            direction of overlap       length of this node (not overlap length)
						                (in the view of query)
	p->v : |___________31___________|__________1___________|
				        tn             reverse direction of overlap
						              (in the view of target)
	p->ol: overlap length
	**/
	p->ul = (uint64_t)u<<32 | l, p->v = v, p->ol = ql - l, p->del = 0;
	///l is the length of node in string graph (not the overlap length)

	p->strong = h->ml;
	p->el = h->el;
	p->no_l_indel = h->no_l_indel;
	return l;
}



#define asg_arc_len(arc) ((uint32_t)(arc).ul)
#define asg_arc_n(g, v) ((uint32_t)(g)->idx[(v)])
#define asg_arc_a(g, v) (&(g)->arc[(g)->idx[(v)]>>32])

static inline uint32_t asg_get_arc(asg_t *g, uint32_t v, uint32_t w, asg_arc_t* t)
{
	uint32_t i, nv = asg_arc_n(g, v);
	asg_arc_t *av = asg_arc_a(g, v);
	for (i = 0; i < nv; ++i)
	{
		if(av[i].del) continue;
		if(av[i].v == w)
		{
			(*t) = av[i];
			return 1;
		}
	}

	return 0;
}

// append an arc
static inline asg_arc_t *asg_arc_pushp(asg_t *g)
{
	if (g->n_arc == g->m_arc) {
		g->m_arc = g->m_arc? g->m_arc<<1 : 16;
		g->arc = (asg_arc_t*)realloc(g->arc, g->m_arc * sizeof(asg_arc_t));
	}
	return &g->arc[g->n_arc++];
}

// set asg_arc_t::del for v->w
static inline void asg_arc_del(asg_t *g, uint32_t v, uint32_t w, int del)
{
	uint32_t i, nv = asg_arc_n(g, v);
	asg_arc_t *av = asg_arc_a(g, v);
	for (i = 0; i < nv; ++i)
		if (av[i].v == w) av[i].del = !!del;
}

// set asg_arc_t::del and asg_seq_t::del to 1 for sequence s and all its associated arcs
static inline void asg_seq_del(asg_t *g, uint32_t s)
{
	uint32_t k;
	g->seq[s].del = 1;
	for (k = 0; k < 2; ++k) {
		uint32_t i, v = s<<1 | k;
		uint32_t nv = asg_arc_n(g, v);
		asg_arc_t *av = asg_arc_a(g, v);
		for (i = 0; i < nv; ++i) {
			av[i].del = 1;
			asg_arc_del(g, av[i].v^1, v^1, 1);
		}
	}
}

static inline void asg_seq_drop(asg_t *g, uint32_t s)
{
	///s is not at primary
	if(g->seq[s].c == ALTER_LABLE)
	{
		uint32_t k;
		for (k = 0; k < 2; ++k) 
		{
			///two directions of this node
			uint32_t i, v = s<<1 | k;
			uint32_t nv = asg_arc_n(g, v);
			asg_arc_t *av = asg_arc_a(g, v);
			for (i = 0; i < nv; ++i) 
			{
				if(av[i].del) continue;
				///if output node is at primary
				/****************************may have hap bugs********************************/
				///if(g->seq[(av[i].v>>1)].c == PRIMARY_LABLE)
				///if(g->seq[(av[i].v>>1)].c == PRIMARY_LABLE || g->seq[(av[i].v>>1)].c == HAP_LABLE)
				if(g->seq[(av[i].v>>1)].c != ALTER_LABLE)
				{/****************************may have hap bugs********************************/
					av[i].del = 1;
					asg_arc_del(g, av[i].v^1, v^1, 1);
				}
			}
		}
	}
}





/******************
 * Bubble popping *
 ******************/

typedef struct {
	uint32_t p; // the optimal parent vertex
	uint32_t d; // the shortest distance from the initial vertex
	uint32_t c; // max count of positive reads
	uint32_t m; // max count of negative reads
	uint32_t np; // max count of non-positive reads
	uint32_t nc; // max count of reads, no matter positive or negative
	uint32_t r:31, s:1; // r: the number of remaining incoming arc; s: state
	//s: state, s=0, this edge has not been visited, otherwise, s=1
} binfo_t;

typedef struct {
	///all information for each node
	binfo_t *a;
	kvec_t(uint32_t) S; // set of vertices without parents, nodes with all incoming edges visited
	kvec_t(uint32_t) T; // set of tips
	kvec_t(uint32_t) b; // visited vertices
	kvec_t(uint32_t) e; // visited edges/arcs
} buf_t;


typedef struct {
	kvec_t(uint64_t) Nodes; 
	kvec_t(uint64_t) Edges; 
	uint32_t pre_n_seq, seqID; 
} C_graph;

typedef struct {
	kvec_t(uint8_t) a;
	uint32_t i;
} kvec_t_u8_warp;

typedef struct {
	kvec_t(uint32_t) a;
	uint32_t i;
} kvec_t_u32_warp;

typedef struct {
	kvec_t(int32_t) a;
	uint32_t i;
} kvec_t_i32_warp;

typedef struct {
	kvec_t(uint64_t) a;
	uint64_t i;
} kvec_t_u64_warp;

typedef struct {
	kvec_t(asg_arc_t) a;
	uint64_t i;
}kvec_asg_arc_t_warp;

void sort_kvec_t_u64_warp(kvec_t_u64_warp* u_vecs, uint32_t is_descend);
int asg_arc_del_multi(asg_t *g);
int asg_arc_del_asymm(asg_t *g);

typedef struct {
	uint32_t q_pos;
	uint32_t t_pos;
	uint32_t t_id;
	uint32_t is_color;
} Hap_Align;

typedef struct {
	kvec_t(Hap_Align) x;
	uint64_t i;
} Hap_Align_warp;

typedef struct {
	buf_t* b_0;
    uint32_t untigI;
    uint32_t readI;
    uint32_t offset;
} rIdContig;

// count the number of outgoing arcs, including reduced arcs
static inline int count_out_with_del(const asg_t *g, uint32_t v)
{
	uint32_t nv = asg_arc_n(g, v);
	return nv;
}


// count the number of outgoing arcs, including reduced arcs
static inline int count_out_without_del(const asg_t *g, uint32_t v)
{
	uint32_t i, n, nv = asg_arc_n(g, v);
	const asg_arc_t *av = asg_arc_a(g, v);

	for (i = n = 0; i < nv; ++i)
		if (!av[i].del) ++n;
	return n;
}


void build_string_graph_without_clean(
int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
long long n_read, uint64_t* readLen, long long mini_overlap_length, 
long long max_hang_length, long long clean_round, long long gap_fuzz,
float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, 
long long bubble_dist, int read_graph, int write);

void debug_info_of_specfic_read(char* name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int id, char* command);
void collect_abnormal_edges(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long readNum);
void add_overlaps(ma_hit_t_alloc* source_paf, ma_hit_t_alloc* dest_paf, uint64_t* source_index, long long listLen);
void remove_overlaps(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen);
void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf, 
uint64_t* source_index, long long listLen);

#define EvaluateLen(U, id) ((U).a[(id)].start)
#define IsMerge(U, id) ((U).a[(id)].end)
#define kv_reuse(v, rn, rm, r) ((v).n = (rn), (v).m = (rm), (v).a = (r))
#define Get_vis(visit, v, d) (((visit)[(v)>>1])&(((((v)<<(d))&1)+1)))
#define Set_vis(visit, v, d) (((visit)[(v)>>1])|=(((((v)<<(d))&1)+1)))





typedef struct {
	uint64_t len;
	uint32_t* index;
} R_to_U;

void init_R_to_U(R_to_U* x, uint64_t len);
void destory_R_to_U(R_to_U* x);
void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig, uint8_t* flag);
void get_R_to_U(R_to_U* x, uint32_t rID, uint32_t* uID, uint32_t* is_Unitig);
void transfor_R_to_U(R_to_U* x);
void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, kvec_asg_arc_t_warp* edge, int require_equal_nv, int test_tangle);
long long asg_arc_del_simple_circle_untig(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *g, long long circleLen, int is_drop);

typedef struct {
	asg_t* g;

	asg_arc_t *av;
	uint32_t nv;
	uint32_t av_i;

	asg_arc_t* new_edges;
	uint32_t new_edges_n;
	uint32_t new_edges_i;
} Edge_iter;

typedef struct {
	asg_arc_t x;
    uint64_t Off;
    uint64_t weight;
}asg_arc_t_offset;

typedef struct {
	kvec_t(asg_arc_t_offset) a;
	uint64_t i;
}kvec_asg_arc_t_offset;



void init_Edge_iter(asg_t* g, uint32_t v, asg_arc_t* new_edges, uint32_t new_edges_n, Edge_iter* x);
int get_arc_t(Edge_iter* x, asg_arc_t* get);


inline int get_real_length(asg_t *g, uint32_t v, uint32_t* v_s)
{
    uint32_t i, kv = 0;
    for (i = 0, kv = 0; i < asg_arc_n(g, v); i++)
    {
        if(!asg_arc_a(g, v)[i].del)
        {
            if(v_s) v_s[kv] = asg_arc_a(g, v)[i].v;
            kv++;
        }
    }

    return kv;
}

inline uint32_t check_tip(asg_t *sg, uint32_t begNode, uint32_t* endNode, buf_t* b, uint32_t max_ext)
{
    ///cut tip of length <= max_ext
    uint32_t v = begNode, w;
    uint32_t kv;
    uint32_t eLen = 0;
    (*endNode) = (uint32_t)-1;
    b->b.n = 0;
    while (1)
    {
        kv = get_real_length(sg, v, NULL);
        (*endNode) = v;
        eLen++;
        if(b) kv_push(uint32_t, b->b, v);
        if(kv == 0) return END_TIPS;
        if(kv > 1) return MUL_OUTPUT;
        ///if(eLen > max_ext) return LONG_TIPS;
        ///kv must be 1 here
        kv = get_real_length(sg, v, &w);
        ///here this value must be >= 1
        if(get_real_length(sg, w^1, NULL)!=1) return MUL_INPUT;
        v = w;
        if(v == begNode) return LOOP;
        if(eLen >= max_ext) return LONG_TIPS;
    }
}

inline uint32_t get_unitig(asg_t *sg, ma_ug_t *ug, uint32_t begNode, uint32_t* endNode, 
long long* nodeLen, long long* baseLen, long long* max_stop_nodeLen, long long* max_stop_baseLen, 
uint32_t stops_threshold, buf_t* b)
{
	ma_utg_v* u = NULL; 
    uint32_t v = begNode, w, k;
    uint32_t kv, return_flag, n_stops = 0;
	long long pre_baseLen = 0, pre_nodeLen = 0;
	long long cur_baseLen = 0, cur_nodeLen = 0;
	(*max_stop_nodeLen) = (*max_stop_baseLen) = (*nodeLen) = (*baseLen) = 0;
    (*endNode) = (uint32_t)-1;
	if(ug!=NULL) u = &(ug->u);

    while (1)
    {
        kv = get_real_length(sg, v, NULL);
        (*endNode) = v;
		if(u == NULL) 
		{
			(*nodeLen)++;
		}
		else
		{
			(*nodeLen) += EvaluateLen((*u), v>>1);
		}
        if(b) kv_push(uint32_t, b->b, v);
		///means reach the end of a unitig
		if(kv!=1) (*baseLen) += sg->seq[v>>1].len;
		if(kv==0)
		{
			return_flag = END_TIPS;
			break;
			///return END_TIPS;
		} 
		if(kv>1)
		{
			return_flag = MUL_OUTPUT;
			break;
			///return MUL_OUTPUT;
		}
        ///kv must be 1 here
        kv = get_real_length(sg, v, &w);
		///means reach the end of a unitig
        if(get_real_length(sg, w^1, NULL)!=1)
		{

			n_stops++;
			if(n_stops >= stops_threshold)
			{
				(*baseLen) += sg->seq[v>>1].len;
				return_flag = MUL_INPUT;
				break;
				///return MUL_INPUT;
			}
			else
			{
				for (k = 0; k < asg_arc_n(sg, v); k++)
				{
					if(asg_arc_a(sg, v)[k].del) continue;
					///here is just one undeleted edge
					(*baseLen) += asg_arc_len(asg_arc_a(sg, v)[k]);
					break;
				}
			}

			cur_baseLen = (*baseLen) - pre_baseLen;
			pre_baseLen = (*baseLen);
			if(cur_baseLen > (*max_stop_baseLen))
			{
				(*max_stop_baseLen) = cur_baseLen;
			}


			cur_nodeLen = (*nodeLen) - pre_nodeLen;
			pre_nodeLen = (*nodeLen);
			if(cur_nodeLen > (*max_stop_nodeLen))
			{
				(*max_stop_nodeLen) = cur_nodeLen;
			}
		}
		else
		{
			for (k = 0; k < asg_arc_n(sg, v); k++)
			{
				if(asg_arc_a(sg, v)[k].del) continue;
				///here is just one undeleted edge
				(*baseLen) += asg_arc_len(asg_arc_a(sg, v)[k]);
				break;
			}
		}
		

        v = w;
        if(v == begNode)
		{
			return_flag = LOOP;
			break;
			///return LOOP;
		} 
    }




	cur_baseLen = (*baseLen) - pre_baseLen;
	pre_baseLen = (*baseLen);
	if(cur_baseLen > (*max_stop_baseLen))
	{
		(*max_stop_baseLen) = cur_baseLen;
	}


	cur_nodeLen = (*nodeLen) - pre_nodeLen;
	pre_nodeLen = (*nodeLen);
	if(cur_nodeLen > (*max_stop_nodeLen))
	{
		(*max_stop_nodeLen) = cur_nodeLen;
	}

	return return_flag;
}

#define UNAVAILABLE (uint32_t)-1
#define PLOID 0
#define NON_PLOID 1
#define DIFF_HAP_RATE 0.75
#define TRIO_DROP_THRES 0.9
#define TRIO_DROP_LENGTH_THRES 0.8
#define MAX_STOP_RATE 0.6
#define TANGLE_MISSED_THRES 0.6
///if ug == NULL, nsg should be equal to read_sg
inline uint32_t check_different_haps(asg_t *nsg, ma_ug_t *ug, asg_t *read_sg, 
uint32_t v_0, uint32_t v_1, ma_hit_t_alloc* reverse_sources, buf_t* b_0, buf_t* b_1, 
R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
{
	uint32_t vEnd, qn, tn, j, is_Unitig, uId;
	long long ELen_0, ELen_1, tmp, max_stop_nodeLen, max_stop_baseLen;

	b_0->b.n = b_1->b.n = 0;
	if(get_unitig(nsg, ug, v_0, &vEnd, &ELen_0, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 
	stops_threshold, b_0) == LOOP)
	{
		return UNAVAILABLE;
	}
	if(get_unitig(nsg, ug, v_1, &vEnd, &ELen_1, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 
	stops_threshold, b_1) == LOOP)
	{
		return UNAVAILABLE;
	}
	if(ELen_0<=min_edge_length || ELen_1<=min_edge_length) return UNAVAILABLE;

	rIdContig b_max, b_min;
	b_max.b_0 = b_min.b_0 = NULL;
	b_max.offset = b_max.readI = b_max.untigI = 0;
	b_min.offset = b_min.readI = b_min.untigI = 0;

	if(ELen_0<=ELen_1)
	{
        b_min.b_0 = b_0;
        b_max.b_0 = b_1;
    }
    else
    {
        b_min.b_0 = b_1;
        b_max.b_0 = b_0;
    }

	uint32_t max_count = 0, min_count = 0;
	ma_utg_t *node_min = NULL, *node_max = NULL;
	if(ug != NULL)
	{
		/*****************************label all unitigs****************************************/
		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
		{
			node_max = &(ug->u.a[b_max.b_0->b.a[b_max.untigI]>>1]);
			///each read
			for (b_max.readI = 0; b_max.readI < node_max->n; b_max.readI++)
			{
				qn = (node_max->a[b_max.readI]>>33);
				set_R_to_U(ruIndex, qn, (b_max.b_0->b.a[b_max.untigI]>>1), 1, &(read_sg->seq[qn].c));
			}
		}
		/*****************************label all unitigs****************************************/

		///each unitig
		for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
		{
			
			node_min = &(ug->u.a[(b_min.b_0->b.a[b_min.untigI]>>1)]);

			///each read
			for (b_min.readI = 0; b_min.readI < node_min->n; b_min.readI++)
        	{
				qn = node_min->a[b_min.readI]>>33;

				/************************BUG: don't forget****************************/
				if(reverse_sources[qn].length > 0) min_count++;
				///if(reverse_sources[qn].length >= 0) min_count++;
				/************************BUG: don't forget****************************/
				for (j = 0; j < (long long)reverse_sources[qn].length; j++)
				{
					tn = Get_tn(reverse_sources[qn].buffer[j]);
					if(read_sg->seq[tn].del == 1)
					{
						get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
						if(tn == (uint32_t)-1 || is_Unitig == 1 || read_sg->seq[tn].del == 1) continue;
					}

					get_R_to_U(ruIndex, tn, &uId, &is_Unitig);
					if(uId!=(uint32_t)-1 && is_Unitig == 1)
					{
						max_count++;
						break;
					}
				}
			}
		}
		/*****************************label all unitigs****************************************/
		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
		{
			node_max = &(ug->u.a[b_max.b_0->b.a[b_max.untigI]>>1]);
			///each read
			for (b_max.readI = 0; b_max.readI < node_max->n; b_max.readI++)
			{
				qn = (node_max->a[b_max.readI]>>33);
				ruIndex->index[qn] = (uint32_t)-1;
			}
		}
		/*****************************label all unitigs****************************************/
	}
	else
	{
		/*****************************label all reads****************************************/
		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
		{
			qn = (b_max.b_0->b.a[b_max.untigI]>>1);
			set_R_to_U(ruIndex, qn, 1, 1, &(read_sg->seq[qn].c));
		}
		/*****************************label all reads****************************************/

		///each read
		for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
		{
			qn = (b_min.b_0->b.a[b_min.untigI]>>1);

			/************************BUG: don't forget****************************/
			if(reverse_sources[qn].length > 0) min_count++;
			///if(reverse_sources[qn].length >= 0) min_count++;
			/************************BUG: don't forget****************************/

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


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

		/*****************************label all reads****************************************/
		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
		{
			qn = (b_max.b_0->b.a[b_max.untigI]>>1);
			ruIndex->index[qn] = (uint32_t)-1;
		}
		/*****************************label all reads****************************************/
	}

	// if(((v_0==7707) && (v_1==26867))||((v_1==7707) && (v_0==26867)))
	// {
	// 	fprintf(stderr, "******\nv_0>>1: %u, v_0&1: %u, ELen_0: %u\n", v_0>>1, v_0&1, (uint32_t)ELen_0);
	// 	fprintf(stderr, "v_1>>1: %u, v_1&1: %u, ELen_1: %u\n", v_1>>1, v_1&1, (uint32_t)ELen_1);
	// 	fprintf(stderr, "min_count: %u, max_count: %u, DIFF_HAP_RATE: %f\n\n", 
	// 	min_count, max_count, DIFF_HAP_RATE);
	// }
	
	if(min_count == 0) return UNAVAILABLE;
	if(max_count > min_count*DIFF_HAP_RATE) return PLOID;
	return NON_PLOID;
}

inline uint32_t check_different_haps_naive(asg_t *nsg, ma_ug_t *ug, asg_t *read_sg, 
uint32_t v_0, uint32_t v_1, ma_hit_t_alloc* reverse_sources, buf_t* b_0, buf_t* b_1, 
R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
{
	uint32_t vEnd, qn, tn, j, is_Unitig;
	long long ELen_0, ELen_1, tmp, max_stop_nodeLen, max_stop_baseLen;

	b_0->b.n = b_1->b.n = 0;
	if(get_unitig(nsg, ug, v_0, &vEnd, &ELen_0, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 
	stops_threshold, b_0) == LOOP)
	{
		return UNAVAILABLE;
	} 
	if(get_unitig(nsg, ug, v_1, &vEnd, &ELen_1, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 
	stops_threshold, b_1) == LOOP)
	{
		return UNAVAILABLE;
	}

	if(ELen_0<=min_edge_length || ELen_1<=min_edge_length) return UNAVAILABLE;

	rIdContig b_max, b_min;
	b_max.b_0 = b_min.b_0 = NULL;
	b_max.offset = b_max.readI = b_max.untigI = 0;
	b_min.offset = b_min.readI = b_min.untigI = 0;

	if(ELen_0<=ELen_1)
	{
        b_min.b_0 = b_0;
        b_max.b_0 = b_1;
    }
    else
    {
        b_min.b_0 = b_1;
        b_max.b_0 = b_0;
    }

	uint32_t max_count = 0, min_count = 0;
	ma_utg_t *node_min = NULL, *node_max = NULL;

	if(ug != NULL)
	{
		///each unitig
		for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
		{
			
			node_min = &(ug->u.a[(b_min.b_0->b.a[b_min.untigI]>>1)]);

			///each read
			for (b_min.readI = 0; b_min.readI < node_min->n; b_min.readI++)
        	{
				qn = node_min->a[b_min.readI]>>33;

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

					///each unitig
					for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
					{
						node_max = &(ug->u.a[b_max.b_0->b.a[b_max.untigI]>>1]);
						///each read
						for (b_max.readI = 0; b_max.readI < node_max->n; b_max.readI++)
						{
							if(tn == (node_max->a[b_max.readI]>>33))
							{
								max_count++;
								goto end_check_different_haps_ug;
							}
						}
					}
				}

				end_check_different_haps_ug:;
			}
		}
	}
	else
	{
		///each read
		for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
		{
			qn = (b_min.b_0->b.a[b_min.untigI]>>1);

			if(reverse_sources[qn].length > 0) min_count++;

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

				///each read
				for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
                {
					if((b_max.b_0->b.a[b_max.untigI]>>1) == tn)
					{
						max_count++;
						goto end_check_different_haps_non_ug;
					}
				}
			}

			end_check_different_haps_non_ug:;
		}
	}

	// if(((v_0==7707) && (v_1==26867))||((v_1==7707) && (v_0==26867)))
	// {
	// 	fprintf(stderr, "******\nv_0>>1: %u, v_0&1: %u, ELen_0: %u\n", v_0>>1, v_0&1, (uint32_t)ELen_0);
	// 	fprintf(stderr, "v_1>>1: %u, v_1&1: %u, ELen_1: %u\n", v_1>>1, v_1&1, (uint32_t)ELen_1);
	// 	fprintf(stderr, "min_count: %u, max_count: %u, DIFF_HAP_RATE: %f\n\n", 
	// 	min_count, max_count, DIFF_HAP_RATE);
	// }
	
	if(min_count == 0) return UNAVAILABLE;
	if(max_count > min_count*DIFF_HAP_RATE) return PLOID;
	return NON_PLOID;
}



typedef struct {
	uint32_t father_occ;
	uint32_t mother_occ;
	uint32_t ambig_occ;
	uint32_t drop_occ;
	uint32_t total;
} Trio_counter;

typedef struct {
	uint32_t p; // the optimal parent vertex
	uint32_t d; // the shortest distance from the initial vertex
	uint32_t r:31, s:1; // r: the number of remaining incoming arc; s: state
} binfo_s_t;

typedef struct {
	///all information for each node
	binfo_s_t *a;
	kvec_t(uint32_t) S; // set of vertices without parents, nodes with all incoming edges visited
	kvec_t(uint32_t) b; // visited vertices
	kvec_t(uint32_t) e; // visited edges/arcs
} buf_s_t;

typedef struct{
    buf_s_t *b;
    uint32_t n_thres, n_reads;
    asg_t *g;
    uint32_t check_cross;
	uint64_t bub_dist;
} bub_label_t;

void resolve_tangles(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig, 
long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint32_t trio_flag, 
float drop_ratio);
void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
void rescue_contained_reads_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t chainLenThres, uint32_t is_bubble_check, 
uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp* new_rtg_nodes, bub_label_t* b_mask_t);
void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t is_bubble_check, uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t);
void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, 
ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate);
void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
void rescue_wrong_overlaps_to_unitigs(ma_ug_t *i_ug, asg_t *r_g,  ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges, bub_label_t* b_mask_t);
void get_unitig_trio_flag(ma_utg_t* nsu, uint32_t flag, uint32_t* require, uint32_t* non_require, uint32_t* ambigious);
void rescue_missing_overlaps_backward(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t backward_steps, uint32_t is_bubble_check, uint32_t is_primary_check, bub_label_t* b_mask_t);
uint32_t get_edge_from_source(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, 
R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t query, uint32_t target, asg_arc_t* t);
int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio);
void asg_bub_backtrack_primary(asg_t *g, uint32_t v0, buf_t *b);

typedef struct{
    double weight;
    uint32_t uID:31, del:1;
    uint64_t dis;
	uint64_t occ;
	///uint64_t occ:63, scaff:1;
    ///uint32_t enzyme;
} hc_edge;

typedef struct{
    kvec_t(hc_edge) e;
    kvec_t(hc_edge) f;//forbiden
} hc_linkeage;

typedef struct{
	uint64_t beg, end;
}bed_interval;

typedef struct{
	size_t n, m;
	bed_interval* a;
}bed_in;

typedef struct{
    kvec_t(hc_linkeage) a;
    kvec_t(uint64_t) enzymes;
	kvec_t(bed_in) bed;
    ///uint32_t* u_idx;
	///uint64_t r_num;
} hc_links;


typedef struct{
	kvec_t(uint32_t) uIDs;
	kvec_t(uint32_t) iDXs;
	kvec_t(uint32_t) rescue_hom;
	uint32_t* u_idx;
	uint32_t r_num;
	uint32_t chain_num;
	uint32_t l0_chain, l1_chain;
}trans_chain;

typedef struct {
    uint32_t n;
    uint32_t* cov;
    uint64_t* pos_idx;
    ma_hit_t_alloc* reverse_sources;
    ma_sub_t *coverage_cut;
    R_to_U* ruIndex;
	asg_t *read_g;
    int max_hang;
    int min_ovlp;
    kvec_asg_arc_t_offset u_buffer; 
    kvec_t_i32_warp tailIndex;
    kvec_t_i32_warp prevIndex;
	///hc_links* link;
	trans_chain* t_ch;
}hap_cov_t;

typedef struct{
    ///kvec_t(hc_edge) a;
    size_t n, m; 
    hc_edge *a;
}hc_edge_warp;

void init_hc_links(hc_links* link, uint64_t ug_num, uint64_t r_num);
void destory_hc_links(hc_links* link);
uint64_t get_bub_pop_max_dist(asg_t *g, buf_t *b);
uint64_t get_bub_pop_max_dist_advance(asg_t *g, buf_t *b);
int asg_pop_bubble_primary_trio(ma_ug_t *ug, uint64_t* i_max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov);
uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, uint64_t max_dist, buf_t *b, uint32_t positive_flag, 
uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes, hap_cov_t *cov);

void adjust_utg_by_primary(ma_ug_t **ug, asg_t* read_g, float drop_rate,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, 
long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp,
kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link, bub_label_t* b_mask_t);
ma_ug_t* copy_untig_graph(ma_ug_t *src);
ma_ug_t* output_trio_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, 
uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, 
float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench, bub_label_t* b_mask_t);
asg_t* copy_read_graph(asg_t *src);
ma_ug_t *ma_ug_gen(asg_t *g);
void ma_ug_destroy(ma_ug_t *ug);

inline int inter_interval(int a_s, int a_e, int b_s, int b_e, int* i_s, int* i_e)
{
    if(a_s > b_e || b_s > a_e) return 0;
    if(i_s) (*i_s) = a_s >= b_s? a_s : b_s; ///MAX(a_s, b_s);
    if(i_e) (*i_e) = a_e <= b_e? a_e : b_e; ///MIN(a_e, b_e);
    return 1;
}




#define JUNK_COV 5
#define DISCARD_RATE 0.8

#endif